Search in sources :

Example 51 with CoordinateReferenceSystem

use of org.opengis.referencing.crs.CoordinateReferenceSystem 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...
        }
    }
}
Also used : Random(java.util.Random) DefaultEllipsoid(org.apache.sis.referencing.datum.DefaultEllipsoid) CoordinateOperation(org.opengis.referencing.operation.CoordinateOperation) CoordinateReferenceSystem(org.opengis.referencing.crs.CoordinateReferenceSystem) DefaultEllipsoid(org.apache.sis.referencing.datum.DefaultEllipsoid) Ellipsoid(org.opengis.referencing.datum.Ellipsoid) Test(org.junit.Test)

Example 52 with CoordinateReferenceSystem

use of org.opengis.referencing.crs.CoordinateReferenceSystem in project sis by apache.

the class Angle method valueOf.

/**
 * Returns the angular value of the axis having the given direction.
 * This helper method is used for subclass constructors expecting a {@link DirectPosition} argument.
 *
 * @param  position  the position from which to get an angular value.
 * @param  positive  axis direction of positive values.
 * @param  negative  axis direction of negative values.
 * @return angular value in degrees.
 * @throws IllegalArgumentException if the given coordinate it not associated to a CRS,
 *         or if no axis oriented toward the given directions is found, or if that axis
 *         does not use {@linkplain Units#isAngular angular units}.
 */
static double valueOf(final DirectPosition position, final AxisDirection positive, final AxisDirection negative) {
    final CoordinateReferenceSystem crs = position.getCoordinateReferenceSystem();
    if (crs == null) {
        throw new IllegalArgumentException(Errors.format(Errors.Keys.UnspecifiedCRS));
    }
    final CoordinateSystem cs = crs.getCoordinateSystem();
    final int dimension = cs.getDimension();
    IncommensurableException cause = null;
    for (int i = 0; i < dimension; i++) {
        final CoordinateSystemAxis axis = cs.getAxis(i);
        final AxisDirection dir = axis.getDirection();
        final boolean isPositive = dir.equals(positive);
        if (isPositive || dir.equals(negative)) {
            double value = position.getOrdinate(i);
            if (!isPositive)
                value = -value;
            final Unit<?> unit = axis.getUnit();
            if (unit != Units.DEGREE)
                try {
                    value = unit.getConverterToAny(Units.DEGREE).convert(value);
                } catch (IncommensurableException e) {
                    cause = e;
                    break;
                }
            return value;
        }
    }
    throw new IllegalArgumentException(Errors.format(Errors.Keys.IllegalCRSType_1, Classes.getLeafInterfaces(crs.getClass(), CoordinateReferenceSystem.class)[0]), cause);
}
Also used : IncommensurableException(javax.measure.IncommensurableException) CoordinateSystem(org.opengis.referencing.cs.CoordinateSystem) CoordinateSystemAxis(org.opengis.referencing.cs.CoordinateSystemAxis) AxisDirection(org.opengis.referencing.cs.AxisDirection) CoordinateReferenceSystem(org.opengis.referencing.crs.CoordinateReferenceSystem)

Example 53 with CoordinateReferenceSystem

use of org.opengis.referencing.crs.CoordinateReferenceSystem in project sis by apache.

the class ServicesForMetadata method addElements.

/**
 * Initializes a horizontal, vertical and temporal extent with the values inferred from the given envelope.
 *
 * @param  envelope  the source envelope.
 * @param  target    the target extent where to store envelope information.
 * @throws TransformException if a coordinate transformation was required and failed.
 */
@Override
public void addElements(final Envelope envelope, final DefaultExtent target) throws TransformException {
    final CoordinateReferenceSystem crs = envelope.getCoordinateReferenceSystem();
    final SingleCRS horizontalCRS = CRS.getHorizontalComponent(crs);
    final VerticalCRS verticalCRS = CRS.getVerticalComponent(crs, true);
    final TemporalCRS temporalCRS = CRS.getTemporalComponent(crs);
    if (horizontalCRS == null && verticalCRS == null && temporalCRS == null) {
        throw new TransformException(dimensionNotFound(Resources.Keys.MissingSpatioTemporalDimension_1, crs));
    }
    if (horizontalCRS != null) {
        final DefaultGeographicBoundingBox extent = new DefaultGeographicBoundingBox();
        extent.setInclusion(Boolean.TRUE);
        setBounds(envelope, extent);
        target.getGeographicElements().add(extent);
    }
    if (verticalCRS != null) {
        final DefaultVerticalExtent extent = new DefaultVerticalExtent();
        setVerticalExtent(envelope, extent, crs, verticalCRS);
        target.getVerticalElements().add(extent);
    }
    if (temporalCRS != null) {
        final DefaultTemporalExtent extent = new DefaultTemporalExtent();
        setTemporalExtent(envelope, extent, crs, temporalCRS);
        target.getTemporalElements().add(extent);
    }
}
Also used : SingleCRS(org.opengis.referencing.crs.SingleCRS) DefaultTemporalCRS(org.apache.sis.referencing.crs.DefaultTemporalCRS) TemporalCRS(org.opengis.referencing.crs.TemporalCRS) DefaultGeographicBoundingBox(org.apache.sis.metadata.iso.extent.DefaultGeographicBoundingBox) DefaultVerticalExtent(org.apache.sis.metadata.iso.extent.DefaultVerticalExtent) DefaultTemporalExtent(org.apache.sis.metadata.iso.extent.DefaultTemporalExtent) VerticalCRS(org.opengis.referencing.crs.VerticalCRS) TransformException(org.opengis.referencing.operation.TransformException) CoordinateReferenceSystem(org.opengis.referencing.crs.CoordinateReferenceSystem)

Example 54 with CoordinateReferenceSystem

use of org.opengis.referencing.crs.CoordinateReferenceSystem in project sis by apache.

the class ServicesForMetadata method setBounds.

/**
 * Sets a geographic bounding box from the specified envelope.
 * If the envelope contains a CRS which is not geographic, then the bounding box will be transformed
 * to a geographic CRS (without datum shift if possible). Otherwise, the envelope is assumed already
 * in a geographic CRS using (<var>longitude</var>, <var>latitude</var>) axis order.
 *
 * @param  envelope  the source envelope.
 * @param  target    the target bounding box where to store envelope information.
 * @throws TransformException if the given envelope can not be transformed.
 */
@Override
public void setBounds(Envelope envelope, final DefaultGeographicBoundingBox target) throws TransformException {
    final CoordinateReferenceSystem crs = envelope.getCoordinateReferenceSystem();
    GeographicCRS normalizedCRS = ReferencingUtilities.toNormalizedGeographicCRS(crs);
    if (normalizedCRS == null) {
        if (crs != null) {
            normalizedCRS = CommonCRS.defaultGeographic();
        } else if (envelope.getDimension() != 2) {
            throw new TransformException(dimensionNotFound(Resources.Keys.MissingHorizontalDimension_1, crs));
        }
    }
    setGeographicExtent(envelope, target, crs, normalizedCRS);
}
Also used : TransformException(org.opengis.referencing.operation.TransformException) CoordinateReferenceSystem(org.opengis.referencing.crs.CoordinateReferenceSystem) GeographicCRS(org.opengis.referencing.crs.GeographicCRS)

Example 55 with CoordinateReferenceSystem

use of org.opengis.referencing.crs.CoordinateReferenceSystem in project sis by apache.

the class ServicesForMetadata method setBounds.

/**
 * Sets the geographic, vertical and temporal extents with the values inferred from the given envelope.
 * If the given {@code target} has more geographic or vertical extents than needed (0 or 1), then the
 * extraneous extents are removed.
 *
 * @param  envelope  the source envelope.
 * @param  target    the target spatiotemporal extent where to store envelope information.
 * @throws TransformException if no temporal component can be extracted from the given envelope.
 */
@Override
public void setBounds(final Envelope envelope, final DefaultSpatialTemporalExtent target) throws TransformException {
    final CoordinateReferenceSystem crs = envelope.getCoordinateReferenceSystem();
    final SingleCRS horizontalCRS = CRS.getHorizontalComponent(crs);
    final VerticalCRS verticalCRS = CRS.getVerticalComponent(crs, true);
    final TemporalCRS temporalCRS = CRS.getTemporalComponent(crs);
    if (horizontalCRS == null && verticalCRS == null && temporalCRS == null) {
        throw new TransformException(dimensionNotFound(Resources.Keys.MissingSpatioTemporalDimension_1, crs));
    }
    /*
         * Try to set the geographic bounding box first, because this operation may fail with a
         * TransformException while the other operations (vertical and temporal) should not fail.
         * So doing the geographic part first help us to get a "all or nothing" behavior.
         */
    DefaultGeographicBoundingBox box = null;
    boolean useExistingBox = (horizontalCRS != null);
    final Collection<GeographicExtent> spatialExtents = target.getSpatialExtent();
    final Iterator<GeographicExtent> it = spatialExtents.iterator();
    while (it.hasNext()) {
        final GeographicExtent extent = it.next();
        if (extent instanceof GeographicBoundingBox) {
            if (useExistingBox && (extent instanceof DefaultGeographicBoundingBox)) {
                box = (DefaultGeographicBoundingBox) extent;
                useExistingBox = false;
            } else {
                it.remove();
            }
        }
    }
    if (horizontalCRS != null) {
        if (box == null) {
            box = new DefaultGeographicBoundingBox();
            spatialExtents.add(box);
        }
        GeographicCRS normalizedCRS = ReferencingUtilities.toNormalizedGeographicCRS(crs);
        if (normalizedCRS == null) {
            normalizedCRS = CommonCRS.defaultGeographic();
        }
        setGeographicExtent(envelope, box, crs, normalizedCRS);
    }
    /*
         * Other dimensions (vertical and temporal).
         */
    if (verticalCRS != null) {
        VerticalExtent e = target.getVerticalExtent();
        if (!(e instanceof DefaultVerticalExtent)) {
            e = new DefaultVerticalExtent();
            target.setVerticalExtent(e);
        }
        setVerticalExtent(envelope, (DefaultVerticalExtent) e, crs, verticalCRS);
    } else {
        target.setVerticalExtent(null);
    }
    if (temporalCRS != null) {
        setTemporalExtent(envelope, target, crs, temporalCRS);
    } else {
        target.setExtent(null);
    }
}
Also used : SingleCRS(org.opengis.referencing.crs.SingleCRS) DefaultVerticalExtent(org.apache.sis.metadata.iso.extent.DefaultVerticalExtent) TransformException(org.opengis.referencing.operation.TransformException) DefaultVerticalExtent(org.apache.sis.metadata.iso.extent.DefaultVerticalExtent) VerticalExtent(org.opengis.metadata.extent.VerticalExtent) GeographicBoundingBox(org.opengis.metadata.extent.GeographicBoundingBox) DefaultGeographicBoundingBox(org.apache.sis.metadata.iso.extent.DefaultGeographicBoundingBox) GeographicExtent(org.opengis.metadata.extent.GeographicExtent) DefaultTemporalCRS(org.apache.sis.referencing.crs.DefaultTemporalCRS) TemporalCRS(org.opengis.referencing.crs.TemporalCRS) DefaultGeographicBoundingBox(org.apache.sis.metadata.iso.extent.DefaultGeographicBoundingBox) VerticalCRS(org.opengis.referencing.crs.VerticalCRS) CoordinateReferenceSystem(org.opengis.referencing.crs.CoordinateReferenceSystem) GeographicCRS(org.opengis.referencing.crs.GeographicCRS)

Aggregations

CoordinateReferenceSystem (org.opengis.referencing.crs.CoordinateReferenceSystem)210 Test (org.junit.Test)80 MathTransform (org.opengis.referencing.operation.MathTransform)32 FactoryException (org.opengis.referencing.FactoryException)25 CoordinateOperation (org.opengis.referencing.operation.CoordinateOperation)24 ReferencedEnvelope (org.geotools.geometry.jts.ReferencedEnvelope)23 Geometry (com.vividsolutions.jts.geom.Geometry)21 TransformException (org.opengis.referencing.operation.TransformException)21 DependsOnMethod (org.apache.sis.test.DependsOnMethod)19 CoordinateSystem (org.opengis.referencing.cs.CoordinateSystem)13 Geometry (org.locationtech.jts.geom.Geometry)11 FactoryException (org.opengis.util.FactoryException)11 SimpleFeature (org.opengis.feature.simple.SimpleFeature)9 DirectPosition (org.opengis.geometry.DirectPosition)9 GeographicCRS (org.opengis.referencing.crs.GeographicCRS)9 VerticalCRS (org.opengis.referencing.crs.VerticalCRS)9 CoordinateSystemAxis (org.opengis.referencing.cs.CoordinateSystemAxis)9 ArrayList (java.util.ArrayList)8 GeometryType (org.opengis.feature.type.GeometryType)8 RevFeatureType (org.locationtech.geogig.api.RevFeatureType)7