Search in sources :

Example 1 with BursaWolfParameters

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

the class GeocentricAffine method createParameters.

/**
 * Returns the parameters for creating a datum shift operation.
 * The operation method will be one of the {@code GeocentricAffine} subclasses,
 * unless the specified {@code method} argument is {@link DatumShiftMethod#NONE}.
 * If no single operation method can be used, then this method returns {@code null}.
 *
 * <p>This method does <strong>not</strong> change the coordinate system type.
 * The source and target coordinate systems can be both {@code EllipsoidalCS} or both {@code CartesianCS}.
 * Any other type or mix of types (e.g. a {@code EllipsoidalCS} source and {@code CartesianCS} target)
 * will cause this method to return {@code null}. In such case, it is caller's responsibility to apply
 * the datum shift itself in Cartesian geocentric coordinates.</p>
 *
 * @param  sourceCS    the source coordinate system. Only the type and number of dimensions is checked.
 * @param  targetCS    the target coordinate system. Only the type and number of dimensions is checked.
 * @param  datumShift  the datum shift as a matrix, or {@code null} if there is no datum shift information.
 * @param  method      the preferred datum shift method. Note that {@code createParameters(…)} may overwrite.
 * @return the parameter values, or {@code null} if no single operation method can be found.
 */
public static ParameterValueGroup createParameters(final CoordinateSystem sourceCS, final CoordinateSystem targetCS, final Matrix datumShift, DatumShiftMethod method) {
    final boolean isEllipsoidal = (sourceCS instanceof EllipsoidalCS);
    if (!(isEllipsoidal ? (targetCS instanceof EllipsoidalCS) : (targetCS instanceof CartesianCS && sourceCS instanceof CartesianCS))) {
        // Coordinate systems are not two EllipsoidalCS or two CartesianCS.
        return null;
    }
    @SuppressWarnings("null") int dimension = sourceCS.getDimension();
    if (dimension != targetCS.getDimension()) {
        // Any value greater than 3 means "mismatched dimensions" for this method.
        dimension = 4;
    }
    if (method == DatumShiftMethod.NONE) {
        if (dimension <= 3) {
            return Affine.identity(dimension);
        } else if (isEllipsoidal) {
            final ParameterDescriptorGroup descriptor;
            switch(sourceCS.getDimension()) {
                case 2:
                    descriptor = Geographic2Dto3D.PARAMETERS;
                    break;
                case 3:
                    descriptor = Geographic3Dto2D.PARAMETERS;
                    break;
                default:
                    return null;
            }
            return descriptor.createValue();
        } else {
            return null;
        }
    }
    /*
         * Try to convert the matrix into (tX, tY, tZ, rX, rY, rZ, dS) parameters.
         * The matrix may not be convertible, in which case we will let the caller
         * uses the matrix directly in Cartesian geocentric coordinates.
         */
    final BursaWolfParameters parameters = new BursaWolfParameters(null, null);
    if (datumShift != null)
        try {
            parameters.setPositionVectorTransformation(datumShift, BURSAWOLF_TOLERANCE);
        } catch (IllegalArgumentException e) {
            log(Loggers.COORDINATE_OPERATION, "createParameters", e);
            return null;
        }
    else {
        /*
             * If there is no datum shift parameters (not to be confused with identity), then those parameters
             * are assumed unknown. Using the most accurate methods would give a false impression of accuracy,
             * so we use the fastest method instead. Since all parameter values are zero, Apache SIS should use
             * the AbridgedMolodenskyTransform2D optimization.
             */
        method = DatumShiftMethod.ABRIDGED_MOLODENSKY;
    }
    final boolean isTranslation = parameters.isTranslation();
    final ParameterDescriptorGroup descriptor;
    /*
         * Following "if" blocks are ordered from most accurate to less accurate datum shift method
         * supported by GeocentricAffine subclasses (except NONE which has already been handled).
         * Special cases:
         *
         *   - If the datum shift is applied between geocentric CRS, then the Molodensky approximations do not apply
         *     as they are designed for transformations between geographic CRS only. User preference is then ignored.
         *
         *   - Molodensky methods are approximations for datum shifts having only translation terms in their Bursa-Wolf
         *     parameters. If there is also a scale or rotation terms, then we can not use Molodensky methods. The user
         *     preference is then ignored.
         */
    if (!isEllipsoidal) {
        method = DatumShiftMethod.GEOCENTRIC_DOMAIN;
        descriptor = isTranslation ? GeocentricTranslation.PARAMETERS : PositionVector7Param.PARAMETERS;
    } else if (!isTranslation) {
        method = DatumShiftMethod.GEOCENTRIC_DOMAIN;
        descriptor = (dimension >= 3) ? PositionVector7Param3D.PARAMETERS : PositionVector7Param2D.PARAMETERS;
    } else
        switch(method) {
            case GEOCENTRIC_DOMAIN:
                {
                    descriptor = (dimension >= 3) ? GeocentricTranslation3D.PARAMETERS : GeocentricTranslation2D.PARAMETERS;
                    break;
                }
            case MOLODENSKY:
                {
                    descriptor = Molodensky.PARAMETERS;
                    break;
                }
            case ABRIDGED_MOLODENSKY:
                {
                    descriptor = AbridgedMolodensky.PARAMETERS;
                    break;
                }
            default:
                throw new AssertionError(method);
        }
    /*
         * Following lines will set all Bursa-Wolf parameter values (scale, translation
         * and rotation terms). In the particular case of Molodensky method, we have an
         * additional parameter for the number of source and target dimensions (2 or 3).
         */
    final Parameters values = createParameters(descriptor, parameters, isTranslation);
    switch(method) {
        case MOLODENSKY:
        case ABRIDGED_MOLODENSKY:
            {
                if (dimension <= 3) {
                    values.getOrCreate(Molodensky.DIMENSION).setValue(dimension);
                }
                break;
            }
    }
    return values;
}
Also used : CartesianCS(org.opengis.referencing.cs.CartesianCS) Parameters(org.apache.sis.parameter.Parameters) BursaWolfParameters(org.apache.sis.referencing.datum.BursaWolfParameters) ParameterDescriptorGroup(org.opengis.parameter.ParameterDescriptorGroup) EllipsoidalCS(org.opengis.referencing.cs.EllipsoidalCS) BursaWolfParameters(org.apache.sis.referencing.datum.BursaWolfParameters)

Example 2 with BursaWolfParameters

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

the class GeocentricAffine method asDatumShift.

/**
 * Given a transformation chain, conditionally replaces the affine transform elements by an alternative object
 * showing the Bursa-Wolf parameters. The replacement is applied if and only if the affine transform is a scale,
 * translation or rotation in the geocentric domain.
 *
 * <p>This method is invoked only by {@code ConcatenatedTransform.getPseudoSteps()} for the need of WKT formatting.
 * The purpose of this method is very similar to the purpose of {@code AbstractMathTransform.beforeFormat(List, int,
 * boolean)} except that we need to perform the {@code forDatumShift(…)} work only after {@code beforeFormat(…)}
 * finished its work for all {@code ContextualParameters}, including the {@code EllipsoidToCentricTransform}'s one.</p>
 *
 * @param  transforms  the full chain of concatenated transforms.
 */
public static void asDatumShift(final List<Object> transforms) {
    for (int i = transforms.size() - 2; --i >= 0; ) {
        if (isOperation(GeographicToGeocentric.NAME, transforms.get(i)) && isOperation(GeocentricToGeographic.NAME, transforms.get(i + 2))) {
            final Object step = transforms.get(i + 1);
            if (step instanceof LinearTransform) {
                final BursaWolfParameters parameters = new BursaWolfParameters(null, null);
                try {
                    parameters.setPositionVectorTransformation(((LinearTransform) step).getMatrix(), BURSAWOLF_TOLERANCE);
                } catch (IllegalArgumentException e) {
                    /*
                         * Should not occur, except sometime on inverse transform of relatively complex datum shifts
                         * (more than just translation terms). We can fallback on formatting the full matrix.
                         */
                    log(Loggers.WKT, "asDatumShift", e);
                    continue;
                }
                final boolean isTranslation = parameters.isTranslation();
                final Parameters values = createParameters(isTranslation ? GeocentricTranslation.PARAMETERS : PositionVector7Param.PARAMETERS, parameters, isTranslation);
                transforms.set(i + 1, new FormattableObject() {

                    @Override
                    protected String formatTo(final Formatter formatter) {
                        WKTUtilities.appendParamMT(values, formatter);
                        return WKTKeywords.Param_MT;
                    }
                });
            }
        }
    }
}
Also used : Parameters(org.apache.sis.parameter.Parameters) BursaWolfParameters(org.apache.sis.referencing.datum.BursaWolfParameters) Formatter(org.apache.sis.io.wkt.Formatter) FormattableObject(org.apache.sis.io.wkt.FormattableObject) BursaWolfParameters(org.apache.sis.referencing.datum.BursaWolfParameters) LinearTransform(org.apache.sis.referencing.operation.transform.LinearTransform) FormattableObject(org.apache.sis.io.wkt.FormattableObject)

Example 3 with BursaWolfParameters

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

the class GeocentricAffine method createMathTransform.

/**
 * Creates a math transform from the specified group of parameter values.
 * The default implementation creates an affine transform, but some subclasses
 * will wrap that affine operation into Geographic/Geocentric conversions.
 *
 * @param  factory  the factory to use for creating concatenated transforms.
 * @param  values   the group of parameter values.
 * @return the created math transform.
 * @throws FactoryException if a transform can not be created.
 */
@Override
@SuppressWarnings("fallthrough")
public MathTransform createMathTransform(final MathTransformFactory factory, final ParameterValueGroup values) throws FactoryException {
    final BursaWolfParameters parameters = new BursaWolfParameters(null, null);
    final Parameters pv = Parameters.castOrWrap(values);
    boolean reverseRotation = false;
    switch(getType()) {
        default:
            throw new AssertionError();
        // Fall through
        case FRAME_ROTATION:
            reverseRotation = true;
        case SEVEN_PARAM:
            parameters.rX = pv.doubleValue(RX);
            parameters.rY = pv.doubleValue(RY);
            parameters.rZ = pv.doubleValue(RZ);
            parameters.dS = pv.doubleValue(DS);
        case // Fall through
        TRANSLATION:
            // Fall through
            parameters.tX = pv.doubleValue(TX);
            parameters.tY = pv.doubleValue(TY);
            parameters.tZ = pv.doubleValue(TZ);
    }
    if (reverseRotation) {
        parameters.reverseRotation();
    }
    return MathTransforms.linear(parameters.getPositionVectorTransformation(null));
}
Also used : Parameters(org.apache.sis.parameter.Parameters) BursaWolfParameters(org.apache.sis.referencing.datum.BursaWolfParameters) BursaWolfParameters(org.apache.sis.referencing.datum.BursaWolfParameters)

Example 4 with BursaWolfParameters

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

the class EPSGDataAccess method createBursaWolfParameters.

/**
 * Returns Bursa-Wolf parameters for a geodetic datum. If the specified datum has no conversion informations,
 * then this method returns {@code null}.
 *
 * <p>This method is for compatibility with <cite>Well Known Text</cite> (WKT) version 1 formatting.
 * That legacy format had a {@code TOWGS84} element which needs the information provided by this method.
 * Note that {@code TOWGS84} is a deprecated element as of WKT 2 (ISO 19162).</p>
 *
 * @param  meridian  the source datum prime meridian, used for discarding any target datum using a different meridian.
 * @param  code      the EPSG code of the source {@link GeodeticDatum}.
 * @return an array of Bursa-Wolf parameters, or {@code null}.
 */
private BursaWolfParameters[] createBursaWolfParameters(final PrimeMeridian meridian, final Integer code) throws SQLException, FactoryException {
    /*
         * We do not provide TOWGS84 information for WGS84 itself or for any other datum on our list of target datum,
         * in order to avoid infinite recursivity. The 'ensureNonRecursive' call is an extra safety check which should
         * never fail, unless TARGET_CRS and TARGET_DATUM values do not agree with database content.
         */
    if (code == BursaWolfInfo.TARGET_DATUM) {
        return null;
    }
    final List<BursaWolfInfo> bwInfos = new ArrayList<>();
    try (ResultSet result = executeQuery("BursaWolfParametersSet", "SELECT COORD_OP_CODE," + " COORD_OP_METHOD_CODE," + " TARGET_CRS_CODE," + " AREA_OF_USE_CODE" + " FROM [Coordinate_Operation]" + // Do not put spaces around "=" - SQLTranslator searches for this exact match.
    " WHERE DEPRECATED=0" + " AND TARGET_CRS_CODE = " + BursaWolfInfo.TARGET_CRS + " AND COORD_OP_METHOD_CODE >= " + BursaWolfInfo.MIN_METHOD_CODE + " AND COORD_OP_METHOD_CODE <= " + BursaWolfInfo.MAX_METHOD_CODE + " AND SOURCE_CRS_CODE IN " + "(SELECT COORD_REF_SYS_CODE FROM [Coordinate Reference System] WHERE DATUM_CODE = ?)" + " ORDER BY TARGET_CRS_CODE, COORD_OP_ACCURACY, COORD_OP_CODE DESC", code)) {
        while (result.next()) {
            final BursaWolfInfo info = new BursaWolfInfo(// Operation
            getInteger(code, result, 1), // Method
            getInteger(code, result, 2), // Target datum
            getInteger(code, result, 3), // Domain of validity
            getInteger(code, result, 4));
            if (info.target != code) {
                // Paranoiac check.
                bwInfos.add(info);
            }
        }
    }
    int size = bwInfos.size();
    if (size == 0) {
        return null;
    }
    /*
         * Sort the infos in preference order. The "ORDER BY" clause above was not enough;
         * we also need to take the "Supersession" table in account. Once the sorting is done,
         * keep only one Bursa-Wolf parameters for each datum.
         */
    if (size > 1) {
        final BursaWolfInfo[] codes = bwInfos.toArray(new BursaWolfInfo[size]);
        sort("Coordinate_Operation", codes);
        bwInfos.clear();
        BursaWolfInfo.filter(owner, codes, bwInfos);
        size = bwInfos.size();
    }
    /*
         * Now, iterate over the results and fetch the parameter values for each BursaWolfParameters object.
         */
    final BursaWolfParameters[] parameters = new BursaWolfParameters[size];
    final Locale locale = getLocale();
    int count = 0;
    for (int i = 0; i < size; i++) {
        final BursaWolfInfo info = bwInfos.get(i);
        final GeodeticDatum datum;
        // See comment at the begining of this method.
        ensureNoCycle(BursaWolfParameters.class, code);
        try {
            datum = owner.createGeodeticDatum(String.valueOf(info.target));
        } finally {
            endOfRecursivity(BursaWolfParameters.class, code);
        }
        /*
             * Accept only Bursa-Wolf parameters between datum that use the same prime meridian.
             * This is for avoiding ambiguity about whether longitude rotation should be applied
             * before or after the datum change. This check is useless for EPSG dataset 8.9 since
             * all datum seen by this method use Greenwich. But we nevertheless perform this check
             * as a safety for future evolution or customized EPSG dataset.
             */
        if (!equalsIgnoreMetadata(meridian, datum.getPrimeMeridian())) {
            continue;
        }
        final BursaWolfParameters bwp = new BursaWolfParameters(datum, info.getDomainOfValidity(owner));
        try (ResultSet result = executeQuery("BursaWolfParameters", "SELECT PARAMETER_CODE," + " PARAMETER_VALUE," + " UOM_CODE" + " FROM [Coordinate_Operation Parameter Value]" + " WHERE COORD_OP_CODE = ?" + " AND COORD_OP_METHOD_CODE = ?", info.operation, info.method)) {
            while (result.next()) {
                BursaWolfInfo.setBursaWolfParameter(bwp, getInteger(info.operation, result, 1), getDouble(info.operation, result, 2), owner.createUnit(getString(info.operation, result, 3)), locale);
            }
        }
        if (info.isFrameRotation()) {
            // Coordinate frame rotation (9607): same as 9606,
            // except for the sign of rotation parameters.
            bwp.reverseRotation();
        }
        parameters[count++] = bwp;
    }
    return ArraysExt.resize(parameters, count);
}
Also used : Locale(java.util.Locale) ArrayList(java.util.ArrayList) ResultSet(java.sql.ResultSet) DefaultGeodeticDatum(org.apache.sis.referencing.datum.DefaultGeodeticDatum) BursaWolfParameters(org.apache.sis.referencing.datum.BursaWolfParameters)

Example 5 with BursaWolfParameters

use of org.apache.sis.referencing.datum.BursaWolfParameters 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);
}
Also used : DefaultGeodeticDatum(org.apache.sis.referencing.datum.DefaultGeodeticDatum) ParameterValueGroup(org.opengis.parameter.ParameterValueGroup) BursaWolfParameters(org.apache.sis.referencing.datum.BursaWolfParameters) Test(org.junit.Test) DependsOnMethod(org.apache.sis.test.DependsOnMethod)

Aggregations

BursaWolfParameters (org.apache.sis.referencing.datum.BursaWolfParameters)10 DefaultGeodeticDatum (org.apache.sis.referencing.datum.DefaultGeodeticDatum)5 Parameters (org.apache.sis.parameter.Parameters)4 Test (org.junit.Test)3 DefaultGeographicCRS (org.apache.sis.referencing.crs.DefaultGeographicCRS)2 DependsOnMethod (org.apache.sis.test.DependsOnMethod)2 ParameterValueGroup (org.opengis.parameter.ParameterValueGroup)2 CartesianCS (org.opengis.referencing.cs.CartesianCS)2 EllipsoidalCS (org.opengis.referencing.cs.EllipsoidalCS)2 ResultSet (java.sql.ResultSet)1 ArrayList (java.util.ArrayList)1 Locale (java.util.Locale)1 Unit (javax.measure.Unit)1 Angle (javax.measure.quantity.Angle)1 FormattableObject (org.apache.sis.io.wkt.FormattableObject)1 Formatter (org.apache.sis.io.wkt.Formatter)1 UnavailableFactoryException (org.apache.sis.referencing.factory.UnavailableFactoryException)1 LinearTransform (org.apache.sis.referencing.operation.transform.LinearTransform)1 Identifier (org.opengis.metadata.Identifier)1 GeneralParameterValue (org.opengis.parameter.GeneralParameterValue)1