Search in sources :

Example 11 with CartesianCS

use of org.opengis.referencing.cs.CartesianCS in project sis by apache.

the class CRSTest method testSuggestCommonTarget.

/**
 * Tests {@link CRS#suggestCommonTarget(GeographicBoundingBox, CoordinateReferenceSystem...)}.
 *
 * @since 0.8
 */
@Test
public void testSuggestCommonTarget() {
    /*
         * Prepare 4 CRS with different datum (so we can more easily differentiate them in the assertions) and
         * different domain of validity. CRS[1] is given a domain large enough for all CRS except the last one.
         */
    final Map<String, Object> properties = new HashMap<>(4);
    final CartesianCS cs = (CartesianCS) StandardDefinitions.createCoordinateSystem(Constants.EPSG_PROJECTED_CS);
    final ProjectedCRS[] crs = new ProjectedCRS[4];
    for (int i = 0; i < crs.length; i++) {
        final CommonCRS baseCRS;
        final double ymin, ymax;
        switch(i) {
            case 0:
                baseCRS = CommonCRS.WGS84;
                ymin = 2;
                ymax = 4;
                break;
            case 1:
                baseCRS = CommonCRS.WGS72;
                ymin = 1;
                ymax = 4;
                break;
            case 2:
                baseCRS = CommonCRS.SPHERE;
                ymin = 2;
                ymax = 3;
                break;
            case 3:
                baseCRS = CommonCRS.NAD27;
                ymin = 3;
                ymax = 5;
                break;
            default:
                throw new AssertionError(i);
        }
        properties.put(DefaultProjectedCRS.NAME_KEY, "CRS #" + i);
        properties.put(DefaultProjectedCRS.DOMAIN_OF_VALIDITY_KEY, new DefaultExtent(null, new DefaultGeographicBoundingBox(-1, +1, ymin, ymax), null, null));
        crs[i] = new DefaultProjectedCRS(properties, baseCRS.geographic(), HardCodedConversions.MERCATOR, cs);
    }
    // Exclude the last CRS only.
    final ProjectedCRS[] overlappingCRS = Arrays.copyOf(crs, 3);
    /*
         * Test between the 3 overlapping CRS without region of interest. We expect the CRS having a domain
         * of validity large enough for all CRS; this is the second CRS created in above 'switch' statement.
         */
    assertSame("Expected CRS with widest domain of validity.", crs[1], CRS.suggestCommonTarget(null, overlappingCRS));
    /*
         * If we specify a smaller region of interest, we should get the CRS having the smallest domain of validity that
         * cover the ROI. Following lines gradually increase the ROI size and verify that we get CRS for larger domain.
         */
    final DefaultGeographicBoundingBox regionOfInterest = new DefaultGeographicBoundingBox(-1, +1, 2.1, 2.9);
    assertSame("Expected best fit for [2.1 … 2.9]°N", crs[2], CRS.suggestCommonTarget(regionOfInterest, overlappingCRS));
    regionOfInterest.setNorthBoundLatitude(3.1);
    assertSame("Expected best fit for [2.1 … 3.1]°N", crs[0], CRS.suggestCommonTarget(regionOfInterest, overlappingCRS));
    regionOfInterest.setSouthBoundLatitude(1.9);
    assertSame("Expected best fit for [1.9 … 3.1]°N", crs[1], CRS.suggestCommonTarget(regionOfInterest, overlappingCRS));
    /*
         * All above tests returned one of the CRS in the given array. Test now a case where none of those CRS
         * have a domain of validity wide enough, so suggestCommonTarget(…) need to search among the base CRS.
         */
    assertSame("Expected a GeodeticCRS since none of the ProjectedCRS have a domain of validity wide enough.", crs[0].getBaseCRS(), CRS.suggestCommonTarget(null, crs));
    /*
         * With the same domain of validity than above, suggestCommonTarget(…) should not need to fallback on the
         * base CRS anymore.
         */
    assertSame("Expected best fit for [1.9 … 3.1]°N", crs[1], CRS.suggestCommonTarget(regionOfInterest, crs));
}
Also used : CartesianCS(org.opengis.referencing.cs.CartesianCS) HashMap(java.util.HashMap) DefaultGeographicBoundingBox(org.apache.sis.metadata.iso.extent.DefaultGeographicBoundingBox) ProjectedCRS(org.opengis.referencing.crs.ProjectedCRS) DefaultProjectedCRS(org.apache.sis.referencing.crs.DefaultProjectedCRS) DefaultExtent(org.apache.sis.metadata.iso.extent.DefaultExtent) DefaultProjectedCRS(org.apache.sis.referencing.crs.DefaultProjectedCRS) Test(org.junit.Test)

Example 12 with CartesianCS

use of org.opengis.referencing.cs.CartesianCS in project sis by apache.

the class CoordinateSystemTransform method create.

/**
 * Implementation of {@link DefaultMathTransformFactory#createCoordinateSystemChange(CoordinateSystem,
 * CoordinateSystem, Ellipsoid)}, defined here for reducing the {@code DefaultMathTransformFactory}
 * weight in the common case where the conversions handled by this class are not needed.
 */
static MathTransform create(final MathTransformFactory factory, final CoordinateSystem source, final CoordinateSystem target) throws FactoryException {
    int passthrough = 0;
    CoordinateSystemTransform kernel = null;
    if (source instanceof CartesianCS) {
        if (target instanceof SphericalCS) {
            kernel = CartesianToSpherical.INSTANCE;
        } else if (target instanceof PolarCS) {
            kernel = CartesianToPolar.INSTANCE;
        } else if (target instanceof CylindricalCS) {
            kernel = CartesianToPolar.INSTANCE;
            passthrough = 1;
        }
    } else if (target instanceof CartesianCS) {
        if (source instanceof SphericalCS) {
            kernel = SphericalToCartesian.INSTANCE;
        } else if (source instanceof PolarCS) {
            kernel = PolarToCartesian.INSTANCE;
        } else if (source instanceof CylindricalCS) {
            kernel = PolarToCartesian.INSTANCE;
            passthrough = 1;
        }
    }
    Exception cause = null;
    try {
        if (kernel == null) {
            return factory.createAffineTransform(CoordinateSystems.swapAndScaleAxes(source, target));
        } else if (source.getDimension() == kernel.getSourceDimensions() + passthrough && target.getDimension() == kernel.getTargetDimensions() + passthrough) {
            final MathTransform tr = (passthrough == 0) ? kernel.completeTransform(factory) : kernel.passthrough(factory);
            final MathTransform before = factory.createAffineTransform(CoordinateSystems.swapAndScaleAxes(source, CoordinateSystems.replaceAxes(source, AxesConvention.NORMALIZED)));
            final MathTransform after = factory.createAffineTransform(CoordinateSystems.swapAndScaleAxes(CoordinateSystems.replaceAxes(target, AxesConvention.NORMALIZED), target));
            return factory.createConcatenatedTransform(before, factory.createConcatenatedTransform(tr, after));
        }
    } catch (IllegalArgumentException | IncommensurableException e) {
        cause = e;
    }
    throw new OperationNotFoundException(Resources.format(Resources.Keys.CoordinateOperationNotFound_2, WKTUtilities.toType(CoordinateSystem.class, source.getClass()), WKTUtilities.toType(CoordinateSystem.class, target.getClass())), cause);
}
Also used : CartesianCS(org.opengis.referencing.cs.CartesianCS) SphericalCS(org.opengis.referencing.cs.SphericalCS) PolarCS(org.opengis.referencing.cs.PolarCS) CylindricalCS(org.opengis.referencing.cs.CylindricalCS) OperationNotFoundException(org.opengis.referencing.operation.OperationNotFoundException) MathTransform(org.opengis.referencing.operation.MathTransform) IncommensurableException(javax.measure.IncommensurableException) OperationNotFoundException(org.opengis.referencing.operation.OperationNotFoundException) IncommensurableException(javax.measure.IncommensurableException) FactoryException(org.opengis.util.FactoryException)

Example 13 with CartesianCS

use of org.opengis.referencing.cs.CartesianCS in project sis by apache.

the class DefaultMathTransformFactory method createCoordinateSystemChange.

/**
 * Creates a math transform that represent a change of coordinate system. If exactly one argument is
 * an {@linkplain org.apache.sis.referencing.cs.DefaultEllipsoidalCS ellipsoidal coordinate systems},
 * then the {@code ellipsoid} argument is mandatory. In all other cases (including the case where both
 * coordinate systems are ellipsoidal), the ellipsoid argument is ignored and can be {@code null}.
 *
 * <div class="note"><b>Design note:</b>
 * this method does not accept separated ellipsoid arguments for {@code source} and {@code target} because
 * this method should not be used for datum shifts. If the two given coordinate systems are ellipsoidal,
 * then they are assumed to use the same ellipsoid. If different ellipsoids are desired, then a
 * {@linkplain #createParameterizedTransform parameterized transform} like <cite>"Molodensky"</cite>,
 * <cite>"Geocentric translations"</cite>, <cite>"Coordinate Frame Rotation"</cite> or
 * <cite>"Position Vector transformation"</cite> should be used instead.</div>
 *
 * @param  source     the source coordinate system.
 * @param  target     the target coordinate system.
 * @param  ellipsoid  the ellipsoid of {@code EllipsoidalCS}, or {@code null} if none.
 * @return a conversion from the given source to the given target coordinate system.
 * @throws FactoryException if the conversion can not be created.
 *
 * @since 0.8
 */
public MathTransform createCoordinateSystemChange(final CoordinateSystem source, final CoordinateSystem target, final Ellipsoid ellipsoid) throws FactoryException {
    ArgumentChecks.ensureNonNull("source", source);
    ArgumentChecks.ensureNonNull("target", target);
    if (ellipsoid != null) {
        final boolean isEllipsoidalSource = (source instanceof EllipsoidalCS);
        if (isEllipsoidalSource != (target instanceof EllipsoidalCS)) {
            /*
                 * For now we support only conversion between EllipsoidalCS and CartesianCS.
                 * But future Apache SIS versions could add support for conversions between
                 * EllipsoidalCS and SphericalCS or other coordinate systems.
                 */
            if ((isEllipsoidalSource ? target : source) instanceof CartesianCS) {
                final Context context = new Context();
                final EllipsoidalCS cs;
                final String operation;
                if (isEllipsoidalSource) {
                    operation = GeographicToGeocentric.NAME;
                    context.setSource(cs = (EllipsoidalCS) source, ellipsoid);
                    context.setTarget(target);
                } else {
                    operation = GeocentricToGeographic.NAME;
                    context.setSource(source);
                    context.setTarget(cs = (EllipsoidalCS) target, ellipsoid);
                }
                final ParameterValueGroup pg = getDefaultParameters(operation);
                // Apache SIS specific parameter.
                if (cs.getDimension() < 3)
                    pg.parameter("dim").setValue(2);
                return createParameterizedTransform(pg, context);
            }
        }
    }
    return CoordinateSystemTransform.create(this, source, target);
// No need to use unique(…) here.
}
Also used : CartesianCS(org.opengis.referencing.cs.CartesianCS) DefaultParameterValueGroup(org.apache.sis.parameter.DefaultParameterValueGroup) ParameterValueGroup(org.opengis.parameter.ParameterValueGroup) EllipsoidalCS(org.opengis.referencing.cs.EllipsoidalCS)

Example 14 with CartesianCS

use of org.opengis.referencing.cs.CartesianCS in project sis by apache.

the class Proj4 method definition.

/**
 * Infers a {@literal Proj.4} definition from the given projected, geographic or geocentric coordinate reference system.
 * This method does not need the Proj.4 native library; it can be used in a pure Java application.
 * However the returned definition string may differ depending on whether the Proj.4 library is available or not.
 *
 * @param  crs  the coordinate reference system for which to create a Proj.4 definition.
 * @return the definition of the given CRS in a Proj.4 format.
 * @throws FactoryException if the Proj.4 definition string can not be created from the given CRS.
 */
public static String definition(final CoordinateReferenceSystem crs) throws FactoryException {
    ArgumentChecks.ensureNonNull("crs", crs);
    /*
         * If the given CRS object is associated to a Proj.4 structure, let Proj.4 formats itself
         * the definition string. Note that this operation may fail if there is no Proj.4 library
         * in the current system, or no JNI bindings to that library.
         */
    try {
        for (final Identifier id : crs.getIdentifiers()) {
            if (id instanceof PJ) {
                return ((PJ) id).getCode();
            }
        }
    } catch (UnsatisfiedLinkError e) {
        // Thrown the first time that we try to use the library.
        Logging.unexpectedException(Logging.getLogger(Modules.GDAL), Proj4.class, "definition", e);
    } catch (NoClassDefFoundError e) {
        // Thrown on all attempts after the first one.
        Logging.recoverableException(Logging.getLogger(Modules.GDAL), Proj4.class, "definition", e);
    }
    /*
         * If we found no Proj.4 structure, formats the definition string ourself. The string may differ from
         * what Proj.4 would have given. In particular, we do not provide "+init=" or "+datum=" parameter.
         * But the definition should still be semantically equivalent.
         */
    final String method;
    final GeodeticDatum datum;
    final ParameterValueGroup parameters;
    final CoordinateSystem cs = crs.getCoordinateSystem();
    if (crs instanceof GeodeticCRS) {
        if (cs instanceof EllipsoidalCS) {
            method = "latlon";
        } else if (cs instanceof CartesianCS) {
            method = "geocent";
        } else {
            throw new FactoryException(Errors.format(Errors.Keys.UnsupportedCoordinateSystem_1, cs.getClass()));
        }
        datum = ((GeodeticCRS) crs).getDatum();
        parameters = null;
    } else if (crs instanceof ProjectedCRS) {
        Projection c = ((ProjectedCRS) crs).getConversionFromBase();
        datum = ((ProjectedCRS) crs).getDatum();
        method = name(c.getMethod());
        parameters = c.getParameterValues();
    } else {
        throw new FactoryException(Errors.format(Errors.Keys.UnsupportedType_1, crs.getClass()));
    }
    /*
         * Append the map projection parameters. Those parameters may include axis lengths (a and b),
         * but not necessarily. If axis lengths are specified, then we will ignore the Ellipsoid instance
         * associated to the CRS.
         */
    final StringBuilder definition = new StringBuilder(100);
    definition.append(Proj4Factory.PROJ_PARAM).append(method);
    boolean hasSemiMajor = false;
    boolean hasSemiMinor = false;
    if (parameters != null) {
        definition.append(Proj4Factory.STANDARD_OPTIONS);
        for (final GeneralParameterValue parameter : parameters.values()) {
            if (parameter instanceof ParameterValue<?>) {
                final ParameterValue<?> pv = (ParameterValue<?>) parameter;
                final Object value;
                Unit<?> unit = pv.getUnit();
                if (unit != null) {
                    unit = Units.isAngular(unit) ? Units.DEGREE : unit.getSystemUnit();
                    // Always in metres or degrees.
                    value = pv.doubleValue(unit);
                } else {
                    value = pv.getValue();
                    if (value == null) {
                        continue;
                    }
                }
                final String pn = name(parameter.getDescriptor());
                hasSemiMajor |= pn.equals("a");
                hasSemiMinor |= pn.equals("b");
                definition.append(" +").append(pn).append('=').append(value);
            }
        }
    }
    /*
         * Append datum information: axis lengths if they were not part of the parameters, then prime meridian.
         */
    final Ellipsoid ellipsoid = datum.getEllipsoid();
    if (!hasSemiMajor)
        definition.append(" +a=").append(ellipsoid.getSemiMajorAxis());
    if (!hasSemiMinor)
        definition.append(" +b=").append(ellipsoid.getSemiMinorAxis());
    final PrimeMeridian pm = datum.getPrimeMeridian();
    if (pm != null) {
        double lon = pm.getGreenwichLongitude();
        final Unit<Angle> unit = pm.getAngularUnit();
        if (unit != null) {
            lon = unit.getConverterTo(Units.DEGREE).convert(lon);
        }
        definition.append(" +pm=").append(lon);
    }
    /*
         * Appends axis directions. This method always format a vertical direction (up or down)
         * even if the coordinate system is two-dimensional, because Proj.4 seems to require it.
         * Also extract axis units in the process.
         */
    // Horizontal at index 0, vertical at index 1.
    final Unit<?>[] units = new Unit<?>[2];
    boolean validCS = true;
    definition.append(' ').append(Proj4Factory.AXIS_ORDER_PARAM);
    final int dimension = Math.min(cs.getDimension(), 3);
    boolean hasVertical = false;
    for (int i = 0; i < dimension; i++) {
        final CoordinateSystemAxis axis = cs.getAxis(i);
        final AxisDirection dir = axis.getDirection();
        int unitIndex = 0;
        if (!AxisDirections.isCardinal(dir)) {
            if (!AxisDirections.isVertical(dir)) {
                throw new FactoryException(Errors.format(Errors.Keys.UnsupportedAxisDirection_1, dir));
            }
            hasVertical = true;
            unitIndex = 1;
        }
        final Unit<?> old = units[unitIndex];
        units[unitIndex] = axis.getUnit();
        validCS &= (old == null || old.equals(units[unitIndex]));
        definition.appendCodePoint(Character.toLowerCase(dir.name().codePointAt(0)));
    }
    if (!hasVertical && dimension < 3) {
        // Add a UP direction if not already present.
        definition.append('u');
    }
    /*
         * Append units of measurement, then verify the coordinate system validity.
         */
    for (int i = 0; i < units.length; i++) {
        final Unit<?> unit = units[i];
        if (unit != null && !unit.equals(Units.DEGREE) && !unit.equals(Units.METRE)) {
            validCS &= Units.isLinear(unit);
            definition.append(" +");
            // "+vto_meter" parameter.
            if (i == 1)
                definition.append('v');
            definition.append("to_meter=").append(Units.toStandardUnit(unit));
        }
    }
    /*
         * Append the "+towgs84" element if any. This is the last piece of information.
         * Note that the use of a "+towgs84" parameter is an "early binding" approach,
         * which is usually not recommended. But Proj4 works that way.
         */
    if (validCS) {
        if (datum instanceof DefaultGeodeticDatum) {
            for (final BursaWolfParameters bwp : ((DefaultGeodeticDatum) datum).getBursaWolfParameters()) {
                if (Utilities.equalsIgnoreMetadata(CommonCRS.WGS84.datum(), bwp.getTargetDatum())) {
                    definition.append(" +towgs84=").append(bwp.tX).append(',').append(bwp.tY).append(',').append(bwp.tZ);
                    if (!bwp.isTranslation()) {
                        definition.append(',').append(bwp.rX).append(',').append(bwp.rY).append(',').append(bwp.rZ).append(',').append(bwp.dS);
                    }
                    break;
                }
            }
        }
        return definition.toString();
    }
    /*
         * If we reach this point, we detected a coordinate system that we can not format as a
         * Proj.4 definition string. Format an error message with axis directions and units.
         */
    definition.setLength(0);
    definition.append('(');
    for (int i = 0; i < units.length; i++) {
        final CoordinateSystemAxis axis = cs.getAxis(i);
        if (i != 0)
            definition.append(", ");
        definition.append(axis.getUnit()).append(' ').append(Types.getCodeName(axis.getDirection()));
    }
    throw new FactoryException(Errors.format(Errors.Keys.IllegalCoordinateSystem_1, definition.append(')')));
}
Also used : ParameterValueGroup(org.opengis.parameter.ParameterValueGroup) UnavailableFactoryException(org.apache.sis.referencing.factory.UnavailableFactoryException) FactoryException(org.opengis.util.FactoryException) CoordinateSystem(org.opengis.referencing.cs.CoordinateSystem) Projection(org.opengis.referencing.operation.Projection) CoordinateSystemAxis(org.opengis.referencing.cs.CoordinateSystemAxis) Unit(javax.measure.Unit) Identifier(org.opengis.metadata.Identifier) AxisDirection(org.opengis.referencing.cs.AxisDirection) EllipsoidalCS(org.opengis.referencing.cs.EllipsoidalCS) CartesianCS(org.opengis.referencing.cs.CartesianCS) GeneralParameterValue(org.opengis.parameter.GeneralParameterValue) ParameterValue(org.opengis.parameter.ParameterValue) GeneralParameterValue(org.opengis.parameter.GeneralParameterValue) GeodeticDatum(org.opengis.referencing.datum.GeodeticDatum) DefaultGeodeticDatum(org.apache.sis.referencing.datum.DefaultGeodeticDatum) GeodeticCRS(org.opengis.referencing.crs.GeodeticCRS) PrimeMeridian(org.opengis.referencing.datum.PrimeMeridian) ProjectedCRS(org.opengis.referencing.crs.ProjectedCRS) Angle(javax.measure.quantity.Angle) DefaultGeodeticDatum(org.apache.sis.referencing.datum.DefaultGeodeticDatum) IdentifiedObject(org.opengis.referencing.IdentifiedObject) BursaWolfParameters(org.apache.sis.referencing.datum.BursaWolfParameters) Ellipsoid(org.opengis.referencing.datum.Ellipsoid)

Aggregations

CartesianCS (org.opengis.referencing.cs.CartesianCS)14 EllipsoidalCS (org.opengis.referencing.cs.EllipsoidalCS)7 CoordinateSystem (org.opengis.referencing.cs.CoordinateSystem)6 ProjectedCRS (org.opengis.referencing.crs.ProjectedCRS)5 GeographicCRS (org.opengis.referencing.crs.GeographicCRS)4 Angle (javax.measure.quantity.Angle)3 GeodeticCRS (org.opengis.referencing.crs.GeodeticCRS)3 SphericalCS (org.opengis.referencing.cs.SphericalCS)3 GeodeticDatum (org.opengis.referencing.datum.GeodeticDatum)3 FactoryException (org.opengis.util.FactoryException)3 HashMap (java.util.HashMap)2 Convention (org.apache.sis.io.wkt.Convention)2 BursaWolfParameters (org.apache.sis.referencing.datum.BursaWolfParameters)2 SimpleInternationalString (org.apache.sis.util.iso.SimpleInternationalString)2 ParameterValueGroup (org.opengis.parameter.ParameterValueGroup)2 IdentifiedObject (org.opengis.referencing.IdentifiedObject)2 VerticalCRS (org.opengis.referencing.crs.VerticalCRS)2 CoordinateSystemAxis (org.opengis.referencing.cs.CoordinateSystemAxis)2 PrimeMeridian (org.opengis.referencing.datum.PrimeMeridian)2 InternationalString (org.opengis.util.InternationalString)2