Examples with GeodeticCRS org.opengis.referencing.crs.GeodeticCRS used on opensource projects

Example 1 with GeodeticCRS

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

the class CRS method suggestCommonTarget.

/**
 * Suggests a coordinate reference system which could be a common target for coordinate operations having the
 * given sources. This method compares the {@linkplain #getGeographicBoundingBox(CoordinateReferenceSystem)
 * domain of validity} of all given CRSs. If a CRS has a domain of validity that contains the domain of all other
 * CRS, than that CRS is returned. Otherwise this method verifies if a {@linkplain GeneralDerivedCRS#getBaseCRS()
 * base CRS} (usually a {@linkplain org.apache.sis.referencing.crs.DefaultGeographicCRS geographic CRS} instance)
 * would be suitable. If no suitable CRS is found, then this method returns {@code null}.
 *
 * <div class="note"><b>Use case:</b>
 * before to test if two arbitrary envelopes {@linkplain GeneralEnvelope#intersects(Envelope) intersect} each other,
 * they need to be {@linkplain Envelopes#transform(Envelope, CoordinateReferenceSystem) transformed} in the same CRS.
 * However if one CRS is a Transverse Mercator projection while the other CRS is a world-wide geographic CRS, then
 * attempts to use the Transverse Mercator projection as the common CRS is likely to fail since the geographic envelope
 * may span an area far outside the projection domain of validity. This {@code suggestCommonTarget(…)} method can used
 * for choosing a common CRS which is less likely to fail.</div>
 *
 * @param  regionOfInterest  the geographic area for which the coordinate operations will be applied,
 *                           or {@code null} if unknown.
 * @param  sourceCRS         the coordinate reference systems for which a common target CRS is desired.
 * @return a CRS that may be used as a common target for all the given source CRS in the given region of interest,
 *         or {@code null} if this method did not find a common target CRS. The returned CRS may be different than
 *         all given CRS.
 *
 * @since 0.8
 */
public static CoordinateReferenceSystem suggestCommonTarget(GeographicBoundingBox regionOfInterest, CoordinateReferenceSystem... sourceCRS) {
    CoordinateReferenceSystem bestCRS = null;
    /*
         * Compute the union of the domain of validity of all CRS. If a CRS does not specify a domain of validity,
         * then assume that the CRS is valid for the whole world if the CRS is geodetic or return null otherwise.
         * Opportunistically remember the domain of validity of each CRS in this loop since we will need them later.
         */
    boolean worldwide = false;
    DefaultGeographicBoundingBox domain = null;
    final GeographicBoundingBox[] domains = new GeographicBoundingBox[sourceCRS.length];
    for (int i = 0; i < sourceCRS.length; i++) {
        final CoordinateReferenceSystem crs = sourceCRS[i];
        final GeographicBoundingBox bbox = getGeographicBoundingBox(crs);
        if (bbox == null) {
            /*
                 * If no domain of validity is specified and we can not fallback
                 * on some knowledge about what the CRS is, abandon.
                 */
            if (!(crs instanceof GeodeticCRS)) {
                return null;
            }
            /*
                 * Geodetic CRS (geographic or geocentric) can generally be presumed valid in a worldwide area.
                 * The 'worldwide' flag is a little optimization for remembering that we do not need to compute
                 * the union anymore, but we still need to continue the loop for fetching all bounding boxes.
                 */
            // Fallback to be used if we don't find anything better.
            bestCRS = crs;
            worldwide = true;
        } else {
            domains[i] = bbox;
            if (!worldwide) {
                if (domain == null) {
                    domain = new DefaultGeographicBoundingBox(bbox);
                } else {
                    domain.add(bbox);
                }
            }
        }
    }
    /*
         * At this point we got the union of the domain of validity of all CRS. We are interested only in the
         * part that intersect the region of interest. If the union is whole world, we do not need to compute
         * the intersection; we can just leave the region of interest unchanged.
         */
    if (domain != null && !worldwide) {
        if (regionOfInterest != null) {
            domain.intersect(regionOfInterest);
        }
        regionOfInterest = domain;
        domain = null;
    }
    /*
         * Iterate again over the domain of validity of all CRS.  For each domain of validity, compute the area
         * which is inside the domain or interest and the area which is outside. The "best CRS" will be the one
         * which comply with the following rules, in preference order:
         *
         *   1) The CRS which is valid over the largest area of the region of interest.
         *   2) If two CRS are equally good according rule 1, then the CRS with the smallest "outside area".
         *
         * Example: given two source CRS, a geographic one and a projected one:
         *
         *   - If the projected CRS contains fully the region of interest, then it will be returned.
         *     The preference is given to the projected CRS because geometric are likely to be more
         *     accurate in that space. Furthermore forward conversions from geographic to projected
         *     CRS are usually faster than inverse conversions.
         *
         *   - Otherwise (i.e. if the region of interest is likely to be wider than the projected CRS
         *     domain of validity), then the geographic CRS will be returned.
         */
    // NaN if 'regionOfInterest' is null.
    final double roiArea = Extents.area(regionOfInterest);
    double maxInsideArea = 0;
    double minOutsideArea = Double.POSITIVE_INFINITY;
    boolean tryDerivedCRS = false;
    do {
        for (int i = 0; i < domains.length; i++) {
            final GeographicBoundingBox bbox = domains[i];
            if (bbox != null) {
                double insideArea = Extents.area(bbox);
                double outsideArea = 0;
                if (regionOfInterest != null) {
                    if (domain == null) {
                        domain = new DefaultGeographicBoundingBox(bbox);
                    } else {
                        domain.setBounds(bbox);
                    }
                    domain.intersect(regionOfInterest);
                    final double area = insideArea;
                    insideArea = Extents.area(domain);
                    outsideArea = area - insideArea;
                }
                if (insideArea > maxInsideArea || (insideArea == maxInsideArea && outsideArea < minOutsideArea)) {
                    maxInsideArea = insideArea;
                    minOutsideArea = outsideArea;
                    bestCRS = sourceCRS[i];
                }
            }
        }
        /*
             * If the best CRS does not cover fully the region of interest, then we will redo the check again
             * but using base CRS instead. For example if the list of source CRS had some projected CRS, we
             * will try with the geographic CRS on which those projected CRS are based.
             */
        if (maxInsideArea < roiArea) {
            // Do not try twice.
            if (tryDerivedCRS)
                break;
            final CoordinateReferenceSystem[] derivedCRS = new CoordinateReferenceSystem[sourceCRS.length];
            for (int i = 0; i < derivedCRS.length; i++) {
                GeographicBoundingBox bbox = null;
                final CoordinateReferenceSystem crs = sourceCRS[i];
                if (crs instanceof GeneralDerivedCRS) {
                    final CoordinateReferenceSystem baseCRS = ((GeneralDerivedCRS) crs).getBaseCRS();
                    bbox = getGeographicBoundingBox(baseCRS);
                    if (bbox == null && bestCRS == null && baseCRS instanceof GeodeticCRS) {
                        // Fallback to be used if we don't find anything better.
                        bestCRS = baseCRS;
                    }
                    tryDerivedCRS = true;
                    derivedCRS[i] = baseCRS;
                }
                domains[i] = bbox;
            }
            sourceCRS = derivedCRS;
        } else {
            break;
        }
    } while (tryDerivedCRS);
    return bestCRS;
}

Example 2 with GeodeticCRS

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

the class DefaultCompoundCRS method verify.

/**
 * Verifies that the given array does not contain duplicated horizontal or vertical components.
 * Verifies also that if there is an horizontal component, then there is no ellipsoidal height
 * defined separately.
 *
 * @param  properties  the user-specified properties, for determining the locale of error messages.
 * @param  components  the components to verify.
 */
private static void verify(final Map<String, ?> properties, final CoordinateReferenceSystem[] components) {
    int allTypes = 0;
    // 0 for false, 1 for true.
    int isProjected = 0;
    boolean isEllipsoidalHeight = false;
    for (final CoordinateReferenceSystem component : components) {
        final int type;
        if (component instanceof GeodeticCRS) {
            // Must match the number used in Resources.Keys.DuplicatedSpatialComponents_1.
            type = 1;
        } else if (component instanceof ProjectedCRS) {
            isProjected = 1;
            // Intentionally same number than for GeographicCRS case.
            type = 1;
        } else if (component instanceof VerticalCRS) {
            isEllipsoidalHeight = ReferencingUtilities.isEllipsoidalHeight(((VerticalCRS) component).getDatum());
            // Must match the number used in Resources.Keys.DuplicatedSpatialComponents_1.
            type = 2;
        } else {
            // Skip other types. In particular, we allow 2 temporal CRS (used in meteorology).
            continue;
        }
        if (allTypes == (allTypes |= type)) {
            throw new IllegalArgumentException(Resources.forProperties(properties).getString(Resources.Keys.DuplicatedSpatialComponents_1, type));
        }
    }
    if (isEllipsoidalHeight && ((allTypes & 1) != 0)) {
        throw new IllegalArgumentException(Resources.forProperties(properties).getString(Resources.Keys.EllipsoidalHeightNotAllowed_1, isProjected));
    }
}

Example 3 with GeodeticCRS

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

the class CRS method horizontalCode.

/**
 * If the given CRS would quality as horizontal except for its number of dimensions, returns that number.
 * Otherwise returns 0. The number of dimensions can only be 2 or 3.
 */
private static int horizontalCode(final CoordinateReferenceSystem crs) {
    /*
         * In order to determine if the CRS is geographic, checking the CoordinateSystem type is more reliable
         * then checking if the CRS implements the GeographicCRS interface.  This is because the GeographicCRS
         * type did not existed in ISO 19111:2007, so a CRS could be standard-compliant without implementing
         * the GeographicCRS interface.
         */
    boolean isEngineering = false;
    final boolean isGeodetic = (crs instanceof GeodeticCRS);
    if (isGeodetic || crs instanceof ProjectedCRS || (isEngineering = (crs instanceof EngineeringCRS))) {
        final CoordinateSystem cs = crs.getCoordinateSystem();
        final int dim = cs.getDimension();
        if ((dim & ~1) == 2 && (!isGeodetic || (cs instanceof EllipsoidalCS))) {
            if (isEngineering) {
                int n = 0;
                for (int i = 0; i < dim; i++) {
                    if (AxisDirections.isCompass(cs.getAxis(i).getDirection()))
                        n++;
                }
                // If we don't have exactly 2 east, north, etc. directions, consider as non-horizontal.
                if (n != 2)
                    return 0;
            }
            return dim;
        }
    }
    return 0;
}

Example 4 with GeodeticCRS

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

the class SubTypes method castOrCopy.

/**
 * Returns a SIS implementation for the given coordinate reference system.
 *
 * @see AbstractCRS#castOrCopy(CoordinateReferenceSystem)
 */
static AbstractCRS castOrCopy(final CoordinateReferenceSystem object) {
    if (object instanceof DerivedCRS) {
        return DefaultDerivedCRS.castOrCopy((DerivedCRS) object);
    }
    if (object instanceof ProjectedCRS) {
        return DefaultProjectedCRS.castOrCopy((ProjectedCRS) object);
    }
    if (object instanceof GeodeticCRS) {
        if (object instanceof GeographicCRS) {
            return DefaultGeographicCRS.castOrCopy((GeographicCRS) object);
        }
        if (object instanceof GeocentricCRS) {
            return DefaultGeocentricCRS.castOrCopy((GeocentricCRS) object);
        }
        /*
             * The GeographicCRS and GeocentricCRS types are not part of ISO 19111.
             * ISO uses a single type, GeodeticCRS, for both of them and infer the
             * geographic or geocentric type from the coordinate system. We do this
             * check here for instantiating the most appropriate SIS type, but only
             * if we need to create a new object anyway (see below for rational).
             */
        if (object instanceof DefaultGeodeticCRS) {
            /*
                 * Result of XML unmarshalling — keep as-is. We avoid creating a new object because it
                 * would break object identities specified in GML document by the xlink:href attribute.
                 * However we may revisit this policy in the future. See SC_CRS.setElement(AbstractCRS).
                 */
            return (DefaultGeodeticCRS) object;
        }
        final Map<String, ?> properties = IdentifiedObjects.getProperties(object);
        final GeodeticDatum datum = ((GeodeticCRS) object).getDatum();
        final CoordinateSystem cs = object.getCoordinateSystem();
        if (cs instanceof EllipsoidalCS) {
            return new DefaultGeographicCRS(properties, datum, (EllipsoidalCS) cs);
        }
        if (cs instanceof SphericalCS) {
            return new DefaultGeocentricCRS(properties, datum, (SphericalCS) cs);
        }
        if (cs instanceof CartesianCS) {
            return new DefaultGeocentricCRS(properties, datum, (CartesianCS) cs);
        }
    }
    if (object instanceof VerticalCRS) {
        return DefaultVerticalCRS.castOrCopy((VerticalCRS) object);
    }
    if (object instanceof TemporalCRS) {
        return DefaultTemporalCRS.castOrCopy((TemporalCRS) object);
    }
    if (object instanceof EngineeringCRS) {
        return DefaultEngineeringCRS.castOrCopy((EngineeringCRS) object);
    }
    if (object instanceof ImageCRS) {
        return DefaultImageCRS.castOrCopy((ImageCRS) object);
    }
    if (object instanceof CompoundCRS) {
        return DefaultCompoundCRS.castOrCopy((CompoundCRS) object);
    }
    /*
         * Intentionally check for AbstractCRS after the interfaces because user may have defined his own
         * subclass implementing the interface. If we were checking for AbstractCRS before the interfaces,
         * the returned instance could have been a user subclass without the JAXB annotations required
         * for XML marshalling.
         */
    if (object == null || object instanceof AbstractCRS) {
        return (AbstractCRS) object;
    }
    return new AbstractCRS(object);
}

Example 5 with GeodeticCRS

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

the class CoordinateOperationRegistry method propagateVertical.

/**
 * Appends a vertical axis in the source CRS of the first step {@code forward = true} or in
 * the target CRS of the last step {@code forward = false} of the given operations chain.
 *
 * @param  source3D    the potentially three-dimensional source CRS
 * @param  target3D    the potentially three-dimensional target CRS
 * @param  operations  the chain of operations in which to add a vertical axis.
 * @param  forward     {@code true} for adding the vertical axis at the beginning, or
 *                     {@code false} for adding the vertical axis at the end.
 * @return {@code true} on success.
 * @throws IllegalArgumentException if the operation method can not have the desired number of dimensions.
 */
private boolean propagateVertical(final CoordinateReferenceSystem source3D, final CoordinateReferenceSystem target3D, final ListIterator<CoordinateOperation> operations, final boolean forward) throws IllegalArgumentException, FactoryException {
    while (forward ? operations.hasNext() : operations.hasPrevious()) {
        final CoordinateOperation op = forward ? operations.next() : operations.previous();
        /*
             * We will accept to increase the number of dimensions only for operations between geographic CRS.
             * We do not increase the number of dimensions for operations between other kind of CRS because we
             * would not know which value to give to the new dimension.
             */
        CoordinateReferenceSystem sourceCRS, targetCRS;
        if (!((sourceCRS = op.getSourceCRS()) instanceof GeodeticCRS && (targetCRS = op.getTargetCRS()) instanceof GeodeticCRS && sourceCRS.getCoordinateSystem() instanceof EllipsoidalCS && targetCRS.getCoordinateSystem() instanceof EllipsoidalCS)) {
            break;
        }
        /*
             * We can process mostly linear operations, otherwise it is hard to know how to add a dimension.
             * Examples of linear operations are:
             *
             *   - Longitude rotation (EPSG:9601). Note that this is a transformation rather than a conversion.
             *   - Geographic3D to 2D conversion (EPSG:9659).
             *
             * However there is a few special cases where we may be able to add a dimension in a non-linear operation.
             * We can attempt those special cases by just giving the same parameters to the math transform factory
             * together with the desired CRS. Examples of such special cases are:
             *
             *   - Geocentric translations (geog2D domain)
             *   - Coordinate Frame Rotation (geog2D domain)
             *   - Position Vector transformation (geog2D domain)
             */
        Matrix matrix = MathTransforms.getMatrix(op.getMathTransform());
        if (matrix == null) {
            if (SubTypes.isSingleOperation(op)) {
                final MathTransformFactory mtFactory = factorySIS.getMathTransformFactory();
                if (mtFactory instanceof DefaultMathTransformFactory) {
                    if (forward)
                        sourceCRS = toGeodetic3D(sourceCRS, source3D);
                    else
                        targetCRS = toGeodetic3D(targetCRS, target3D);
                    final MathTransform mt;
                    try {
                        mt = ((DefaultMathTransformFactory) mtFactory).createParameterizedTransform(((SingleOperation) op).getParameterValues(), ReferencingUtilities.createTransformContext(sourceCRS, targetCRS, null));
                    } catch (InvalidGeodeticParameterException e) {
                        log(null, e);
                        break;
                    }
                    operations.set(recreate(op, sourceCRS, targetCRS, mt, mtFactory.getLastMethodUsed()));
                    return true;
                }
            }
            break;
        }
        /*
             * We can process only one of the following cases:
             *
             *   - Replace a 2D → 2D operation by a 3D → 3D one (i.e. add a passthrough operation).
             *   - Usually remove (or otherwise edit) the operation that change the number of dimensions
             *     between the 2D and 3D cases.
             */
        final int numRow = matrix.getNumRow();
        final int numCol = matrix.getNumCol();
        // 2D → 2D operation.
        final boolean is2D = (numCol == 3 && numRow == 3);
        if (!(is2D || (// 2D → 3D operation.
        forward ? // 2D → 3D operation.
        (numCol == 3 && numRow == 4) : // 3D → 2D operation.
        (numCol == 4 && numRow == 3)))) {
            break;
        }
        matrix = Matrices.resizeAffine(matrix, 4, 4);
        if (matrix.isIdentity()) {
            operations.remove();
        } else {
            /*
                 * If we can not just remove the operation, build a new one with the expected number of dimensions.
                 * The new operation will inherit the same properties except the identifiers, since it is no longer
                 * conform to the definition provided by the authority.
                 */
            final MathTransform mt = factorySIS.getMathTransformFactory().createAffineTransform(matrix);
            operations.set(recreate(op, toGeodetic3D(sourceCRS, source3D), toGeodetic3D(targetCRS, target3D), mt, null));
        }
        /*
             * If we processed the operation that change the number of dimensions, we are done.
             */
        if (!is2D) {
            return true;
        }
    }
    return false;
}