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;
}
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));
}
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));
}
}
}
}
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);
}
}
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);
}
Aggregations