Search in sources :

Example 1 with Angle

use of javax.measure.quantity.Angle in project sis by apache.

the class FranceGeocentricInterpolationTest method testGridAsFloats.

/**
 * Tests a small grid file with interpolations in geocentric coordinates as {@code float} values.
 *
 * <p>This method is part of a chain.
 * The next method is {@link #testGridAsShorts(DatumShiftGridFile)}.</p>
 *
 * @return the loaded grid with values as {@code float}.
 * @throws URISyntaxException if the URL to the test file can not be converted to a path.
 * @throws IOException if an error occurred while loading the grid.
 * @throws FactoryException if an error occurred while computing the grid.
 * @throws TransformException if an error occurred while computing the envelope.
 */
@TestStep
private static DatumShiftGridFile<Angle, Length> testGridAsFloats() throws URISyntaxException, IOException, FactoryException, TransformException {
    final Path file = getResource(TEST_FILE);
    final DatumShiftGridFile.Float<Angle, Length> grid;
    try (BufferedReader in = Files.newBufferedReader(file)) {
        grid = FranceGeocentricInterpolation.load(in, file);
    }
    assertEquals("cellPrecision", 0.005, grid.getCellPrecision(), STRICT);
    assertEquals("getCellMean", 168.2587, grid.getCellMean(0), 0.0001);
    assertEquals("getCellMean", 58.7163, grid.getCellMean(1), 0.0001);
    assertEquals("getCellMean", -320.1801, grid.getCellMean(2), 0.0001);
    verifyGrid(grid);
    return grid;
}
Also used : Path(java.nio.file.Path) Angle(javax.measure.quantity.Angle) Length(javax.measure.quantity.Length) BufferedReader(java.io.BufferedReader) TestStep(org.apache.sis.test.TestStep)

Example 2 with Angle

use of javax.measure.quantity.Angle in project sis by apache.

the class NADCONTest method testNADCON.

/**
 * Implementation of {@link #testLoader()} and {@link #testNADCON(Path, Path)}.
 *
 * @param  xmin  westmost longitude.
 * @param  xmax  eastmost longitude.
 * @param  ymin  southmost latitude.
 * @param  ymax  northmost latitude.
 */
private static void testNADCON(final Path latitudeShifts, final Path longitudeShifts, final double xmin, final double xmax, final double ymin, final double ymax) throws IOException, FactoryException, TransformException {
    final DatumShiftGridFile<Angle, Angle> grid = NADCON.getOrLoad(latitudeShifts, longitudeShifts);
    assertInstanceOf("Should not be compressed.", DatumShiftGridFile.Float.class, grid);
    assertEquals("coordinateUnit", Units.DEGREE, grid.getCoordinateUnit());
    assertEquals("translationUnit", Units.DEGREE, grid.getTranslationUnit());
    assertEquals("translationDimensions", 2, grid.getTranslationDimensions());
    assertTrue("isCellValueRatio", grid.isCellValueRatio());
    assertTrue("cellPrecision", grid.getCellPrecision() > 0);
    /*
         * Verify the envelope and the conversion between geographic coordinates and grid indices.
         * The cells are expected to have the same size (0.25°) in longitudes and latitudes.
         */
    final double cellSize = 0.25;
    final Envelope envelope = grid.getDomainOfValidity();
    assertEquals("xmin", xmin - cellSize / 2, envelope.getMinimum(0), STRICT);
    assertEquals("xmax", xmax + cellSize / 2, envelope.getMaximum(0), STRICT);
    assertEquals("ymin", ymin - cellSize / 2, envelope.getMinimum(1), STRICT);
    assertEquals("ymax", ymax + cellSize / 2, envelope.getMaximum(1), STRICT);
    assertMatrixEquals("coordinateToGrid", new Matrix3(cellSize, 0, xmin, 0, cellSize, ymin, 0, 0, 1), grid.getCoordinateToGrid().inverse().getMatrix(), STRICT);
    /*
         * Test the Meades Ranch station. If we were using the complete Conus files, we would obtain
         * after conversion the grid indices listed on the left side. But since we are using a sub-set
         * of the grid, we rather obtain the grid indices on the right side.
         *
         *   gridX ≈ 129.83277 in official file,   ≈ 4.83277 in the test file.
         *   gridY ≈  76.89632 in official file,   ≈ 6.89632 in the test file.
         */
    final double[] position = samplePoint(1);
    final double[] expected = samplePoint(2);
    final double[] indices = new double[position.length];
    final double[] vector = new double[2];
    grid.getCoordinateToGrid().transform(position, 0, indices, 0, 1);
    grid.interpolateInCell(indices[0], indices[1], vector);
    vector[0] *= cellSize * DatumShiftGridLoader.DEGREES_TO_SECONDS;
    vector[1] *= cellSize * DatumShiftGridLoader.DEGREES_TO_SECONDS;
    assertArrayEquals("interpolateInCell", expected, vector, 0.5E-5);
    // Same test than above, but let DatumShiftGrid do the conversions for us.
    expected[0] /= DatumShiftGridLoader.DEGREES_TO_SECONDS;
    expected[1] /= DatumShiftGridLoader.DEGREES_TO_SECONDS;
    assertArrayEquals("interpolateAt", expected, grid.interpolateAt(position), ANGULAR_TOLERANCE);
    assertSame("Grid should be cached.", grid, NADCON.getOrLoad(latitudeShifts, longitudeShifts));
}
Also used : Angle(javax.measure.quantity.Angle) Envelope(org.opengis.geometry.Envelope) Matrix3(org.apache.sis.referencing.operation.matrix.Matrix3)

Example 3 with Angle

use of javax.measure.quantity.Angle in project sis by apache.

the class CRSBuilder method verify.

/**
 * Verifies if the user-defined conversion created from GeoTIFF values
 * matches the given conversion created from the EPSG geodetic dataset.
 * This method does not verify the EPSG code of the given conversion.
 *
 * @param  projection  the conversion created from the EPSG geodetic dataset.
 */
private void verify(final Conversion projection, final Unit<Angle> angularUnit, final Unit<Length> linearUnit) throws FactoryException {
    final Unit<Angle> azimuthUnit = createUnit(GeoKeys.AzimuthUnits, (short) 0, Angle.class, Units.DEGREE);
    final String type = getAsString(GeoKeys.CoordTrans);
    if (type != null) {
        /*
             * Compare the name of the map projection declared in the GeoTIFF file with the name
             * of the projection used by the EPSG geodetic dataset.
             */
        final OperationMethod method = projection.getMethod();
        if (!IdentifiedObjects.isHeuristicMatchForName(method, type)) {
            Identifier expected = IdentifiedObjects.getIdentifier(method, Citations.GEOTIFF);
            if (expected == null) {
                expected = IdentifiedObjects.getIdentifier(method, null);
            }
            warning(Resources.Keys.NotTheEpsgValue_5, IdentifiedObjects.getIdentifierOrName(projection), expected.getCode(), GeoKeys.name(GeoKeys.CoordTrans), type, "");
        }
        /*
             * Compare the parameter values with the ones declared in the EPSG geodetic dataset.
             */
        final ParameterValueGroup parameters = projection.getParameterValues();
        for (final short key : remainingKeys()) {
            final Unit<?> unit;
            switch(GeoKeys.unitOf(key)) {
                case GeoKeys.RATIO:
                    unit = Units.UNITY;
                    break;
                case GeoKeys.LINEAR:
                    unit = linearUnit;
                    break;
                case GeoKeys.ANGULAR:
                    unit = angularUnit;
                    break;
                case GeoKeys.AZIMUTH:
                    unit = azimuthUnit;
                    break;
                default:
                    continue;
            }
            try {
                verify(projection, parameters.parameter("GeoTIFF:" + key).doubleValue(unit), key, unit);
            } catch (ParameterNotFoundException e) {
                warning(Resources.Keys.UnexpectedParameter_2, type, GeoKeys.name(key));
            }
        }
    }
}
Also used : Identifier(org.opengis.metadata.Identifier) Angle(javax.measure.quantity.Angle) ParameterValueGroup(org.opengis.parameter.ParameterValueGroup) ParameterNotFoundException(org.opengis.parameter.ParameterNotFoundException) OperationMethod(org.opengis.referencing.operation.OperationMethod)

Example 4 with Angle

use of javax.measure.quantity.Angle in project sis by apache.

the class DefaultGeodeticCRS method formatTo.

/**
 * Formats this CRS as a <cite>Well Known Text</cite> {@code GeodeticCRS[…]} element.
 * More information about the WKT format is documented in subclasses.
 *
 * @return {@code "GeodeticCRS"} (WKT 2) or {@code "GeogCS"}/{@code "GeocCS"} (WKT 1).
 */
@Override
protected String formatTo(final Formatter formatter) {
    WKTUtilities.appendName(this, formatter, null);
    CoordinateSystem cs = getCoordinateSystem();
    final Convention convention = formatter.getConvention();
    final boolean isWKT1 = (convention.majorVersion() == 1);
    final boolean isGeographicWKT1 = isWKT1 && (cs instanceof EllipsoidalCS);
    if (isGeographicWKT1 && cs.getDimension() == 3) {
        /*
             * Version 1 of WKT format did not have three-dimensional GeographicCRS. Instead, such CRS were formatted
             * as a CompoundCRS made of a two-dimensional GeographicCRS with a VerticalCRS for the ellipsoidal height.
             * Note that such compound is illegal in WKT 2 and ISO 19111 standard, as ellipsoidal height shall not be
             * separated from the geographic component. So we perform this separation only at WKT 1 formatting time.
             */
        SingleCRS first = CRS.getHorizontalComponent(this);
        SingleCRS second = CRS.getVerticalComponent(this, true);
        if (first != null && second != null) {
            // Should not be null, but we are paranoiac.
            if (AxisDirection.UP.equals(AxisDirections.absolute(cs.getAxis(0).getDirection()))) {
                // It is very unusual to have VerticalCRS first, but our code tries to be robust.
                final SingleCRS t = first;
                first = second;
                second = t;
            }
            formatter.newLine();
            formatter.append(WKTUtilities.toFormattable(first));
            formatter.newLine();
            formatter.append(WKTUtilities.toFormattable(second));
            formatter.newLine();
            return WKTKeywords.Compd_CS;
        }
    }
    /*
         * Unconditionally format the datum element, followed by the prime meridian.
         * The prime meridian is part of datum according ISO 19111, but is formatted
         * as a sibling (rather than a child) element in WKT for historical reasons.
         */
    // Gives subclasses a chance to override.
    final GeodeticDatum datum = getDatum();
    formatter.newLine();
    formatter.append(WKTUtilities.toFormattable(datum));
    formatter.newLine();
    final PrimeMeridian pm = datum.getPrimeMeridian();
    final Unit<Angle> angularUnit = AxisDirections.getAngularUnit(cs, null);
    if (// Really this specific enum, not Convention.isSimplified().
    convention != Convention.WKT2_SIMPLIFIED || ReferencingUtilities.getGreenwichLongitude(pm, Units.DEGREE) != 0) {
        final Unit<Angle> oldUnit = formatter.addContextualUnit(angularUnit);
        formatter.indent(1);
        formatter.append(WKTUtilities.toFormattable(pm));
        formatter.indent(-1);
        formatter.newLine();
        formatter.restoreContextualUnit(angularUnit, oldUnit);
    }
    /*
         * Get the coordinate system to format. This will also determine the units to write and the keyword to
         * return in WKT 1 format. Note that for the WKT 1 format, we need to replace the coordinate system by
         * an instance conform to the legacy conventions.
         *
         * We can not delegate the work below to subclasses,  because XML unmarshalling of a geodetic CRS will
         * NOT create an instance of a subclass (because the distinction between geographic and geocentric CRS
         * is not anymore in ISO 19111:2007).
         */
    final boolean isBaseCRS;
    if (isWKT1) {
        if (!isGeographicWKT1) {
            // If not geographic, then presumed geocentric.
            if (cs instanceof CartesianCS) {
                cs = Legacy.forGeocentricCRS((CartesianCS) cs, true);
            } else {
                // SphericalCS was not supported in WKT 1.
                formatter.setInvalidWKT(cs, null);
            }
        }
        isBaseCRS = false;
    } else {
        isBaseCRS = isBaseCRS(formatter);
    }
    /*
         * Format the coordinate system, except if this CRS is the base CRS of an AbstractDerivedCRS in WKT 2 format.
         * This is because ISO 19162 omits the coordinate system definition of enclosed base CRS in order to simplify
         * the WKT. The 'formatCS(…)' method may write axis unit before or after the axes depending on whether we are
         * formatting WKT version 1 or 2 respectively.
         *
         * Note that even if we do not format the CS, we may still write the units if we are formatting in "simplified"
         * mode (as opposed to the more verbose mode). This looks like the opposite of what we would expect, but this is
         * because formatting the unit here allow us to avoid repeating the unit in projection parameters when this CRS
         * is part of a ProjectedCRS. Note however that in such case, the units to format are the angular units because
         * the linear units will be formatted in the enclosing PROJCS[…] element.
         */
    if (!isBaseCRS || convention == Convention.INTERNAL) {
        // Will also format the axes unit.
        formatCS(formatter, cs, ReferencingUtilities.getUnit(cs), isWKT1);
    } else if (convention.isSimplified()) {
        formatter.append(formatter.toContextualUnit(angularUnit));
    }
    /*
         * For WKT 1, the keyword depends on the subclass: "GeogCS" for GeographicCRS or "GeocCS" for GeocentricCRS.
         * However we can not rely on the subclass for choosing the keyword, because after XML unmarhaling we only
         * have a GeodeticCRS. We need to make the choice in this base class. The CS type is a sufficient criterion.
         */
    if (isWKT1) {
        return isGeographicWKT1 ? WKTKeywords.GeogCS : WKTKeywords.GeocCS;
    } else {
        return isBaseCRS ? WKTKeywords.BaseGeodCRS : formatter.shortOrLong(WKTKeywords.GeodCRS, WKTKeywords.GeodeticCRS);
    }
}
Also used : SingleCRS(org.opengis.referencing.crs.SingleCRS) CartesianCS(org.opengis.referencing.cs.CartesianCS) Convention(org.apache.sis.io.wkt.Convention) Angle(javax.measure.quantity.Angle) CoordinateSystem(org.opengis.referencing.cs.CoordinateSystem) EllipsoidalCS(org.opengis.referencing.cs.EllipsoidalCS) GeodeticDatum(org.opengis.referencing.datum.GeodeticDatum) PrimeMeridian(org.opengis.referencing.datum.PrimeMeridian)

Example 5 with Angle

use of javax.measure.quantity.Angle in project sis by apache.

the class DefaultProjectedCRS method formatTo.

/**
 * Formats the inner part of the <cite>Well Known Text</cite> (WKT) representation of this CRS.
 *
 * <div class="note"><b>Example:</b> Well-Known Text (version 2)
 * of a projected coordinate reference system using the Lambert Conformal method.
 *
 * {@preformat wkt
 *   ProjectedCRS[“NTF (Paris) / Lambert zone II”,
 *     BaseGeodCRS[“NTF (Paris)”,
 *       Datum[“Nouvelle Triangulation Francaise”,
 *         Ellipsoid[“NTF”, 6378249.2, 293.4660212936269, LengthUnit[“metre”, 1]]],
 *         PrimeMeridian[“Paris”, 2.5969213, AngleUnit[“grad”, 0.015707963267948967]]],
 *     Conversion[“Lambert zone II”,
 *       Method[“Lambert Conic Conformal (1SP)”, Id[“EPSG”, 9801, Citation[“IOGP”]]],
 *       Parameter[“Latitude of natural origin”, 52.0, AngleUnit[“grad”, 0.015707963267948967], Id[“EPSG”, 8801]],
 *       Parameter[“Longitude of natural origin”, 0.0, AngleUnit[“degree”, 0.017453292519943295], Id[“EPSG”, 8802]],
 *       Parameter[“Scale factor at natural origin”, 0.99987742, ScaleUnit[“unity”, 1], Id[“EPSG”, 8805]],
 *       Parameter[“False easting”, 600000.0, LengthUnit[“metre”, 1], Id[“EPSG”, 8806]],
 *       Parameter[“False northing”, 2200000.0, LengthUnit[“metre”, 1], Id[“EPSG”, 8807]]],
 *     CS[“Cartesian”, 2],
 *       Axis[“Easting (E)”, east, Order[1]],
 *       Axis[“Northing (N)”, north, Order[2]],
 *       LengthUnit[“metre”, 1],
 *     Id[“EPSG”, 27572, Citation[“IOGP”], URI[“urn:ogc:def:crs:EPSG::27572”]]]
 * }
 *
 * <p>Same coordinate reference system using WKT 1.</p>
 *
 * {@preformat wkt
 *   PROJCS[“NTF (Paris) / Lambert zone II”,
 *     GEOGCS[“NTF (Paris)”,
 *       DATUM[“Nouvelle Triangulation Francaise”,
 *         SPHEROID[“NTF”, 6378249.2, 293.4660212936269]],
 *         PRIMEM[“Paris”, 2.33722917],
 *       UNIT[“degree”, 0.017453292519943295],
 *       AXIS[“Longitude”, EAST],
 *       AXIS[“Latitude”, NORTH]],
 *     PROJECTION[“Lambert_Conformal_Conic_1SP”, AUTHORITY[“EPSG”, “9801”]],
 *     PARAMETER[“latitude_of_origin”, 46.8],
 *     PARAMETER[“central_meridian”, 0.0],
 *     PARAMETER[“scale_factor”, 0.99987742],
 *     PARAMETER[“false_easting”, 600000.0],
 *     PARAMETER[“false_northing”, 2200000.0],
 *     UNIT[“metre”, 1],
 *     AXIS[“Easting”, EAST],
 *     AXIS[“Northing”, NORTH],
 *     AUTHORITY[“EPSG”, “27572”]]
 * }
 * </div>
 *
 * @return {@code "ProjectedCRS"} (WKT 2) or {@code "ProjCS"} (WKT 1).
 *
 * @see <a href="http://docs.opengeospatial.org/is/12-063r5/12-063r5.html#57">WKT 2 specification §9</a>
 */
@Override
protected String formatTo(final Formatter formatter) {
    if (super.getConversionFromBase() == null) {
        /*
             * Should never happen except temporarily at construction time, or if the user invoked the copy constructor
             * with an invalid Conversion. Delegates to the super-class method for avoiding a NullPointerException.
             * That method returns 'null', which will cause the WKT to be declared invalid.
             */
        return super.formatTo(formatter);
    }
    WKTUtilities.appendName(this, formatter, null);
    final Convention convention = formatter.getConvention();
    final boolean isWKT1 = (convention.majorVersion() == 1);
    final CartesianCS cs = getCoordinateSystem();
    final GeographicCRS baseCRS = getBaseCRS();
    final Unit<?> lengthUnit = ReferencingUtilities.getUnit(cs);
    final Unit<Angle> angularUnit = AxisDirections.getAngularUnit(baseCRS.getCoordinateSystem(), null);
    final Unit<Angle> oldAngle = formatter.addContextualUnit(angularUnit);
    final Unit<?> oldLength = formatter.addContextualUnit(lengthUnit);
    /*
         * Format the enclosing base CRS. Note that WKT 1 formats a full GeographicCRS while WKT 2 formats only
         * the datum with the prime meridian (no coordinate system) and uses a different keyword ("BaseGeodCRS"
         * instead of "GeodeticCRS"). The DefaultGeodeticCRS.formatTo(Formatter) method detects when the CRS to
         * format is part of an enclosing ProjectedCRS and will adapt accordingly.
         */
    formatter.newLine();
    formatter.append(toFormattable(baseCRS));
    formatter.newLine();
    final Parameters p = new Parameters(this);
    final boolean isBaseCRS;
    if (isWKT1) {
        // Format outside of any "Conversion" element.
        p.append(formatter);
        isBaseCRS = false;
    } else {
        // Format inside a "Conversion" element.
        formatter.append(p);
        isBaseCRS = isBaseCRS(formatter);
    }
    /*
         * In WKT 2 format, the coordinate system axes are written only if this projected CRS is not the base CRS
         * of another derived CRS.
         */
    if (!isBaseCRS || convention == Convention.INTERNAL) {
        formatCS(formatter, cs, lengthUnit, isWKT1);
    }
    formatter.restoreContextualUnit(lengthUnit, oldLength);
    formatter.restoreContextualUnit(angularUnit, oldAngle);
    return isWKT1 ? WKTKeywords.ProjCS : isBaseCRS ? WKTKeywords.BaseProjCRS : formatter.shortOrLong(WKTKeywords.ProjCRS, WKTKeywords.ProjectedCRS);
}
Also used : CartesianCS(org.opengis.referencing.cs.CartesianCS) Convention(org.apache.sis.io.wkt.Convention) AxesConvention(org.apache.sis.referencing.cs.AxesConvention) Angle(javax.measure.quantity.Angle) GeographicCRS(org.opengis.referencing.crs.GeographicCRS)

Aggregations

Angle (javax.measure.quantity.Angle)19 Length (javax.measure.quantity.Length)7 Path (java.nio.file.Path)3 Unit (javax.measure.Unit)3 Convention (org.apache.sis.io.wkt.Convention)3 CartesianCS (org.opengis.referencing.cs.CartesianCS)3 GeodeticDatum (org.opengis.referencing.datum.GeodeticDatum)3 FactoryException (org.opengis.util.FactoryException)3 BufferedReader (java.io.BufferedReader)2 BigDecimal (java.math.BigDecimal)2 Matrix3 (org.apache.sis.referencing.operation.matrix.Matrix3)2 Test (org.junit.Test)2 Envelope (org.opengis.geometry.Envelope)2 Identifier (org.opengis.metadata.Identifier)2 ParameterValueGroup (org.opengis.parameter.ParameterValueGroup)2 GeographicCRS (org.opengis.referencing.crs.GeographicCRS)2 CoordinateSystem (org.opengis.referencing.cs.CoordinateSystem)2 EllipsoidalCS (org.opengis.referencing.cs.EllipsoidalCS)2 NoninvertibleTransformException (org.opengis.referencing.operation.NoninvertibleTransformException)2 GeographicCoordinateSystem (com.revolsys.geometry.cs.GeographicCoordinateSystem)1