Examples with CoordinateSystem org.opengis.referencing.cs.CoordinateSystem used on opensource projects

Example 11 with CoordinateSystem

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

the class DefaultGeodeticCRS method formatTo.

/**
 * Formats this CRS as a <cite>Well Known Text</cite> {@code GeodeticCRS[…]} element.
 * More information about the WKT format is documented in subclasses.
 *
 * @return {@code "GeodeticCRS"} (WKT 2) or {@code "GeogCS"}/{@code "GeocCS"} (WKT 1).
 */
@Override
protected String formatTo(final Formatter formatter) {
    WKTUtilities.appendName(this, formatter, null);
    CoordinateSystem cs = getCoordinateSystem();
    final Convention convention = formatter.getConvention();
    final boolean isWKT1 = (convention.majorVersion() == 1);
    final boolean isGeographicWKT1 = isWKT1 && (cs instanceof EllipsoidalCS);
    if (isGeographicWKT1 && cs.getDimension() == 3) {
        /*
             * Version 1 of WKT format did not have three-dimensional GeographicCRS. Instead, such CRS were formatted
             * as a CompoundCRS made of a two-dimensional GeographicCRS with a VerticalCRS for the ellipsoidal height.
             * Note that such compound is illegal in WKT 2 and ISO 19111 standard, as ellipsoidal height shall not be
             * separated from the geographic component. So we perform this separation only at WKT 1 formatting time.
             */
        SingleCRS first = CRS.getHorizontalComponent(this);
        SingleCRS second = CRS.getVerticalComponent(this, true);
        if (first != null && second != null) {
            // Should not be null, but we are paranoiac.
            if (AxisDirection.UP.equals(AxisDirections.absolute(cs.getAxis(0).getDirection()))) {
                // It is very unusual to have VerticalCRS first, but our code tries to be robust.
                final SingleCRS t = first;
                first = second;
                second = t;
            }
            formatter.newLine();
            formatter.append(WKTUtilities.toFormattable(first));
            formatter.newLine();
            formatter.append(WKTUtilities.toFormattable(second));
            formatter.newLine();
            return WKTKeywords.Compd_CS;
        }
    }
    /*
         * Unconditionally format the datum element, followed by the prime meridian.
         * The prime meridian is part of datum according ISO 19111, but is formatted
         * as a sibling (rather than a child) element in WKT for historical reasons.
         */
    // Gives subclasses a chance to override.
    final GeodeticDatum datum = getDatum();
    formatter.newLine();
    formatter.append(WKTUtilities.toFormattable(datum));
    formatter.newLine();
    final PrimeMeridian pm = datum.getPrimeMeridian();
    final Unit<Angle> angularUnit = AxisDirections.getAngularUnit(cs, null);
    if (// Really this specific enum, not Convention.isSimplified().
    convention != Convention.WKT2_SIMPLIFIED || ReferencingUtilities.getGreenwichLongitude(pm, Units.DEGREE) != 0) {
        final Unit<Angle> oldUnit = formatter.addContextualUnit(angularUnit);
        formatter.indent(1);
        formatter.append(WKTUtilities.toFormattable(pm));
        formatter.indent(-1);
        formatter.newLine();
        formatter.restoreContextualUnit(angularUnit, oldUnit);
    }
    /*
         * Get the coordinate system to format. This will also determine the units to write and the keyword to
         * return in WKT 1 format. Note that for the WKT 1 format, we need to replace the coordinate system by
         * an instance conform to the legacy conventions.
         *
         * We can not delegate the work below to subclasses,  because XML unmarshalling of a geodetic CRS will
         * NOT create an instance of a subclass (because the distinction between geographic and geocentric CRS
         * is not anymore in ISO 19111:2007).
         */
    final boolean isBaseCRS;
    if (isWKT1) {
        if (!isGeographicWKT1) {
            // If not geographic, then presumed geocentric.
            if (cs instanceof CartesianCS) {
                cs = Legacy.forGeocentricCRS((CartesianCS) cs, true);
            } else {
                // SphericalCS was not supported in WKT 1.
                formatter.setInvalidWKT(cs, null);
            }
        }
        isBaseCRS = false;
    } else {
        isBaseCRS = isBaseCRS(formatter);
    }
    /*
         * Format the coordinate system, except if this CRS is the base CRS of an AbstractDerivedCRS in WKT 2 format.
         * This is because ISO 19162 omits the coordinate system definition of enclosed base CRS in order to simplify
         * the WKT. The 'formatCS(…)' method may write axis unit before or after the axes depending on whether we are
         * formatting WKT version 1 or 2 respectively.
         *
         * Note that even if we do not format the CS, we may still write the units if we are formatting in "simplified"
         * mode (as opposed to the more verbose mode). This looks like the opposite of what we would expect, but this is
         * because formatting the unit here allow us to avoid repeating the unit in projection parameters when this CRS
         * is part of a ProjectedCRS. Note however that in such case, the units to format are the angular units because
         * the linear units will be formatted in the enclosing PROJCS[…] element.
         */
    if (!isBaseCRS || convention == Convention.INTERNAL) {
        // Will also format the axes unit.
        formatCS(formatter, cs, ReferencingUtilities.getUnit(cs), isWKT1);
    } else if (convention.isSimplified()) {
        formatter.append(formatter.toContextualUnit(angularUnit));
    }
    /*
         * For WKT 1, the keyword depends on the subclass: "GeogCS" for GeographicCRS or "GeocCS" for GeocentricCRS.
         * However we can not rely on the subclass for choosing the keyword, because after XML unmarhaling we only
         * have a GeodeticCRS. We need to make the choice in this base class. The CS type is a sufficient criterion.
         */
    if (isWKT1) {
        return isGeographicWKT1 ? WKTKeywords.GeogCS : WKTKeywords.GeocCS;
    } else {
        return isBaseCRS ? WKTKeywords.BaseGeodCRS : formatter.shortOrLong(WKTKeywords.GeodCRS, WKTKeywords.GeodeticCRS);
    }
}

Example 12 with CoordinateSystem

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

the class DefaultCompoundCS method getAxes.

/**
 * Returns all axes in the given sequence of components.
 */
@SuppressWarnings("ForLoopReplaceableByForEach")
private static CoordinateSystemAxis[] getAxes(final CoordinateSystem[] components) {
    int count = 0;
    for (int i = 0; i < components.length; i++) {
        count += components[i].getDimension();
    }
    final CoordinateSystemAxis[] axis = new CoordinateSystemAxis[count];
    count = 0;
    for (final CoordinateSystem c : components) {
        final int dim = c.getDimension();
        for (int j = 0; j < dim; j++) {
            axis[count++] = c.getAxis(j);
        }
    }
    assert count == axis.length;
    return axis;
}

Example 13 with CoordinateSystem

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

the class EPSGFactoryTest method testCreateByName.

/**
 * Tests the creation of CRS using name instead of primary key.
 *
 * @throws FactoryException if an error occurred while querying the factory.
 *
 * @see #testProjectedByName()
 */
@Test
public void testCreateByName() throws FactoryException {
    final EPSGFactory factory = TestFactorySource.factory;
    assumeNotNull(factory);
    assertSame(factory.createUnit("9002"), factory.createUnit("foot"));
    assertNotSame(factory.createUnit("9001"), factory.createUnit("foot"));
    /*
         * Test a name with colons.
         */
    final CoordinateSystem cs = factory.createCoordinateSystem("Ellipsoidal 2D CS. Axes: latitude, longitude. Orientations: north, east. UoM: degree");
    assertEpsgNameAndIdentifierEqual("Ellipsoidal 2D CS. Axes: latitude, longitude. Orientations: north, east. UoM: degree", 6422, cs);
    /*
         * Tests with a unknown name. The exception should be NoSuchAuthorityCodeException
         * (some previous version wrongly threw a SQLException when using HSQL database).
         */
    try {
        factory.createGeographicCRS("WGS83");
        fail("Should not find a geographic CRS named “WGS83” (the actual name is “WGS 84”).");
    } catch (NoSuchAuthorityCodeException e) {
        // This is the expected exception.
        assertEquals("WGS83", e.getAuthorityCode());
    }
}

Example 14 with CoordinateSystem

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

the class AbstractEnvelope method toSimpleEnvelopes.

/**
 * Returns this envelope as an array of simple (without wraparound) envelopes.
 * The length of the returned array depends on the number of dimensions where a
 * {@linkplain org.opengis.referencing.cs.RangeMeaning#WRAPAROUND wraparound} range is found.
 * Typically, wraparound occurs only in the range of longitude values, when the range crosses
 * the anti-meridian (a.k.a. date line). However this implementation will take in account any
 * axis having wraparound {@linkplain CoordinateSystemAxis#getRangeMeaning() range meaning}.
 *
 * <p>Special cases:</p>
 *
 * <ul>
 *   <li>If this envelope {@linkplain #isEmpty() is empty}, then this method returns an empty array.</li>
 *   <li>If this envelope does not have any wraparound behavior, then this method returns {@code this}
 *       in an array of length 1. This envelope is <strong>not</strong> cloned.</li>
 *   <li>If this envelope crosses the <cite>anti-meridian</cite> (a.k.a. <cite>date line</cite>)
 *       then this method represents this envelope as two separated simple envelopes.
 *   <li>While uncommon, the envelope could theoretically crosses the limit of other axis having
 *       wraparound range meaning. If wraparound occur along <var>n</var> axes, then this method
 *       represents this envelope as 2ⁿ separated simple envelopes.
 * </ul>
 *
 * @return a representation of this envelope as an array of non-empty envelope.
 *
 * @see Envelope2D#toRectangles()
 * @see GeneralEnvelope#simplify()
 *
 * @since 0.4
 */
@SuppressWarnings("ReturnOfCollectionOrArrayField")
public Envelope[] toSimpleEnvelopes() {
    // A bitmask of the dimensions having a "wrap around" behavior.
    long isWrapAround = 0;
    CoordinateReferenceSystem crs = null;
    final int dimension = getDimension();
    for (int i = 0; i != dimension; i++) {
        // Do not use getSpan(i).
        final double span = getUpper(i) - getLower(i);
        if (!(span > 0)) {
            // Use '!' for catching NaN.
            if (!isNegative(span)) {
                // Span is positive zero.
                return EMPTY;
            }
            if (crs == null) {
                crs = getCoordinateReferenceSystem();
            }
            if (!isWrapAround(crs, i)) {
                return EMPTY;
            }
            if (i >= Long.SIZE) {
                // a CRS is unusual enough for not being worth to make the distinction in the error message.
                throw new IllegalStateException(Errors.format(Errors.Keys.ExcessiveListSize_2, "axis", dimension));
            }
            isWrapAround |= (1L << i);
        }
    }
    /*
         * The number of simple envelopes is 2ⁿ where n is the number of wraparound found. In most
         * cases, isWrapAround == 0 so we have an array of length 1 containing only this envelope.
         */
    final int bitCount = Long.bitCount(isWrapAround);
    if (bitCount >= Integer.SIZE - 1) {
        // Should be very unusual, but let be paranoiac.
        throw new IllegalStateException(Errors.format(Errors.Keys.ExcessiveListSize_2, "wraparound", bitCount));
    }
    final Envelope[] envelopes = new Envelope[1 << bitCount];
    if (envelopes.length == 1) {
        envelopes[0] = this;
    } else {
        /*
             * Need to create at least 2 envelopes. Instantiate now all envelopes with ordinate values
             * initialized to a copy of this envelope. We will write directly in their internal arrays later.
             */
        double[] c = new double[dimension * 2];
        for (int i = 0; i < dimension; i++) {
            c[i] = getLower(i);
            c[i + dimension] = getUpper(i);
        }
        final double[][] ordinates = new double[envelopes.length][];
        for (int i = 0; i < envelopes.length; i++) {
            final GeneralEnvelope envelope = new GeneralEnvelope(i == 0 ? c : c.clone());
            envelope.crs = crs;
            envelopes[i] = envelope;
            ordinates[i] = envelope.ordinates;
        }
        /*
             * Assign the minimum and maximum ordinate values in the dimension where a wraparound has been found.
             * The 'for' loop below iterates only over the 'i' values for which the 'isWrapAround' bit is set to 1.
             */
        // For identifying whether we need to set the lower or the upper ordinate.
        int mask = 1;
        @SuppressWarnings("null") final CoordinateSystem // Should not be null at this point.
        cs = crs.getCoordinateSystem();
        for (int i; (i = Long.numberOfTrailingZeros(isWrapAround)) != Long.SIZE; isWrapAround &= ~(1L << i)) {
            final CoordinateSystemAxis axis = cs.getAxis(i);
            final double min = axis.getMinimumValue();
            final double max = axis.getMaximumValue();
            for (int j = 0; j < ordinates.length; j++) {
                c = ordinates[j];
                if ((j & mask) == 0) {
                    c[i + dimension] = max;
                } else {
                    c[i] = min;
                }
            }
            mask <<= 1;
        }
    }
    return envelopes;
}

Example 15 with CoordinateSystem

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

the class Shapes2D method transform.

/**
 * Transforms a rectangular envelope using the given coordinate operation.
 * The transformation is only approximative: the returned envelope may be bigger
 * than the smallest possible bounding box, but should not be smaller in most cases.
 *
 * <p>This method can handle the case where the rectangle contains the North or South pole,
 * or when it cross the ±180° longitude.</p>
 *
 * @param  operation    the operation to use. Source and target dimension must be 2.
 * @param  envelope     the rectangle to transform (may be {@code null}).
 * @param  destination  the destination rectangle (may be {@code envelope}).
 *         If {@code null}, a new rectangle will be created and returned.
 * @return {@code destination}, or a new rectangle if {@code destination} was non-null and {@code envelope} was null.
 * @throws TransformException if a transform failed.
 *
 * @see #transform(MathTransform2D, Rectangle2D, Rectangle2D)
 * @see Envelopes#transform(CoordinateOperation, Envelope)
 */
@SuppressWarnings("null")
public static Rectangle2D transform(final CoordinateOperation operation, final Rectangle2D envelope, Rectangle2D destination) throws TransformException {
    ArgumentChecks.ensureNonNull("operation", operation);
    if (envelope == null) {
        return null;
    }
    final MathTransform transform = operation.getMathTransform();
    if (!(transform instanceof MathTransform2D)) {
        throw new MismatchedDimensionException(Errors.format(Errors.Keys.IllegalPropertyValueClass_3, "transform", MathTransform2D.class, MathTransform.class));
    }
    MathTransform2D mt = (MathTransform2D) transform;
    final double[] center = new double[2];
    destination = transform(mt, envelope, destination, center);
    /*
         * If the source envelope crosses the expected range of valid coordinates, also projects
         * the range bounds as a safety. See the comments in transform(Envelope, ...).
         */
    final CoordinateReferenceSystem sourceCRS = operation.getSourceCRS();
    if (sourceCRS != null) {
        final CoordinateSystem cs = sourceCRS.getCoordinateSystem();
        if (cs != null && cs.getDimension() == 2) {
            // Paranoiac check.
            CoordinateSystemAxis axis = cs.getAxis(0);
            double min = envelope.getMinX();
            double max = envelope.getMaxX();
            Point2D.Double pt = null;
            for (int i = 0; i < 4; i++) {
                if (i == 2) {
                    axis = cs.getAxis(1);
                    min = envelope.getMinY();
                    max = envelope.getMaxY();
                }
                final double v = (i & 1) == 0 ? axis.getMinimumValue() : axis.getMaximumValue();
                if (!(v > min && v < max)) {
                    continue;
                }
                if (pt == null) {
                    pt = new Point2D.Double();
                }
                if ((i & 2) == 0) {
                    pt.x = v;
                    pt.y = envelope.getCenterY();
                } else {
                    pt.x = envelope.getCenterX();
                    pt.y = v;
                }
                destination.add(mt.transform(pt, pt));
            }
        }
    }
    /*
         * Now take the target CRS in account.
         */
    final CoordinateReferenceSystem targetCRS = operation.getTargetCRS();
    if (targetCRS == null) {
        return destination;
    }
    final CoordinateSystem targetCS = targetCRS.getCoordinateSystem();
    if (targetCS == null || targetCS.getDimension() != 2) {
        // It should be an error, but we keep this method tolerant.
        return destination;
    }
    /*
         * Checks for singularity points. See the Envelopes.transform(CoordinateOperation, Envelope)
         * method for comments about the algorithm. The code below is the same algorithm adapted for
         * the 2D case and the related objects (Point2D, Rectangle2D, etc.).
         *
         * The 'border' variable in the loop below is used in order to compress 2 dimensions
         * and 2 extremums in a single loop, in this order: (xmin, xmax, ymin, ymax).
         */
    TransformException warning = null;
    Point2D sourcePt = null;
    Point2D targetPt = null;
    // A bitmask for each (dimension, extremum) pairs.
    int includedBoundsValue = 0;
    for (int border = 0; border < 4; border++) {
        // 2 dimensions and 2 extremums compacted in a flag.
        // The dimension index being examined.
        final int dimension = border >>> 1;
        final CoordinateSystemAxis axis = targetCS.getAxis(dimension);
        if (axis == null) {
            // Should never be null, but check as a paranoiac safety.
            continue;
        }
        final double extremum = (border & 1) == 0 ? axis.getMinimumValue() : axis.getMaximumValue();
        if (Double.isInfinite(extremum) || Double.isNaN(extremum)) {
            continue;
        }
        if (targetPt == null) {
            try {
                mt = mt.inverse();
            } catch (NoninvertibleTransformException exception) {
                Envelopes.recoverableException(Shapes2D.class, exception);
                return destination;
            }
            targetPt = new Point2D.Double();
        }
        switch(dimension) {
            case 0:
                targetPt.setLocation(extremum, center[1]);
                break;
            case 1:
                targetPt.setLocation(center[0], extremum);
                break;
            default:
                throw new AssertionError(border);
        }
        try {
            sourcePt = mt.transform(targetPt, sourcePt);
        } catch (TransformException exception) {
            if (warning == null) {
                warning = exception;
            } else {
                warning.addSuppressed(exception);
            }
            continue;
        }
        if (envelope.contains(sourcePt)) {
            destination.add(targetPt);
            includedBoundsValue |= (1 << border);
        }
    }
    /*
         * Iterate over all dimensions of type "WRAPAROUND" for which minimal or maximal axis
         * values have not yet been included in the envelope. We could inline this check inside
         * the above loop, but we don't in order to have a chance to exclude the dimensions for
         * which the point have already been added.
         *
         * See transform(CoordinateOperation, Envelope) for more comments about the algorithm.
         */
    if (includedBoundsValue != 0) {
        /*
             * Bits mask transformation:
             *   1) Swaps the two dimensions               (YyXx  →  XxYy)
             *   2) Insert a space between each bits       (XxYy  →  X.x.Y.y.)
             *   3) Fill the space with duplicated values  (X.x.Y.y.  →  XXxxYYyy)
             *
             * In terms of bit positions 1,2,4,8 (not bit values), we have:
             *
             *   8421  →  22881144
             *   i.e. (ymax, ymin, xmax, xmin)  →  (xmax², ymax², xmin², ymin²)
             *
             * Now look at the last part: (xmin², ymin²). The next step is to perform a bitwise
             * AND operation in order to have only both of the following conditions:
             *
             *   Borders not yet added to the envelope: ~(ymax, ymin, xmax, xmin)
             *   Borders in which a singularity exists:  (xmin, xmin, ymin, ymin)
             *
             * The same operation is repeated on the next 4 bits for (xmax, xmax, ymax, ymax).
             */
        int toTest = ((includedBoundsValue & 1) << 3) | ((includedBoundsValue & 4) >>> 1) | ((includedBoundsValue & 2) << 6) | ((includedBoundsValue & 8) << 2);
        // Duplicate the bit values.
        toTest |= (toTest >>> 1);
        toTest &= ~(includedBoundsValue | (includedBoundsValue << 4));
        /*
             * Forget any axes that are not of kind "WRAPAROUND". Then get the final
             * bit pattern indicating which points to test. Iterate over that bits.
             */
        if ((toTest & 0x33333333) != 0 && !CoordinateOperations.isWrapAround(targetCS.getAxis(0)))
            toTest &= 0xCCCCCCCC;
        if ((toTest & 0xCCCCCCCC) != 0 && !CoordinateOperations.isWrapAround(targetCS.getAxis(1)))
            toTest &= 0x33333333;
        while (toTest != 0) {
            final int border = Integer.numberOfTrailingZeros(toTest);
            final int bitMask = 1 << border;
            // Clear now the bit, for the next iteration.
            toTest &= ~bitMask;
            final int dimensionToAdd = (border >>> 1) & 1;
            final CoordinateSystemAxis toAdd = targetCS.getAxis(dimensionToAdd);
            final CoordinateSystemAxis added = targetCS.getAxis(dimensionToAdd ^ 1);
            double x = (border & 1) == 0 ? toAdd.getMinimumValue() : toAdd.getMaximumValue();
            double y = (border & 4) == 0 ? added.getMinimumValue() : added.getMaximumValue();
            if (dimensionToAdd != 0) {
                final double t = x;
                x = y;
                y = t;
            }
            targetPt.setLocation(x, y);
            try {
                sourcePt = mt.transform(targetPt, sourcePt);
            } catch (TransformException exception) {
                if (warning == null) {
                    warning = exception;
                } else {
                    warning.addSuppressed(exception);
                }
                continue;
            }
            if (envelope.contains(sourcePt)) {
                destination.add(targetPt);
            }
        }
    }
    /*
         * At this point we finished envelope transformation. Verify if some ordinates need to be "wrapped around"
         * as a result of the coordinate operation.   This is usually the longitude axis where the source CRS uses
         * the [-180 … +180]° range and the target CRS uses the [0 … 360]° range, or the converse. In such case we
         * set the rectangle to the full range (we do not use the mechanism documented in Envelope2D) because most
         * Rectangle2D implementations do not support spanning the anti-meridian. This results in larger rectangle
         * than what would be possible with GeneralEnvelope or Envelope2D, but we try to limit the situation where
         * this expansion is applied.
         */
    final Set<Integer> wrapAroundChanges;
    if (operation instanceof AbstractCoordinateOperation) {
        wrapAroundChanges = ((AbstractCoordinateOperation) operation).getWrapAroundChanges();
    } else {
        wrapAroundChanges = CoordinateOperations.wrapAroundChanges(sourceCRS, targetCS);
    }
    for (int dim : wrapAroundChanges) {
        // Empty in the vast majority of cases.
        final CoordinateSystemAxis axis = targetCS.getAxis(dim);
        final double minimum = axis.getMinimumValue();
        final double maximum = axis.getMaximumValue();
        final double o1, o2;
        if (dim == 0) {
            o1 = destination.getMinX();
            o2 = destination.getMaxX();
        } else {
            o1 = destination.getMinY();
            o2 = destination.getMaxY();
        }
        if (o1 < minimum || o2 > maximum) {
            final double span = maximum - minimum;
            if (dim == 0) {
                destination.setRect(minimum, destination.getY(), span, destination.getHeight());
            } else {
                destination.setRect(destination.getX(), minimum, destination.getWidth(), span);
            }
        }
    }
    if (warning != null) {
        Envelopes.recoverableException(Shapes2D.class, warning);
    }
    return destination;
}