Examples with CoordinateOperation org.opengis.referencing.operation.CoordinateOperation used on opensource projects

Example 21 with CoordinateOperation

use of org.opengis.referencing.operation.CoordinateOperation in project sis by apache.

the class CoordinateOperationTest method testGeocentricTransform.

/**
 * Tests a "geographic to geocentric" conversion.
 *
 * @throws FactoryException if an error occurred while creating a test CRS.
 * @throws TransformException if an error occurred while testing a coordinate conversion.
 */
@Test
public void testGeocentricTransform() throws FactoryException, TransformException {
    final Random random = new Random(661597560);
    /*
         * Gets the math transform from WGS84 to a geocentric transform.
         */
    final Ellipsoid ellipsoid = CommonCRS.WGS84.ellipsoid();
    final CoordinateReferenceSystem sourceCRS = AbstractCRS.castOrCopy(CommonCRS.WGS84.geographic3D()).forConvention(AxesConvention.RIGHT_HANDED);
    final CoordinateReferenceSystem targetCRS = CommonCRS.WGS84.geocentric();
    final CoordinateOperation operation = opFactory.createOperation(sourceCRS, targetCRS);
    transform = operation.getMathTransform();
    final int dimension = transform.getSourceDimensions();
    assertEquals("Source dimension", 3, dimension);
    assertEquals("Target dimension", 3, transform.getTargetDimensions());
    assertSame("Inverse transform", transform, transform.inverse().inverse());
    validate();
    /*
         * Constructs an array of random points. The first 8 points
         * are initialized to know values. Other points are left random.
         */
    final double[] cartesianDistance = new double[4];
    final double[] orthodromicDistance = new double[4];
    // Must be divisible by 3.
    final double[] array0 = new double[900];
    for (int i = 0; i < array0.length; i++) {
        final int range;
        switch(i % 3) {
            // Longitude
            case 0:
                range = 360;
                break;
            // Latitidue
            case 1:
                range = 180;
                break;
            // Altitude
            case 2:
                range = 10000;
                break;
            // Should not happen
            default:
                range = 0;
                break;
        }
        array0[i] = random.nextDouble() * range - (range / 2);
    }
    // 24°N 35°E 8km
    array0[0] = 35.0;
    // 24°N 35°E 8km
    array0[1] = 24.0;
    // 24°N 35°E 8km
    array0[2] = 8000;
    // … about 80 km away
    array0[3] = 34.8;
    // … about 80 km away
    array0[4] = 24.7;
    // … about 80 km away
    array0[5] = 5000;
    cartesianDistance[0] = 80284.00;
    // Not really exact.
    orthodromicDistance[0] = 80302.99;
    array0[6] = 0;
    array0[7] = 0.0;
    array0[8] = 0;
    // Antipodes; distance should be 2*6378.137 km
    array0[9] = 180;
    // Antipodes; distance should be 2*6378.137 km
    array0[10] = 0.0;
    // Antipodes; distance should be 2*6378.137 km
    array0[11] = 0;
    cartesianDistance[1] = ellipsoid.getSemiMajorAxis() * 2;
    orthodromicDistance[1] = ellipsoid.getSemiMajorAxis() * PI;
    array0[12] = 0;
    array0[13] = -90;
    array0[14] = 0;
    // Antipodes; distance should be 2*6356.752 km
    array0[15] = 180;
    // Antipodes; distance should be 2*6356.752 km
    array0[16] = +90;
    // Antipodes; distance should be 2*6356.752 km
    array0[17] = 0;
    cartesianDistance[2] = ellipsoid.getSemiMinorAxis() * 2;
    orthodromicDistance[2] = 20003931.46;
    array0[18] = 95;
    array0[19] = -38;
    array0[20] = 0;
    // Antipodes
    array0[21] = -85;
    // Antipodes
    array0[22] = +38;
    // Antipodes
    array0[23] = 0;
    cartesianDistance[3] = 12740147.19;
    orthodromicDistance[3] = 20003867.86;
    /*
         * Transforms all points, and then inverse transform them. The resulting
         * array2 should be equal to array0 except for rounding errors. We tolerate
         * maximal error of 0.1 second in longitude or latitude and 1 cm in height.
         */
    final double[] array1 = new double[array0.length];
    final double[] array2 = new double[array0.length];
    transform.transform(array0, 0, array1, 0, array0.length / dimension);
    transform.inverse().transform(array1, 0, array2, 0, array1.length / dimension);
    for (int i = 0; i < array0.length; ) {
        assertEquals("Longitude", array2[i], array0[i], 0.1 / 3600);
        i++;
        assertEquals("Latitude", array2[i], array0[i], 0.1 / 3600);
        i++;
        assertEquals("Height", array2[i], array0[i], 0.01);
        i++;
    }
    /*
         * Compares the distances between "special" points with expected distances.
         * This tests the ellipsoid orthodromic distance computation as well.
         * We require a precision of 10 centimetres.
         */
    for (int i = 0; i < array0.length / 6; i++) {
        final int base = i * 6;
        final double cartesian = MathFunctions.magnitude(array1[base + 0] - array1[base + 3], array1[base + 1] - array1[base + 4], array1[base + 2] - array1[base + 5]);
        if (i < cartesianDistance.length) {
            assertEquals("Cartesian distance", cartesianDistance[i], cartesian, 0.1);
        }
        /*
             * Compares with orthodromic distance. Distance is computed using an ellipsoid
             * at the maximal altitude (i.e. the length of semi-major axis is increased to
             * fit the maximal altitude).
             */
        try {
            final double altitude = max(array0[base + 2], array0[base + 5]);
            final DefaultEllipsoid ellip = DefaultEllipsoid.createFlattenedSphere(Collections.singletonMap(Ellipsoid.NAME_KEY, "Temporary"), ellipsoid.getSemiMajorAxis() + altitude, ellipsoid.getInverseFlattening(), ellipsoid.getAxisUnit());
            double orthodromic = ellip.orthodromicDistance(array0[base + 0], array0[base + 1], array0[base + 3], array0[base + 4]);
            orthodromic = hypot(orthodromic, array0[base + 2] - array0[base + 5]);
            if (i < orthodromicDistance.length) {
                assertEquals("Orthodromic distance", orthodromicDistance[i], orthodromic, 0.1);
            }
            assertTrue("Distance consistency", cartesian <= orthodromic);
        } catch (ArithmeticException exception) {
        // Orthodromic distance computation didn't converge. Ignore...
        }
    }
}

Example 22 with CoordinateOperation

use of org.opengis.referencing.operation.CoordinateOperation in project sis by apache.

the class PositionalAccuracyConstant method getLinearAccuracy.

/**
 * Convenience method returning the accuracy in meters for the specified operation.
 * This method tries each of the following procedures and returns the first successful one:
 *
 * <ul>
 *   <li>If at least one {@link QuantitativeResult} is found with a linear unit, then the largest
 *       accuracy estimate is converted to {@linkplain Units#METRE metres} and returned.</li>
 *   <li>Otherwise, if the operation is a {@link Conversion}, then returns 0 since a conversion
 *       is by definition accurate up to rounding errors.</li>
 *   <li>Otherwise, if the operation is a {@link Transformation}, then checks if the datum shift
 *       were applied with the help of Bursa-Wolf parameters. This procedure looks for SIS-specific
 *       {@link #DATUM_SHIFT_APPLIED} and {@link #DATUM_SHIFT_OMITTED DATUM_SHIFT_OMITTED} constants.</li>
 *   <li>Otherwise, if the operation is a {@link ConcatenatedOperation}, returns the sum of the accuracy
 *       of all components. This is a conservative scenario where we assume that errors cumulate linearly.
 *       Note that this is not necessarily the "worst case" scenario since the accuracy could be worst
 *       if the math transforms are highly non-linear.</li>
 * </ul>
 *
 * If the above is modified, please update {@code AbstractCoordinateOperation.getLinearAccuracy()} javadoc.
 *
 * @param  operation  the operation to inspect for accuracy.
 * @return the accuracy estimate (always in meters), or NaN if unknown.
 *
 * @see org.apache.sis.referencing.operation.AbstractCoordinateOperation#getLinearAccuracy()
 */
public static double getLinearAccuracy(final CoordinateOperation operation) {
    double accuracy = Double.NaN;
    final Collection<PositionalAccuracy> accuracies = operation.getCoordinateOperationAccuracy();
    for (final PositionalAccuracy metadata : accuracies) {
        for (final Result result : metadata.getResults()) {
            if (result instanceof QuantitativeResult) {
                final QuantitativeResult quantity = (QuantitativeResult) result;
                final Collection<? extends Record> records = quantity.getValues();
                if (records != null) {
                    final Unit<?> unit = quantity.getValueUnit();
                    if (Units.isLinear(unit)) {
                        final Unit<Length> unitOfLength = unit.asType(Length.class);
                        for (final Record record : records) {
                            for (final Object value : record.getAttributes().values()) {
                                if (value instanceof Number) {
                                    double v = ((Number) value).doubleValue();
                                    v = unitOfLength.getConverterTo(Units.METRE).convert(v);
                                    if (v >= 0 && !(v <= accuracy)) {
                                        // '!' is for replacing the NaN value.
                                        accuracy = v;
                                    }
                                }
                            }
                        }
                    }
                }
            }
        }
    }
    if (Double.isNaN(accuracy)) {
        /*
             * No quantitative (linear) accuracy were found. If the coordinate operation is actually
             * a conversion, the accuracy is up to rounding error (i.e. conceptually 0) by definition.
             */
        if (operation instanceof Conversion) {
            return 0;
        }
        /*
             * If the coordinate operation is actually a transformation, checks if Bursa-Wolf parameters
             * were available for the datum shift. This is SIS-specific. See field javadoc for a rational
             * about the return values chosen.
             */
        if (operation instanceof Transformation) {
            if (accuracies.contains(DATUM_SHIFT_APPLIED)) {
                return DATUM_SHIFT_ACCURACY;
            }
            if (accuracies.contains(DATUM_SHIFT_OMITTED)) {
                return UNKNOWN_ACCURACY;
            }
        }
        /*
             * If the coordinate operation is a compound of other coordinate operations, returns the sum of their accuracy,
             * skipping unknown ones. Making the sum is a conservative approach (not exactly the "worst case" scenario,
             * since it could be worst if the transforms are highly non-linear).
             */
        if (operation instanceof ConcatenatedOperation) {
            for (final CoordinateOperation op : ((ConcatenatedOperation) operation).getOperations()) {
                final double candidate = Math.abs(getLinearAccuracy(op));
                if (!Double.isNaN(candidate)) {
                    if (Double.isNaN(accuracy)) {
                        accuracy = candidate;
                    } else {
                        accuracy += candidate;
                    }
                }
            }
        }
    }
    return accuracy;
}

Example 23 with CoordinateOperation

use of org.opengis.referencing.operation.CoordinateOperation in project sis by apache.

the class Formatter method appendComplement.

/**
 * Appends the optional complementary attributes common to many {@link IdentifiedObject} subtypes.
 * Those attributes are {@code ANCHOR}, {@code SCOPE}, {@code AREA}, {@code BBOX}, {@code VERTICALEXTENT},
 * {@code TIMEEXTENT}, {@code ID} (previously known as {@code AUTHORITY}) and {@code REMARKS},
 * and have a special treatment: they are written by {@link #append(FormattableObject)}
 * after the {@code formatTo(Formatter)} method returned.
 *
 * <p>The {@code ID[<name>,<code>,…]} element is normally written only for the root element
 * (unless the convention is {@code INTERNAL}), but there is various exceptions to this rule.
 * If formatted, the {@code ID} element will be by default on the same line than the enclosing
 * element (e.g. {@code SPHEROID["Clarke 1866", …, ID["EPSG", 7008]]}). Other example:</p>
 *
 * {@preformat text
 *   PROJCS["NAD27 / Idaho Central",
 *     GEOGCS[...etc...],
 *     ...etc...
 *     ID["EPSG", 26769]]
 * }
 *
 * For non-internal conventions, all elements other than {@code ID[…]} are formatted
 * only for {@link CoordinateOperation} and root {@link ReferenceSystem} instances,
 * with an exception for remarks of {@code ReferenceSystem} embedded inside {@code CoordinateOperation}.
 * Those restrictions are our interpretation of the following ISO 19162 requirement:
 *
 * <blockquote>(…snip…) {@code <scope extent identifier remark>} is a collection of four optional attributes
 * which may be applied to a coordinate reference system, a coordinate operation or a boundCRS. (…snip…)
 * Identifier (…snip…) may also be utilised for components of these objects although this is not recommended
 * except for coordinate operation methods (including map projections) and parameters. (…snip…)
 * A {@code <remark>} can be included within the descriptions of source and target CRS embedded within
 * a coordinate transformation as well as within the coordinate transformation itself.</blockquote>
 */
@SuppressWarnings("null")
private void appendComplement(final IdentifiedObject object, final FormattableObject parent, final FormattableObject gp) {
    isComplement = true;
    // Whether to format ID[…] elements.
    final boolean showIDs;
    // Whether we shall limit to a single ID[…] element.
    final boolean filterID;
    // Whether to format any element other than ID[…] and Remarks[…].
    final boolean showOthers;
    // Whether to format Remarks[…].
    final boolean showRemarks;
    if (convention == Convention.INTERNAL) {
        showIDs = true;
        filterID = false;
        showOthers = true;
        showRemarks = true;
    } else {
        /*
             * Except for the special cases of OperationMethod and Parameters, ISO 19162 recommends to format the
             * ID only for the root element.  But Apache SIS adds an other exception to this rule by handling the
             * components of CompoundCRS as if they were root elements. The reason is that users often create their
             * own CompoundCRS from standard components, for example by adding a time axis to some standard CRS like
             * "WGS84". The resulting CompoundCRS usually have no identifier. Then the users often need to extract a
             * particular component of a CompoundCRS, most often the horizontal part, and will need its identifier
             * for example in a Web Map Service (WMS). Those ID are lost if we do not format them here.
             */
        if (parent == null || parent instanceof CompoundCRS) {
            showIDs = true;
        } else if (gp instanceof CoordinateOperation && !(parent instanceof IdentifiedObject)) {
            // "SourceCRS[…]" and "TargetCRS[…]" sub-elements in CoordinateOperation.
            showIDs = true;
        } else if (convention == Convention.WKT2_SIMPLIFIED) {
            showIDs = false;
        } else {
            showIDs = (object instanceof OperationMethod) || (object instanceof GeneralParameterDescriptor);
        }
        if (convention.majorVersion() == 1) {
            filterID = true;
            showOthers = false;
            showRemarks = false;
        } else {
            filterID = (parent != null);
            if (object instanceof CoordinateOperation) {
                showOthers = !(parent instanceof ConcatenatedOperation);
                showRemarks = showOthers;
            } else if (object instanceof ReferenceSystem) {
                showOthers = (parent == null);
                showRemarks = (parent == null) || (gp instanceof CoordinateOperation);
            } else {
                // Mandated by ISO 19162.
                showOthers = false;
                showRemarks = false;
            }
        }
    }
    if (showOthers) {
        appendForSubtypes(object);
    }
    if (showIDs) {
        Collection<ReferenceIdentifier> identifiers = object.getIdentifiers();
        if (identifiers != null) {
            // Paranoiac check
            if (filterID) {
                for (final ReferenceIdentifier id : identifiers) {
                    if (Citations.identifierMatches(authority, id.getAuthority())) {
                        identifiers = Collections.singleton(id);
                        break;
                    }
                }
            }
            for (ReferenceIdentifier id : identifiers) {
                if (!(id instanceof FormattableObject)) {
                    id = ImmutableIdentifier.castOrCopy(id);
                }
                append((FormattableObject) id);
                if (filterID)
                    break;
            }
        }
    }
    if (showRemarks) {
        appendOnNewLine(WKTKeywords.Remark, object.getRemarks(), ElementKind.REMARKS);
    }
    isComplement = false;
}

Example 24 with CoordinateOperation

use of org.opengis.referencing.operation.CoordinateOperation in project sis by apache.

the class DefaultConcatenatedOperation method initialize.

/**
 * Initializes the {@link #sourceCRS}, {@link #targetCRS} and {@link #operations} fields.
 * If the source or target CRS is already non-null (which may happen on JAXB unmarshalling),
 * leaves that CRS unchanged.
 *
 * @param  properties   the properties specified at construction time, or {@code null} if unknown.
 * @param  operations   the operations to concatenate.
 * @param  mtFactory    the math transform factory to use, or {@code null} for not performing concatenation.
 * @throws FactoryException if the factory can not concatenate the math transforms.
 */
private void initialize(final Map<String, ?> properties, final CoordinateOperation[] operations, final MathTransformFactory mtFactory) throws FactoryException {
    final List<CoordinateOperation> flattened = new ArrayList<>(operations.length);
    final CoordinateReferenceSystem crs = initialize(properties, operations, flattened, mtFactory, (sourceCRS == null), (coordinateOperationAccuracy == null), (domainOfValidity == null));
    if (targetCRS == null) {
        targetCRS = crs;
    }
    /*
         * At this point we should have flattened.size() >= 2, except if some operations
         * were omitted because their associated math transform were identity operation.
         */
    this.operations = UnmodifiableArrayList.wrap(flattened.toArray(new SingleOperation[flattened.size()]));
}

Example 25 with CoordinateOperation

use of org.opengis.referencing.operation.CoordinateOperation in project sis by apache.

the class SubOperationInfo method create.

/**
 * Searches in given list of source components for an operation capable to convert or transform coordinates
 * to the given target CRS. If no such operation can be found, then this method returns {@code null}.
 *
 * @param  caller        the object which is inferring a coordinate operation.
 * @param  sourceIsUsed  flags for keeping trace of which source has been used.
 * @param  sources       all components of the source CRS.
 * @param  target        one component of the target CRS.
 * @return information about a coordinate operation from a source CRS to the given target CRS, or {@code null}.
 * @throws FactoryException if an error occurred while grabbing a coordinate operation.
 */
static SubOperationInfo create(final CoordinateOperationFinder caller, final boolean[] sourceIsUsed, final List<? extends SingleCRS> sources, final SingleCRS target) throws FactoryException {
    OperationNotFoundException failure = null;
    final Class<?> targetType = type(target);
    for (final Class<?>[] sourceTypes : COMPATIBLE_TYPES) {
        if (sourceTypes[0].isAssignableFrom(targetType)) {
            for (final Class<?> sourceType : sourceTypes) {
                int startAtDimension;
                int endAtDimension = 0;
                for (int i = 0; i < sourceIsUsed.length; i++) {
                    final SingleCRS source = sources.get(i);
                    startAtDimension = endAtDimension;
                    endAtDimension += source.getCoordinateSystem().getDimension();
                    if (!sourceIsUsed[i] && sourceType.isAssignableFrom(type(source))) {
                        final CoordinateOperation operation;
                        try {
                            operation = caller.createOperation(source, target);
                        } catch (OperationNotFoundException exception) {
                            if (failure == null) {
                                failure = exception;
                            } else {
                                failure.addSuppressed(exception);
                            }
                            continue;
                        }
                        /*
                             * Found an operation.  Exclude the source component from the list because each source
                             * should be used at most once by SourceComponent. Note that the same source may still
                             * be used again in another context if that source is also an interpolation CRS.
                             *
                             * EXAMPLE: consider a coordinate operation from (GeodeticCRS₁, VerticalCRS₁) source
                             * to (GeodeticCRS₂, VerticalCRS₂) target.  The source GeodeticCRS₁ should be mapped
                             * to exactly one target component (which is GeodeticCRS₂)  and  VerticalCRS₁ mapped
                             * to VerticalCRS₂.  But the operation on vertical coordinates may need GeodeticCRS₁
                             * for doing its work, so GeodeticCRS₁ is needed twice.  However when needed for the
                             * vertical coordinate operation, the GeodeticCRS₁ is used as an "interpolation CRS".
                             * Interpolation CRS are handled in other code paths; it is not the business of this
                             * SourceComponent class to care about them. From the point of view of this class,
                             * GeodeticCRS₁ is used only once.
                             */
                        sourceIsUsed[i] = true;
                        if (failure != null) {
                            Logging.recoverableException(Logging.getLogger(Loggers.COORDINATE_OPERATION), CoordinateOperationFinder.class, "decompose", failure);
                        }
                        return new SubOperationInfo(operation, startAtDimension, endAtDimension);
                    }
                }
            }
        }
    }
    if (failure != null) {
        throw failure;
    }
    return null;
}