Examples with DefaultGeodeticDatum org.apache.sis.referencing.datum.DefaultGeodeticDatum used on opensource projects

Example 1 with DefaultGeodeticDatum

use of org.apache.sis.referencing.datum.DefaultGeodeticDatum in project sis by apache.

the class GeodeticObjectParserTest method testProjectedWithID.

/**
 * Tests the parsing of a projected CRS from a WKT 1 string with authority and Bursa-Wolf parameters.
 *
 * @throws ParseException if the parsing failed.
 */
@Test
@DependsOnMethod("testProjectedCRS")
public void testProjectedWithID() throws ParseException {
    final ProjectedCRS crs = parse(ProjectedCRS.class, "PROJCS[“OSGB 1936 / British National Grid”,\n" + "  GEOGCS[“OSGB 1936”,\n" + "    DATUM[“OSGB_1936”,\n" + "      SPHEROID[“Airy 1830”, 6377563.396, 299.3249646, AUTHORITY[“EPSG”, “7001”]],\n" + "      TOWGS84[375.0, -111.0, 431.0, 0.0, 0.0, 0.0, 0.0],\n" + "      AUTHORITY[“EPSG”, “6277”]],\n" + "      PRIMEM[“Greenwich”,0.0, AUTHORITY[“EPSG”, “8901”]],\n" + "    UNIT[“DMSH”,0.0174532925199433],\n" + "    AXIS[“Lat”,NORTH],AXIS[“Long”,EAST], AUTHORITY[“EPSG”, “4277”]],\n" + "  PROJECTION[“Transverse_Mercator”],\n" + "  PARAMETER[“latitude_of_origin”, 49.0],\n" + "  PARAMETER[“central_meridian”, -2.0],\n" + "  PARAMETER[“scale_factor”, 0.999601272],\n" + "  PARAMETER[“false_easting”, 400000.0],\n" + "  PARAMETER[“false_northing”, -100000.0],\n" + "  UNIT[“metre”, 1.0, AUTHORITY[“EPSG”, “9001”]],\n" + "  AXIS[“E”,EAST],\n" + "  AXIS[“N”,NORTH],\n" + "  AUTHORITY[“EPSG”, “27700”]]");
    assertNameAndIdentifierEqual("OSGB 1936 / British National Grid", 27700, crs);
    assertNameAndIdentifierEqual("OSGB 1936", 4277, crs.getBaseCRS());
    assertNameAndIdentifierEqual("OSGB_1936", 6277, crs.getDatum());
    verifyProjectedCS(crs.getCoordinateSystem(), Units.METRE);
    final ParameterValueGroup param = crs.getConversionFromBase().getParameterValues();
    assertEquals("Transverse Mercator", crs.getConversionFromBase().getMethod().getName().getCode());
    assertEquals("semi_major", 6377563.396, param.parameter("semi_major").doubleValue(), 1E-4);
    assertEquals("semi_minor", 6356256.909, param.parameter("semi_minor").doubleValue(), 1E-3);
    assertEquals("latitude_of_origin", 49.0, param.parameter("latitude_of_origin").doubleValue(), 1E-8);
    assertEquals("central_meridian", -2.0, param.parameter("central_meridian").doubleValue(), 1E-8);
    assertEquals("scale_factor", 0.9996, param.parameter("scale_factor").doubleValue(), 1E-5);
    assertEquals("false_easting", 400000.0, param.parameter("false_easting").doubleValue(), 1E-4);
    assertEquals("false_northing", -100000.0, param.parameter("false_northing").doubleValue(), 1E-4);
    final BursaWolfParameters[] bwp = ((DefaultGeodeticDatum) crs.getDatum()).getBursaWolfParameters();
    assertEquals("BursaWolfParameters", 1, bwp.length);
    assertArrayEquals("BursaWolfParameters", new double[] { 375, -111, 431 }, bwp[0].getValues(), STRICT);
}

Example 2 with DefaultGeodeticDatum

use of org.apache.sis.referencing.datum.DefaultGeodeticDatum in project sis by apache.

the class GeodeticObjectParserTest method testProjectedWithMissingName.

/**
 * Tests parsing a WKT with a missing Geographic CRS name.
 * This should be considered invalid, but happen in practice.
 *
 * <p>The WKT tested in this method contains also some other oddities compared to the usual WKT:</p>
 * <ul>
 *   <li>The prime meridian is declared in the {@code "central_meridian"} projection parameter instead
 *       than in the {@code PRIMEM[…]} element.</li>
 *   <li>Some elements are not in the usual order.</li>
 * </ul>
 *
 * @throws ParseException if the parsing failed.
 */
@Test
@DependsOnMethod("testProjectedCRS")
public void testProjectedWithMissingName() throws ParseException {
    final ProjectedCRS crs = parse(ProjectedCRS.class, "PROJCS[“FRANCE/NTF/Lambert III”," + // Missing name (the purpose of this test).
    "GEOGCS[“”," + // Intentionally misplaced coma.
    "DATUM[“NTF=GR3DF97A”,TOWGS84[-168, -60, 320] ," + "SPHEROID[“Clarke 1880 (IGN)”,6378249.2,293.4660212936269]]," + "PRIMEM[“Greenwich”,0],UNIT[“Degrees”,0.0174532925199433]," + "AXIS[“Long”,East],AXIS[“Lat”,North]]," + "PROJECTION[“Lambert_Conformal_Conic_1SP”]," + "PARAMETER[“latitude_of_origin”,44.1]," + // Paris prime meridian.
    "PARAMETER[“central_meridian”,2.33722917]," + "PARAMETER[“scale_factor”,0.999877499]," + "PARAMETER[“false_easting”,600000]," + "PARAMETER[“false_northing”,200000]," + "UNIT[“Meter”,1]," + "AXIS[“Easting”,East],AXIS[“Northing”,North]]");
    assertNameAndIdentifierEqual("FRANCE/NTF/Lambert III", 0, crs);
    verifyProjectedCS(crs.getCoordinateSystem(), Units.METRE);
    final GeographicCRS geoCRS = crs.getBaseCRS();
    // Inherited the datum name.
    assertNameAndIdentifierEqual("NTF=GR3DF97A", 0, geoCRS);
    final GeodeticDatum datum = geoCRS.getDatum();
    assertNameAndIdentifierEqual("NTF=GR3DF97A", 0, datum);
    assertNameAndIdentifierEqual("Greenwich", 0, datum.getPrimeMeridian());
    assertArrayEquals("BursaWolfParameters", new double[] { -168, -60, 320 }, ((DefaultGeodeticDatum) datum).getBursaWolfParameters()[0].getValues(), STRICT);
    final Ellipsoid ellipsoid = datum.getEllipsoid();
    assertNameAndIdentifierEqual("Clarke 1880 (IGN)", 0, ellipsoid);
    assertEquals("semiMajor", 6378249.2, ellipsoid.getSemiMajorAxis(), STRICT);
    assertEquals("inverseFlattening", 293.4660212936269, ellipsoid.getInverseFlattening(), STRICT);
    final EllipsoidalCS cs = geoCRS.getCoordinateSystem();
    assertEquals("dimension", 2, cs.getDimension());
    assertLongitudeAxisEquals(cs.getAxis(0));
    assertLatitudeAxisEquals(cs.getAxis(1));
    final ParameterValueGroup param = crs.getConversionFromBase().getParameterValues();
    assertEquals("Lambert Conic Conformal (1SP)", param.getDescriptor().getName().getCode());
    assertEquals("semi_major", 6378249.2, param.parameter("semi_major").doubleValue(Units.METRE), STRICT);
    assertEquals("semi_minor", 6356515.0, param.parameter("semi_minor").doubleValue(Units.METRE), 1E-12);
    assertEquals("latitude_of_origin", 44.1, param.parameter("latitude_of_origin").doubleValue(Units.DEGREE), STRICT);
    assertEquals("central_meridian", 2.33722917, param.parameter("central_meridian").doubleValue(Units.DEGREE), STRICT);
    assertEquals("scale_factor", 0.999877499, param.parameter("scale_factor").doubleValue(Units.UNITY), STRICT);
    assertEquals("false_easting", 600000.0, param.parameter("false_easting").doubleValue(Units.METRE), STRICT);
    assertEquals("false_northing", 200000.0, param.parameter("false_northing").doubleValue(Units.METRE), STRICT);
}

Example 3 with DefaultGeodeticDatum

use of org.apache.sis.referencing.datum.DefaultGeodeticDatum in project sis by apache.

the class CoordinateOperationFinder method createOperationStep.

/**
 * Creates an operation between two geodetic (geographic or geocentric) coordinate reference systems.
 * The default implementation can:
 *
 * <ul>
 *   <li>adjust axis order and orientation, for example converting from (<cite>North</cite>, <cite>West</cite>)
 *       axes to (<cite>East</cite>, <cite>North</cite>) axes,</li>
 *   <li>apply units conversion if needed,</li>
 *   <li>perform longitude rotation if needed,</li>
 *   <li>perform datum shift if {@linkplain BursaWolfParameters Bursa-Wolf parameters} are available
 *       for the area of interest.</li>
 * </ul>
 *
 * <p>This method returns only <em>one</em> step for a chain of concatenated operations (to be built by the caller).
 * But a list is returned because the same step may be implemented by different operation methods. Only one element
 * in the returned list should be selected (usually the first one).</p>
 *
 * @param  sourceCRS  input coordinate reference system.
 * @param  targetCRS  output coordinate reference system.
 * @return a coordinate operation from {@code sourceCRS} to {@code targetCRS}.
 * @throws FactoryException if the operation can not be constructed.
 */
@SuppressWarnings("null")
protected List<CoordinateOperation> createOperationStep(final GeodeticCRS sourceCRS, final GeodeticCRS targetCRS) throws FactoryException {
    final GeodeticDatum sourceDatum = sourceCRS.getDatum();
    final GeodeticDatum targetDatum = targetCRS.getDatum();
    Matrix datumShift = null;
    /*
         * If the prime meridian is not the same, we will concatenate a longitude rotation before or after datum shift
         * (that concatenation will be performed by the customized DefaultMathTransformFactory.Context created below).
         * Actually we do not know if the longitude rotation should be before or after datum shift. But this ambiguity
         * can usually be ignored because Bursa-Wolf parameters are always used with source and target prime meridians
         * set to Greenwich in EPSG dataset 8.9.  For safety, the SIS's DefaultGeodeticDatum class ensures that if the
         * prime meridians are not the same, then the target meridian must be Greenwich.
         */
    final DefaultMathTransformFactory.Context context = ReferencingUtilities.createTransformContext(sourceCRS, targetCRS, new MathTransformContext(sourceDatum, targetDatum));
    /*
         * If both CRS use the same datum and the same prime meridian, then the coordinate operation is only axis
         * swapping, unit conversion or change of coordinate system type (Ellipsoidal ↔ Cartesian ↔ Spherical).
         * Otherwise (if the datum are not the same), we will need to perform a scale, translation and rotation
         * in Cartesian space using the Bursa-Wolf parameters. If the user does not require the best accuracy,
         * then the Molodensky approximation may be used for avoiding the conversion step to geocentric CRS.
         */
    Identifier identifier;
    boolean isGeographicToGeocentric = false;
    final CoordinateSystem sourceCS = context.getSourceCS();
    final CoordinateSystem targetCS = context.getTargetCS();
    if (equalsIgnoreMetadata(sourceDatum, targetDatum)) {
        final boolean isGeocentricToGeographic;
        isGeographicToGeocentric = (sourceCS instanceof EllipsoidalCS && targetCS instanceof CartesianCS);
        isGeocentricToGeographic = (sourceCS instanceof CartesianCS && targetCS instanceof EllipsoidalCS);
        /*
             * Above booleans should never be true in same time. If it nevertheless happen (we are paranoiac;
             * maybe a lazy user implemented all interfaces in a single class), do not apply any geographic ↔
             * geocentric conversion. Instead do as if the coordinate system types were the same.
             */
        if (isGeocentricToGeographic ^ isGeographicToGeocentric) {
            identifier = GEOCENTRIC_CONVERSION;
        } else {
            identifier = AXIS_CHANGES;
        }
    } else {
        identifier = ELLIPSOID_CHANGE;
        if (sourceDatum instanceof DefaultGeodeticDatum) {
            datumShift = ((DefaultGeodeticDatum) sourceDatum).getPositionVectorTransformation(targetDatum, areaOfInterest);
            if (datumShift != null) {
                identifier = DATUM_SHIFT;
            }
        }
    }
    /*
         * Conceptually, all transformations below could done by first converting from the source coordinate
         * system to geocentric Cartesian coordinates (X,Y,Z), apply an affine transform represented by the
         * datum shift matrix, then convert from the (X′,Y′,Z′) coordinates to the target coordinate system.
         * However there is two exceptions to this path:
         *
         *   1) In the particular where both the source and target CS are ellipsoidal, we may use the
         *      Molodensky approximation as a shortcut (if the desired accuracy allows).
         *
         *   2) Even if we really go through the XYZ coordinates without Molodensky approximation, there is
         *      at least 9 different ways to name this operation depending on whether the source and target
         *      CRS are geocentric or geographic, 2- or 3-dimensional, whether there is a translation or not,
         *      the rotation sign, etc. We try to use the most specific name if we can find one, and fallback
         *      on an arbitrary name only in last resort.
         */
    final DefaultMathTransformFactory mtFactory = factorySIS.getDefaultMathTransformFactory();
    MathTransform before = null, after = null;
    ParameterValueGroup parameters;
    if (identifier == DATUM_SHIFT || identifier == ELLIPSOID_CHANGE) {
        /*
             * If the transform can be represented by a single coordinate operation, returns that operation.
             * Possible operations are:
             *
             *    - Position Vector transformation (in geocentric, geographic-2D or geographic-3D domains)
             *    - Geocentric translation         (in geocentric, geographic-2D or geographic-3D domains)
             *    - [Abridged] Molodensky          (as an approximation of geocentric translation)
             *    - Identity                       (if the desired accuracy is so large than we can skip datum shift)
             *
             * TODO: if both CS are ellipsoidal but with different number of dimensions, then we should use
             * an intermediate 3D geographic CRS in order to enable the use of Molodensky method if desired.
             */
        final DatumShiftMethod preferredMethod = DatumShiftMethod.forAccuracy(desiredAccuracy);
        parameters = GeocentricAffine.createParameters(sourceCS, targetCS, datumShift, preferredMethod);
        if (parameters == null) {
            /*
                 * Failed to select a coordinate operation. Maybe because the coordinate system types are not the same.
                 * Convert unconditionally to XYZ geocentric coordinates and apply the datum shift in that CS space.
                 *
                 * TODO: operation name should not be "Affine" if 'before' or 'after' transforms are not identity.
                 *       Reminder: the parameter group name here determines the OperationMethod later in this method.
                 */
            if (datumShift != null) {
                parameters = TensorParameters.WKT1.createValueGroup(properties(Constants.AFFINE), datumShift);
            } else {
                // Dimension of geocentric CRS.
                parameters = Affine.identity(3);
            }
            final CoordinateSystem normalized = CommonCRS.WGS84.geocentric().getCoordinateSystem();
            before = mtFactory.createCoordinateSystemChange(sourceCS, normalized, sourceDatum.getEllipsoid());
            after = mtFactory.createCoordinateSystemChange(normalized, targetCS, targetDatum.getEllipsoid());
            context.setSource(normalized);
            context.setTarget(normalized);
        }
    } else if (identifier == GEOCENTRIC_CONVERSION) {
        parameters = (isGeographicToGeocentric ? GeographicToGeocentric.PARAMETERS : GeocentricToGeographic.PARAMETERS).createValue();
    } else {
        /*
             * Coordinate system change (including change in the number of dimensions) without datum shift.
             */
        final int sourceDim = sourceCS.getDimension();
        final int targetDim = targetCS.getDimension();
        if (// sourceDim == 2 or 3.
        (sourceDim & ~1) == 2 && // abs(sourceDim - targetDim) == 1.
        (sourceDim ^ targetDim) == 1 && (sourceCS instanceof EllipsoidalCS) && (targetCS instanceof EllipsoidalCS)) {
            parameters = (sourceDim == 2 ? Geographic2Dto3D.PARAMETERS : Geographic3Dto2D.PARAMETERS).createValue();
        } else {
            /*
                 * TODO: instead than creating parameters for an identity operation, we should create the
                 *       CoordinateOperation directly from the MathTransform created by mtFactory below.
                 *       The intent if to get the correct OperationMethod, which should not be "Affine"
                 *       if there is a CS type change.
                 */
            parameters = Affine.identity(targetDim);
            /*
                 * createCoordinateSystemChange(…) needs the ellipsoid associated to the ellipsoidal coordinate system,
                 * if any. If none or both coordinate systems are ellipsoidal, then the ellipsoid will be ignored (see
                 * createCoordinateSystemChange(…) javadoc for the rational) so it does not matter which one we pick.
                 */
            before = mtFactory.createCoordinateSystemChange(sourceCS, targetCS, (sourceCS instanceof EllipsoidalCS ? sourceDatum : targetDatum).getEllipsoid());
            context.setSource(targetCS);
        }
    }
    /*
         * Transform between differents datums using Bursa Wolf parameters. The Bursa Wolf parameters are used
         * with "standard" geocentric CS, i.e. with X axis towards the prime meridian, Y axis towards East and
         * Z axis toward North, unless the Molodensky approximation is used. The following steps are applied:
         *
         *     source CRS                        →
         *     normalized CRS with source datum  →
         *     normalized CRS with target datum  →
         *     target CRS
         *
         * Those steps may be either explicit with the 'before' and 'after' transform, or implicit with the
         * Context parameter.
         */
    MathTransform transform = mtFactory.createParameterizedTransform(parameters, context);
    final OperationMethod method = mtFactory.getLastMethodUsed();
    if (before != null) {
        transform = mtFactory.createConcatenatedTransform(before, transform);
        if (after != null) {
            transform = mtFactory.createConcatenatedTransform(transform, after);
        }
    }
    return asList(createFromMathTransform(properties(identifier), sourceCRS, targetCRS, transform, method, parameters, null));
}

Example 4 with DefaultGeodeticDatum

use of org.apache.sis.referencing.datum.DefaultGeodeticDatum in project sis by apache.

the class EPSGFactoryTest method testGeographic2D.

/**
 * Tests the "Datum 73" geographic CRS (EPSG:4274), which has a datum different than the WGS84 one.
 *
 * @throws FactoryException if an error occurred while querying the factory.
 */
@Test
@DependsOnMethod("testWGS84")
public void testGeographic2D() throws FactoryException {
    final EPSGFactory factory = TestFactorySource.factory;
    assumeNotNull(factory);
    final GeographicCRS crs = factory.createGeographicCRS("4274");
    assertEpsgNameAndIdentifierEqual("Datum 73", 4274, crs);
    assertEpsgNameAndIdentifierEqual("Datum 73", 6274, crs.getDatum());
    assertAxisDirectionsEqual("EPSG::6422", crs.getCoordinateSystem(), AxisDirection.NORTH, AxisDirection.EAST);
    final BursaWolfParameters[] bwp = ((DefaultGeodeticDatum) crs.getDatum()).getBursaWolfParameters();
    assertTrue("Expected a transformation to WGS84.", bwp.length >= 1);
    assertSame("CRS shall be cached", crs, factory.createCoordinateReferenceSystem("4274"));
}

Example 5 with DefaultGeodeticDatum

use of org.apache.sis.referencing.datum.DefaultGeodeticDatum in project sis by apache.

the class EPSGFactoryTest method testWGS84.

/**
 * Tests the "WGS 84" geographic CRS (EPSG:4326).
 *
 * @throws FactoryException if an error occurred while querying the factory.
 */
@Test
public void testWGS84() throws FactoryException {
    final EPSGFactory factory = TestFactorySource.factory;
    assumeNotNull(factory);
    final GeographicCRS crs = factory.createGeographicCRS("EPSG:4326");
    assertEpsgNameAndIdentifierEqual("WGS 84", 4326, crs);
    assertEpsgNameAndIdentifierEqual("World Geodetic System 1984", 6326, crs.getDatum());
    assertAxisDirectionsEqual("EPSG:6422", crs.getCoordinateSystem(), AxisDirection.NORTH, AxisDirection.EAST);
    final BursaWolfParameters[] bwp = ((DefaultGeodeticDatum) crs.getDatum()).getBursaWolfParameters();
    assertEquals("Expected no Bursa-Wolf parameters.", 0, bwp.length);
    assertSame("CRS shall be cached", crs, factory.createCoordinateReferenceSystem("4326"));
    assertSame("Shall accept \"::\"", crs, factory.createGeographicCRS("EPSG::4326"));
}