| /* |
| * Licensed to the Apache Software Foundation (ASF) under one or more |
| * contributor license agreements. See the NOTICE file distributed with |
| * this work for additional information regarding copyright ownership. |
| * The ASF licenses this file to You under the Apache License, Version 2.0 |
| * (the "License"); you may not use this file except in compliance with |
| * the License. You may obtain a copy of the License at |
| * |
| * http://www.apache.org/licenses/LICENSE-2.0 |
| * |
| * Unless required by applicable law or agreed to in writing, software |
| * distributed under the License is distributed on an "AS IS" BASIS, |
| * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. |
| * See the License for the specific language governing permissions and |
| * limitations under the License. |
| */ |
| package org.apache.sis.geometry; |
| |
| /* |
| * Do not add dependency to java.awt.Rectangle2D in this class, because not every platforms |
| * support Java2D (e.g. Android), or applications that do not need it may want to avoid to |
| * force installation of the Java2D module (e.g. JavaFX/SWT). |
| */ |
| import java.util.Set; |
| import java.util.ConcurrentModificationException; |
| import org.opengis.geometry.Envelope; |
| import org.opengis.geometry.DirectPosition; |
| import org.opengis.geometry.MismatchedDimensionException; |
| import org.opengis.referencing.cs.CoordinateSystem; |
| import org.opengis.referencing.cs.CoordinateSystemAxis; |
| import org.opengis.referencing.crs.CoordinateReferenceSystem; |
| import org.opengis.referencing.operation.CoordinateOperation; |
| import org.opengis.referencing.operation.Matrix; |
| import org.opengis.referencing.operation.MathTransform; |
| import org.opengis.referencing.operation.TransformException; |
| import org.opengis.referencing.operation.OperationNotFoundException; |
| import org.opengis.referencing.operation.NoninvertibleTransformException; |
| import org.opengis.util.FactoryException; |
| import org.opengis.metadata.extent.GeographicBoundingBox; |
| import org.apache.sis.metadata.iso.extent.DefaultGeographicBoundingBox; |
| import org.apache.sis.referencing.CRS; |
| import org.apache.sis.referencing.operation.AbstractCoordinateOperation; |
| import org.apache.sis.referencing.operation.transform.AbstractMathTransform; |
| import org.apache.sis.internal.metadata.ReferencingServices; |
| import org.apache.sis.internal.referencing.CoordinateOperations; |
| import org.apache.sis.internal.referencing.DirectPositionView; |
| import org.apache.sis.internal.system.Loggers; |
| import org.apache.sis.util.logging.Logging; |
| import org.apache.sis.util.resources.Errors; |
| import org.apache.sis.util.ArgumentChecks; |
| import org.apache.sis.util.ComparisonMode; |
| import org.apache.sis.util.Utilities; |
| import org.apache.sis.util.Static; |
| import org.apache.sis.math.MathFunctions; |
| |
| import static org.apache.sis.util.StringBuilders.trimFractionalPart; |
| |
| |
| /** |
| * Transforms envelopes to new Coordinate Reference Systems, and miscellaneous utilities. |
| * |
| * <div class="section">Envelope transformations</div> |
| * All {@code transform(…)} methods in this class take in account the curvature of the transformed shape. |
| * For example the shape of a geographic envelope (figure below on the left side) is not rectangular in a |
| * conic projection (figure below on the right side). In order to get the envelope represented by the red |
| * rectangle, projecting the four corners of the geographic envelope is not sufficient since we would miss |
| * the southerner part. |
| * |
| * <table class="sis"> |
| * <caption>Example of curvature induced by a map projection</caption> |
| * <tr> |
| * <th>Envelope before map projection</th> |
| * <th>Shape of the projected envelope</th> |
| * </tr><tr> |
| * <td><img src="doc-files/GeographicArea.png" alt="Envelope in a geographic CRS"></td> |
| * <td><img src="doc-files/ConicArea.png" alt="Shape of the envelope transformed in a conic projection"></td> |
| * </tr> |
| * </table> |
| * |
| * Apache SIS tries to detect the curvature by transforming intermediate points in addition to the corners. |
| * While optional, it is strongly recommended that all {@code MathTransform} implementations involved in the |
| * operation (directly or indirectly) support {@linkplain MathTransform#derivative(DirectPosition) derivative}, |
| * for more accurate calculation of curve extremum. This is the case of most Apache SIS implementations. |
| * |
| * <p>The {@code transform(…)} methods in this class expect an arbitrary {@link Envelope} with <strong>one</strong> |
| * of the following arguments: {@link MathTransform}, {@link CoordinateOperation} or {@link CoordinateReferenceSystem}. |
| * The recommended method is the one expecting a {@code CoordinateOperation} object, |
| * since it contains sufficient information for handling the cases of envelopes that encompass a pole. |
| * The method expecting a {@code CoordinateReferenceSystem} object is merely a convenience method that |
| * infers the coordinate operation itself, but at the cost of performance if the same operation needs |
| * to be applied on many envelopes.</p> |
| * |
| * @author Martin Desruisseaux (IRD, Geomatys) |
| * @author Johann Sorel (Geomatys) |
| * @version 1.0 |
| * |
| * @see org.apache.sis.metadata.iso.extent.Extents |
| * @see CRS |
| * |
| * @since 0.3 |
| * @module |
| */ |
| public final class Envelopes extends Static { |
| /** |
| * Do not allow instantiation of this class. |
| */ |
| private Envelopes() { |
| } |
| |
| /** |
| * Puts together a list of envelopes, each of them using an independent coordinate reference system. |
| * The dimension of the returned envelope is the sum of the dimension of all components. |
| * If all components have a coordinate reference system, then the returned envelope will |
| * have a compound coordinate reference system. |
| * |
| * @param components the envelopes to aggregate in a single envelope, in the given order. |
| * @return the aggregation of all given envelopes. |
| * @throws FactoryException if the geodetic factory failed to create the compound CRS. |
| * |
| * @see org.apache.sis.referencing.CRS#compound(CoordinateReferenceSystem...) |
| * @see org.apache.sis.referencing.operation.transform.MathTransforms#compound(MathTransform...) |
| * |
| * @since 1.0 |
| */ |
| public static Envelope compound(final Envelope... components) throws FactoryException { |
| ArgumentChecks.ensureNonNull("components", components); |
| int sum = 0; |
| for (int i=0; i<components.length; i++) { |
| final Envelope env = components[i]; |
| ArgumentChecks.ensureNonNullElement("components", i, env); |
| sum += env.getDimension(); |
| } |
| final GeneralEnvelope compound = new GeneralEnvelope(sum); |
| CoordinateReferenceSystem[] crsComponents = null; |
| int firstAffectedCoordinate = 0; |
| for (int i=0; i<components.length; i++) { |
| final Envelope env = components[i]; |
| final int dim = env.getDimension(); |
| compound.subEnvelope(firstAffectedCoordinate, firstAffectedCoordinate += dim).setEnvelope(env); |
| if (i == 0) { |
| final CoordinateReferenceSystem crs = env.getCoordinateReferenceSystem(); |
| if (crs != null) { |
| crsComponents = new CoordinateReferenceSystem[components.length]; |
| crsComponents[0] = crs; |
| } |
| } else if (crsComponents != null) { |
| if ((crsComponents[i] = env.getCoordinateReferenceSystem()) == null) { |
| crsComponents = null; // Not all CRS are non-null. |
| } |
| } |
| } |
| if (firstAffectedCoordinate != sum) { |
| // Should never happen unless the number of dimensions of an envelope changed during iteration. |
| throw new ConcurrentModificationException(); |
| } |
| if (crsComponents != null) { |
| compound.setCoordinateReferenceSystem(CRS.compound(crsComponents)); |
| } |
| return compound; |
| } |
| |
| /** |
| * Computes the union of all given envelopes, transforming them to a common CRS if necessary. |
| * If all envelopes use the same CRS ({@link ComparisonMode#IGNORE_METADATA ignoring metadata}) |
| * or if the CRS of all envelopes is {@code null}, then the {@linkplain GeneralEnvelope#add(Envelope) |
| * union is computed} without transforming any envelope. Otherwise all envelopes are transformed to a |
| * {@linkplain CRS#suggestCommonTarget common CRS} before union is computed. |
| * The CRS of the returned envelope may different than the CRS of all given envelopes. |
| * |
| * @param envelopes the envelopes for which to compute union. Null elements are ignored. |
| * @return union of given envelopes, or {@code null} if the given array does not contain non-null elements. |
| * @throws TransformException if this method can not determine a common CRS, or if a transformation failed. |
| * |
| * @see GeneralEnvelope#add(Envelope) |
| * |
| * @since 1.0 |
| */ |
| public static GeneralEnvelope union(final Envelope... envelopes) throws TransformException { |
| return EnvelopeReducer.UNION.reduce(envelopes); |
| } |
| |
| /** |
| * Finds a mathematical operation from the CRS of the given source envelope to the CRS of the given target envelope. |
| * For non-null georeferenced envelopes, this method is equivalent to the following code with {@code areaOfInterest} |
| * computed as the union of the two envelopes: |
| * |
| * <blockquote><code>return {@linkplain CRS#findOperation(CoordinateReferenceSystem, CoordinateReferenceSystem, |
| * GeographicBoundingBox)} CRS.findOperation(source.getCoordinateReferenceSystem(), target.getCoordinateReferenceSystem(), |
| * <var>areaOfInterest</var>)</code></blockquote> |
| * |
| * If at least one envelope is null or has no CRS, then this method returns {@code null}. |
| * |
| * @param source the source envelope, or {@code null}. |
| * @param target the target envelope, or {@code null}. |
| * @return the mathematical operation from {@code source} CRS to {@code target} CRS, |
| * or {@code null} if at least one argument is null or has no CRS. |
| * @throws OperationNotFoundException if no operation was found between the given pair of CRS. |
| * @throws FactoryException if the operation can not be created for another reason. |
| * |
| * @see CRS#findOperation(CoordinateReferenceSystem, CoordinateReferenceSystem, GeographicBoundingBox) |
| * |
| * @since 1.0 |
| */ |
| public static CoordinateOperation findOperation(final Envelope source, final Envelope target) throws FactoryException { |
| if (source != null && target != null) { |
| final CoordinateReferenceSystem sourceCRS, targetCRS; |
| if ((sourceCRS = source.getCoordinateReferenceSystem()) != null && |
| (targetCRS = target.getCoordinateReferenceSystem()) != null) |
| { |
| /* |
| * Computing the area of interest (AOI) unconditionally would be harmless, but it is useless if the CRS |
| * are the same since the result will be identity anyway. Conversely we could skip AOI computation more |
| * often for example with the following condition instead of !=: |
| * |
| * !Utilities.deepEquals(sourceCRS, targetCRS, ComparisonMode.ALLOW_VARIANT) |
| * |
| * but it would not be very advantageous if testing the condition is almost as long as computing the AOI. |
| * For now we keep != as a very cheap test which will work quite often. |
| */ |
| DefaultGeographicBoundingBox areaOfInterest = null; |
| if (sourceCRS != targetCRS) try { |
| final DefaultGeographicBoundingBox targetAOI; |
| final ReferencingServices converter = ReferencingServices.getInstance(); |
| areaOfInterest = converter.setBounds(source, null, "findOperation"); // Should be first. |
| targetAOI = converter.setBounds(target, null, "findOperation"); |
| if (areaOfInterest == null) { |
| areaOfInterest = targetAOI; |
| } else if (targetAOI != null) { |
| areaOfInterest.add(targetAOI); |
| } |
| } catch (TransformException e) { |
| /* |
| * Note: we may succeed to transform 'source' and fail to transform 'target' to geographic bounding box, |
| * but the opposite is unlikely because 'source' should not have less dimensions than 'target'. |
| */ |
| Logging.recoverableException(Logging.getLogger(Loggers.GEOMETRY), Envelopes.class, "findOperation", e); |
| } |
| return CRS.findOperation(sourceCRS, targetCRS, areaOfInterest); |
| } |
| } |
| return null; |
| } |
| |
| /** |
| * Computes the intersection of all given envelopes, transforming them to a common CRS if necessary. |
| * If all envelopes use the same CRS ({@link ComparisonMode#IGNORE_METADATA ignoring metadata}) |
| * or if the CRS of all envelopes is {@code null}, then the {@linkplain GeneralEnvelope#intersect(Envelope) |
| * intersection is computed} without transforming any envelope. Otherwise all envelopes are transformed to a |
| * {@linkplain CRS#suggestCommonTarget common CRS} before intersection is computed. |
| * The CRS of the returned envelope may different than the CRS of all given envelopes. |
| * |
| * @param envelopes the envelopes for which to compute intersection. Null elements are ignored. |
| * @return intersection of given envelopes, or {@code null} if the given array does not contain non-null elements. |
| * @throws TransformException if this method can not determine a common CRS, or if a transformation failed. |
| * |
| * @see GeneralEnvelope#intersect(Envelope) |
| * |
| * @since 1.0 |
| */ |
| public static GeneralEnvelope intersect(final Envelope... envelopes) throws TransformException { |
| return EnvelopeReducer.INTERSECT.reduce(envelopes); |
| } |
| |
| /** |
| * Invoked when a recoverable exception occurred. Those exceptions must be minor enough |
| * that they can be silently ignored in most cases. |
| */ |
| static void recoverableException(final Class<? extends Static> caller, final TransformException exception) { |
| Logging.recoverableException(Logging.getLogger(Loggers.GEOMETRY), caller, "transform", exception); |
| } |
| |
| /** |
| * A buckle method for calculating derivative and coordinate transformation in a single step, |
| * if the given {@code derivative} argument is {@code true}. |
| * |
| * @param transform the transform to use. |
| * @param srcPts the array containing the source coordinate at offset 0. |
| * @param dstPts the array into which the transformed coordinate is returned. |
| * @param dstOff the offset to the location of the transformed point that is stored in the destination array. |
| * @param derivate {@code true} for computing the derivative, or {@code false} if not needed. |
| * @return the matrix of the transform derivative at the given source position, |
| * or {@code null} if the {@code derivate} argument is {@code false}. |
| * @throws TransformException if the point can not be transformed |
| * or if a problem occurred while calculating the derivative. |
| */ |
| @SuppressWarnings("null") |
| static Matrix derivativeAndTransform(final MathTransform transform, final double[] srcPts, |
| final double[] dstPts, final int dstOff, final boolean derivate) throws TransformException |
| { |
| if (transform instanceof AbstractMathTransform) { |
| return ((AbstractMathTransform) transform).transform(srcPts, 0, dstPts, dstOff, derivate); |
| } |
| // Derivative must be calculated before to transform the coordinate. |
| final Matrix derivative = derivate ? transform.derivative( |
| new DirectPositionView.Double(srcPts, 0, transform.getSourceDimensions())) : null; |
| transform.transform(srcPts, 0, dstPts, dstOff, 1); |
| return derivative; |
| } |
| |
| /** |
| * Transforms the given envelope to the specified CRS. If any argument is null, or if the |
| * {@linkplain GeneralEnvelope#getCoordinateReferenceSystem() envelope CRS} is null or the |
| * same instance than the given target CRS, then the given envelope is returned unchanged. |
| * Otherwise a new transformed envelope is returned. |
| * |
| * <div class="section">Performance tip</div> |
| * If there is many envelopes to transform with the same source and target CRS, then it is more efficient |
| * to get the {@link CoordinateOperation} or {@link MathTransform} instance once and invoke one of the |
| * others {@code transform(…)} methods. |
| * |
| * @param envelope the envelope to transform (may be {@code null}). |
| * @param targetCRS the target CRS (may be {@code null}). |
| * @return a new transformed envelope, or directly {@code envelope} if no change was required. |
| * @throws TransformException if a transformation was required and failed. |
| * |
| * @since 0.5 |
| */ |
| public static Envelope transform(Envelope envelope, final CoordinateReferenceSystem targetCRS) |
| throws TransformException |
| { |
| if (envelope != null && targetCRS != null) { |
| final CoordinateReferenceSystem sourceCRS = envelope.getCoordinateReferenceSystem(); |
| if (sourceCRS != targetCRS) { |
| if (sourceCRS == null) { |
| // Slight optimization: just copy the given Envelope. |
| envelope = new GeneralEnvelope(envelope); |
| ((GeneralEnvelope) envelope).setCoordinateReferenceSystem(targetCRS); |
| } else { |
| // TODO: create an CoordinateOperationContext with the envelope as geographic area. |
| // May require that we optimize the search for CoordinateOperation with non-null context first. |
| // See https://issues.apache.org/jira/browse/SIS-442 |
| final CoordinateOperation operation; |
| try { |
| operation = CoordinateOperations.factory().createOperation(sourceCRS, targetCRS); |
| } catch (FactoryException exception) { |
| throw new TransformException(Errors.format(Errors.Keys.CanNotTransformEnvelope), exception); |
| } |
| envelope = transform(operation, envelope); |
| } |
| assert Utilities.deepEquals(targetCRS, envelope.getCoordinateReferenceSystem(), ComparisonMode.DEBUG); |
| } |
| } |
| return envelope; |
| } |
| |
| /** |
| * Transforms an envelope using the given math transform. |
| * The transformation is only approximated: the returned envelope may be bigger than necessary, |
| * or smaller than required if the bounding box contains a pole. |
| * The coordinate reference system of the returned envelope will be null. |
| * |
| * <div class="section">Limitation</div> |
| * This method can not handle the case where the envelope contains the North or South pole, |
| * or when it crosses the ±180° longitude, because {@link MathTransform} does not carry sufficient information. |
| * For a more robust envelope transformation, use {@link #transform(CoordinateOperation, Envelope)} instead. |
| * |
| * @param transform the transform to use. |
| * @param envelope envelope to transform, or {@code null}. This envelope will not be modified. |
| * @return the transformed envelope, or {@code null} if {@code envelope} was null. |
| * @throws TransformException if a transform failed. |
| * |
| * @see #transform(CoordinateOperation, Envelope) |
| * |
| * @since 0.5 |
| */ |
| public static GeneralEnvelope transform(final MathTransform transform, final Envelope envelope) |
| throws TransformException |
| { |
| ArgumentChecks.ensureNonNull("transform", transform); |
| return (envelope != null) ? transform(transform, envelope, null) : null; |
| } |
| |
| /** |
| * Implementation of {@link #transform(MathTransform, Envelope)} with the opportunity to |
| * save the projected center coordinate. |
| * |
| * @param targetPt after this method call, the center of the source envelope projected to the target CRS. |
| * The length of this array must be the number of target dimensions. |
| * May be {@code null} if this information is not needed. |
| */ |
| @SuppressWarnings("null") |
| private static GeneralEnvelope transform(final MathTransform transform, |
| final Envelope envelope, |
| final double[] targetPt) |
| throws TransformException |
| { |
| if (transform.isIdentity()) { |
| /* |
| * Slight optimization: Just copy the envelope. Note that we need to set the CRS |
| * to null because we don't know what the target CRS was supposed to be. Even if |
| * an identity transform often imply that the target CRS is the same one than the |
| * source CRS, it is not always the case. The metadata may be differents, or the |
| * transform may be a datum shift without Bursa-Wolf parameters, etc. |
| */ |
| final GeneralEnvelope transformed = new GeneralEnvelope(envelope); |
| transformed.setCoordinateReferenceSystem(null); |
| if (targetPt != null) { |
| for (int i=envelope.getDimension(); --i>=0;) { |
| targetPt[i] = transformed.getMedian(i); |
| } |
| } |
| return transformed; |
| } |
| /* |
| * Checks argument validity: envelope and math transform dimensions must be consistent. |
| */ |
| final int sourceDim = transform.getSourceDimensions(); |
| final int targetDim = transform.getTargetDimensions(); |
| if (envelope.getDimension() != sourceDim) { |
| throw new MismatchedDimensionException(Errors.format(Errors.Keys.MismatchedDimension_2, |
| sourceDim, envelope.getDimension())); |
| } |
| /* |
| * Allocates all needed objects. The power of 3 below is because the following `while` loop |
| * uses a `pointIndex` to be interpreted as a number in base 3 (see the comment inside the loop). |
| * The coordinate tuple to transform must be initialized to the minimal coordinate values. |
| * This coordinate will be updated in the `switch` statement inside the `while` loop. |
| */ |
| if (sourceDim >= 20) { // Maximal value supported by Formulas.pow3(int) is 19. |
| throw new IllegalArgumentException(Errors.format(Errors.Keys.ExcessiveNumberOfDimensions_1)); |
| } |
| int pointIndex = 0; |
| boolean isDerivativeSupported = true; |
| GeneralEnvelope transformed = null; |
| final Matrix[] derivatives = new Matrix[Math.toIntExact(MathFunctions.pow(3, sourceDim))]; |
| final double[] coordinates = new double[derivatives.length * targetDim]; |
| final double[] sourcePt = new double[sourceDim]; |
| for (int i=sourceDim; --i>=0;) { |
| sourcePt[i] = envelope.getMinimum(i); |
| } |
| // A window over a single coordinate in the `coordinates` array. |
| final DirectPositionView ordinatesView = new DirectPositionView.Double(coordinates, 0, targetDim); |
| /* |
| * Iterates over every minimal, maximal and median coordinate values (3 points) along each |
| * dimension. The total number of iterations is 3 ^ (number of source dimensions). |
| */ |
| transformPoint: while (true) { |
| /* |
| * Compute the derivative (optional operation). If this operation fails, we will |
| * set a flag to 'false' so we don't try again for all remaining points. We try |
| * to compute the derivative and the transformed point in a single operation if |
| * we can. If we can not, we will compute those two information separately. |
| * |
| * Note that the very last point to be projected must be the envelope center. |
| * There is usually no need to calculate the derivative for that last point, |
| * but we let it does anyway for safety. |
| */ |
| final int offset = pointIndex * targetDim; |
| try { |
| derivatives[pointIndex] = derivativeAndTransform(transform, |
| sourcePt, coordinates, offset, isDerivativeSupported); |
| } catch (TransformException e) { |
| if (!isDerivativeSupported) { |
| throw e; // Derivative were already disabled, so something went wrong. |
| } |
| isDerivativeSupported = false; |
| transform.transform(sourcePt, 0, coordinates, offset, 1); |
| recoverableException(Envelopes.class, e); // Log only if the above call was successful. |
| } |
| /* |
| * The transformed point has been saved for future reuse after the enclosing |
| * 'while' loop. Now add the transformed point to the destination envelope. |
| */ |
| if (transformed == null) { |
| transformed = new GeneralEnvelope(targetDim); |
| for (int i=0; i<targetDim; i++) { |
| final double value = coordinates[offset + i]; |
| transformed.setRange(i, value, value); |
| } |
| } else { |
| ordinatesView.offset = offset; |
| transformed.add(ordinatesView); |
| } |
| /* |
| * Get the next point coordinate. The 'coordinateIndex' variable is an index in base 3 |
| * having a number of digits equals to the number of source dimensions. For example a |
| * 4-D space have indexes ranging from "0000" to "2222" (numbers in base 3). The digits |
| * are then mapped to minimal (0), maximal (1) or central (2) coordinates. The outer loop |
| * stops when the counter roll back to "0000". Note that 'targetPt' must keep the value |
| * of the last projected point, which must be the envelope center identified by "2222" |
| * in the 4-D case. |
| */ |
| int indexBase3 = ++pointIndex; |
| for (int dim=sourceDim; --dim>=0; indexBase3 /= 3) { |
| switch (indexBase3 % 3) { |
| case 0: sourcePt[dim] = envelope.getMinimum(dim); break; // Continue the loop. |
| case 1: sourcePt[dim] = envelope.getMaximum(dim); continue transformPoint; |
| case 2: sourcePt[dim] = envelope.getMedian (dim); continue transformPoint; |
| default: throw new AssertionError(indexBase3); // Should never happen |
| } |
| } |
| break; |
| } |
| assert pointIndex == derivatives.length : pointIndex; |
| /* |
| * At this point we finished to build an envelope from all sampled positions. Now iterate |
| * over all points. For each point, iterate over all line segments from that point to a |
| * neighbor median point. Use the derivate information for approximating the transform |
| * behavior in that area by a cubic curve. We can then find analytically the curve extremum. |
| * |
| * The same technic is applied in transform(MathTransform, Rectangle2D), except that in |
| * the Rectangle2D case the calculation was bundled right inside the main loop in order |
| * to avoid the need for storage. |
| */ |
| DirectPosition temporary = null; |
| final DirectPositionView sourceView = new DirectPositionView.Double(sourcePt, 0, sourceDim); |
| final CurveExtremum extremum = new CurveExtremum(); |
| for (pointIndex=0; pointIndex < derivatives.length; pointIndex++) { |
| final Matrix D1 = derivatives[pointIndex]; |
| if (D1 != null) { |
| int indexBase3 = pointIndex, power3 = 1; |
| for (int i=sourceDim; --i>=0; indexBase3 /= 3, power3 *= 3) { |
| final int digitBase3 = indexBase3 % 3; |
| if (digitBase3 != 2) { // Process only if we are not already located on the median along the dimension i. |
| final int medianIndex = pointIndex + power3 * (2 - digitBase3); |
| final Matrix D2 = derivatives[medianIndex]; |
| if (D2 != null) { |
| final double xmin = envelope.getMinimum(i); |
| final double xmax = envelope.getMaximum(i); |
| final double x2 = envelope.getMedian (i); |
| final double x1 = (digitBase3 == 0) ? xmin : xmax; |
| final int offset1 = targetDim * pointIndex; |
| final int offset2 = targetDim * medianIndex; |
| for (int j=0; j<targetDim; j++) { |
| extremum.resolve(x1, coordinates[offset1 + j], D1.getElement(j,i), |
| x2, coordinates[offset2 + j], D2.getElement(j,i)); |
| boolean isP2 = false; |
| do { // Executed exactly twice, one for each extremum point. |
| final double x = isP2 ? extremum.ex2 : extremum.ex1; |
| if (x > xmin && x < xmax) { |
| final double y = isP2 ? extremum.ey2 : extremum.ey1; |
| if (y < transformed.getMinimum(j) || |
| y > transformed.getMaximum(j)) |
| { |
| /* |
| * At this point, we have determined that adding the extremum point |
| * would expand the envelope. However we will not add that point |
| * directly because its position may not be quite right (since we |
| * used a cubic curve approximation). Instead, we project the point |
| * on the envelope border which is located vis-à-vis the extremum. |
| */ |
| for (int ib3 = pointIndex, dim = sourceDim; --dim >= 0; ib3 /= 3) { |
| final double coordinate; |
| if (dim == i) { |
| coordinate = x; // Position of the extremum. |
| } else switch (ib3 % 3) { |
| case 0: coordinate = envelope.getMinimum(dim); break; |
| case 1: coordinate = envelope.getMaximum(dim); break; |
| case 2: coordinate = envelope.getMedian (dim); break; |
| default: throw new AssertionError(ib3); // Should never happen. |
| } |
| sourcePt[dim] = coordinate; |
| } |
| temporary = transform.transform(sourceView, temporary); |
| transformed.add(temporary); |
| } |
| } |
| } while ((isP2 = !isP2) == true); |
| } |
| } |
| } |
| } |
| derivatives[pointIndex] = null; // Let GC do its job earlier. |
| } |
| } |
| if (targetPt != null) { |
| // Copy the coordinate of the center point. |
| System.arraycopy(coordinates, coordinates.length - targetDim, targetPt, 0, targetDim); |
| } |
| return transformed; |
| } |
| |
| /** |
| * Transforms an envelope using the given coordinate operation. |
| * The transformation is only approximated: 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 envelope contains the North or South pole, |
| * or when it cross the ±180° longitude.</p> |
| * |
| * <div class="note"><b>Note:</b> |
| * If the envelope CRS is non-null, then the caller should ensure that the operation source CRS |
| * is the same than the envelope CRS. In case of mismatch, this method transforms the envelope |
| * to the operation source CRS before to apply the operation. This extra step may cause a lost |
| * of accuracy. In order to prevent this method from performing such pre-transformation (if not desired), |
| * callers can ensure that the envelope CRS is {@code null} before to call this method.</div> |
| * |
| * @param operation the operation to use. |
| * @param envelope envelope to transform, or {@code null}. This envelope will not be modified. |
| * @return the transformed envelope, or {@code null} if {@code envelope} was null. |
| * @throws TransformException if a transform failed. |
| * |
| * @see #transform(MathTransform, Envelope) |
| * |
| * @since 0.5 |
| */ |
| @SuppressWarnings("null") |
| public static GeneralEnvelope transform(final CoordinateOperation operation, Envelope envelope) |
| throws TransformException |
| { |
| ArgumentChecks.ensureNonNull("operation", operation); |
| if (envelope == null) { |
| return null; |
| } |
| boolean isOperationComplete = true; |
| final CoordinateReferenceSystem sourceCRS = operation.getSourceCRS(); |
| if (sourceCRS != null) { |
| final CoordinateReferenceSystem crs = envelope.getCoordinateReferenceSystem(); |
| if (crs != null && !Utilities.equalsIgnoreMetadata(crs, sourceCRS)) { |
| /* |
| * Argument-check: the envelope CRS seems inconsistent with the given operation. |
| * However we need to push the check a little bit further, since 3D-GeographicCRS |
| * are considered not equal to CompoundCRS[2D-GeographicCRS + ellipsoidal height]. |
| * Checking for identity MathTransform is a more powerfull (but more costly) check. |
| * Since we have the MathTransform, perform an opportunist envelope transform if it |
| * happen to be required. |
| */ |
| final MathTransform mt; |
| try { |
| mt = CoordinateOperations.factory().createOperation(crs, sourceCRS).getMathTransform(); |
| } catch (FactoryException e) { |
| throw new TransformException(Errors.format(Errors.Keys.CanNotTransformEnvelope), e); |
| } |
| if (!mt.isIdentity()) { |
| isOperationComplete = false; |
| envelope = transform(mt, envelope); |
| } |
| } |
| } |
| MathTransform mt = operation.getMathTransform(); |
| final double[] centerPt = new double[mt.getTargetDimensions()]; |
| final GeneralEnvelope transformed = transform(mt, envelope, centerPt); |
| /* |
| * If the source envelope crosses the expected range of valid coordinates, also projects |
| * the range bounds as a safety. Example: if the source envelope goes from 150 to 200°E, |
| * some map projections will interpret 200° as if it was -160°, and consequently produce |
| * an envelope which do not include the 180°W extremum. We will add those extremum points |
| * explicitly as a safety. It may leads to bigger than necessary target envelope, but the |
| * contract is to include at least the source envelope, not to return the smallest one. |
| */ |
| if (sourceCRS != null) { |
| final CoordinateSystem cs = sourceCRS.getCoordinateSystem(); |
| if (cs != null) { // Should never be null, but check as a paranoiac safety. |
| DirectPosition sourcePt = null; |
| DirectPosition targetPt = null; |
| final int dimension = cs.getDimension(); |
| for (int i=0; i<dimension; i++) { |
| final CoordinateSystemAxis axis = cs.getAxis(i); |
| if (axis == null) { // Should never be null, but check as a paranoiac safety. |
| continue; |
| } |
| final double min = envelope.getMinimum(i); |
| final double max = envelope.getMaximum(i); |
| final double v1 = axis.getMinimumValue(); |
| final double v2 = axis.getMaximumValue(); |
| final boolean b1 = (v1 > min && v1 < max); |
| final boolean b2 = (v2 > min && v2 < max); |
| if (!b1 && !b2) { |
| continue; |
| } |
| if (sourcePt == null) { |
| sourcePt = new GeneralDirectPosition(dimension); |
| for (int j=0; j<dimension; j++) { |
| sourcePt.setOrdinate(j, envelope.getMedian(j)); |
| } |
| } |
| if (b1) { |
| sourcePt.setOrdinate(i, v1); |
| transformed.add(targetPt = mt.transform(sourcePt, targetPt)); |
| } |
| if (b2) { |
| sourcePt.setOrdinate(i, v2); |
| transformed.add(targetPt = mt.transform(sourcePt, targetPt)); |
| } |
| sourcePt.setOrdinate(i, envelope.getMedian(i)); |
| } |
| } |
| } |
| /* |
| * Now takes the target CRS in account... |
| */ |
| final CoordinateReferenceSystem targetCRS = operation.getTargetCRS(); |
| if (targetCRS == null) { |
| return transformed; |
| } |
| transformed.setCoordinateReferenceSystem(targetCRS); |
| final CoordinateSystem targetCS = targetCRS.getCoordinateSystem(); |
| if (targetCS == null) { |
| // It should be an error, but we keep this method tolerant. |
| return transformed; |
| } |
| /* |
| * Checks for singularity points. For example the south pole is a singularity point in |
| * geographic CRS because is is located at the maximal value allowed by one particular |
| * axis, namely latitude. This point is not a singularity in the stereographic projection, |
| * because axes extends toward infinity in all directions (mathematically) and because the |
| * South pole has nothing special apart being the origin (0,0). |
| * |
| * Algorithm: |
| * |
| * 1) Inspect the target axis, looking if there is any bounds. If bounds are found, get |
| * the coordinates of singularity points and project them from target to source CRS. |
| * |
| * Example: If the transformed envelope above is (80 … 85°S, 10 … 50°W), and if the |
| * latitude in the target CRS is bounded at 90°S, then project (90°S, 30°W) |
| * to the source CRS. Note that the longitude is set to the the center of |
| * the envelope longitude range (more on this below). |
| * |
| * 2) If the singularity point computed above is inside the source envelope, add that |
| * point to the target (transformed) envelope. |
| * |
| * 3) If step #2 added the point, iterate over all other axes. If an other bounded axis |
| * is found and that axis is of kind "WRAPAROUND", test for inclusion the same point |
| * than the point tested at step #1, except for the coordinate of the axis found in |
| * this step. That coordinate is set to the minimal and maximal values of that axis. |
| * |
| * Example: If the above steps found that the point (90°S, 30°W) need to be included, |
| * then this step #3 will also test the points (90°S, 180°W) and (90°S, 180°E). |
| * |
| * NOTE: we test (-180°, centerY), (180°, centerY), (centerX, -90°) and (centerX, 90°) |
| * at step #1 before to test (-180°, -90°), (180°, -90°), (-180°, 90°) and (180°, 90°) |
| * at step #3 because the later may not be supported by every projections. For example |
| * if the target envelope is located between 20°N and 40°N, then a Mercator projection |
| * may fail to transform the (-180°, 90°) coordinate while the (-180°, 30°) coordinate |
| * is a valid point. |
| */ |
| TransformException warning = null; |
| AbstractEnvelope generalEnvelope = null; |
| DirectPosition sourcePt = null; |
| DirectPosition targetPt = null; |
| long includedMinValue = 0; // A bitmask for each dimension. |
| long includedMaxValue = 0; |
| long isWrapAroundAxis = 0; |
| long dimensionBitMask = 1; |
| final int dimension = targetCS.getDimension(); |
| poles: for (int i=0; i<dimension; i++, dimensionBitMask <<= 1) { |
| final CoordinateSystemAxis axis = targetCS.getAxis(i); |
| if (axis == null) { // Should never be null, but check as a paranoiac safety. |
| continue; |
| } |
| boolean testMax = false; // Tells if we are testing the minimal or maximal value. |
| do { |
| final double extremum = testMax ? axis.getMaximumValue() : axis.getMinimumValue(); |
| if (!Double.isFinite(extremum)) { |
| /* |
| * The axis is unbounded. It should always be the case when the target CRS is |
| * a map projection, in which case this loop will finish soon and this method |
| * will do nothing more (no object instantiated, no MathTransform inversed...) |
| */ |
| continue; |
| } |
| if (targetPt == null) { |
| try { |
| mt = mt.inverse(); |
| } catch (NoninvertibleTransformException exception) { |
| /* |
| * If the transform is non invertible, this method can't do anything. This |
| * is not a fatal error because the envelope has already be transformed by |
| * the caller. We lost the check for singularity points performed by this |
| * method, but it make no difference in the common case where the source |
| * envelope didn't contains any of those points. |
| * |
| * Note that this exception is normal if target dimension is smaller than |
| * source dimension, since the math transform can not reconstituate the |
| * lost dimensions. So we don't log any warning in this case. |
| */ |
| if (dimension >= mt.getSourceDimensions()) { |
| warning = exception; |
| } |
| break poles; |
| } |
| targetPt = new GeneralDirectPosition(mt.getSourceDimensions()); |
| for (int j=0; j<dimension; j++) { |
| targetPt.setOrdinate(j, centerPt[j]); |
| } |
| // TODO: avoid the hack below if we provide a contains(DirectPosition) |
| // method in the GeoAPI org.opengis.geometry.Envelope interface. |
| generalEnvelope = AbstractEnvelope.castOrCopy(envelope); |
| } |
| targetPt.setOrdinate(i, extremum); |
| try { |
| sourcePt = mt.transform(targetPt, sourcePt); |
| } catch (TransformException exception) { |
| /* |
| * This exception may be normal. For example if may occur when projecting |
| * the latitude extremums with a cylindrical Mercator projection. Do not |
| * log any message (unless logging level is fine) and try the other points. |
| */ |
| if (warning == null) { |
| warning = exception; |
| } else { |
| warning.addSuppressed(exception); |
| } |
| continue; |
| } |
| if (generalEnvelope.contains(sourcePt)) { |
| transformed.add(targetPt); |
| if (testMax) includedMaxValue |= dimensionBitMask; |
| else includedMinValue |= dimensionBitMask; |
| } |
| } while ((testMax = !testMax) == true); |
| /* |
| * Keep trace of axes of kind WRAPAROUND, except if the two extremum values of that |
| * axis have been included in the envelope (in which case the next step after this |
| * loop doesn't need to be executed for that axis). |
| */ |
| if ((includedMinValue & includedMaxValue & dimensionBitMask) == 0 && CoordinateOperations.isWrapAround(axis)) { |
| isWrapAroundAxis |= dimensionBitMask; |
| } |
| // Restore 'targetPt' to its initial state, which is equals to 'centerPt'. |
| if (targetPt != null) { |
| targetPt.setOrdinate(i, centerPt[i]); |
| } |
| } |
| /* |
| * Step #3 described in the above "Algorithm" section: iterate over all dimensions |
| * of type "WRAPAROUND" for which minimal or maximal axis values have not yet been |
| * included in the envelope. The set of axes is specified by a bitmask computed in |
| * the above loop. We examine only the points that have not already been included |
| * in the envelope. |
| */ |
| final long includedBoundsValue = (includedMinValue | includedMaxValue); |
| if (includedBoundsValue != 0) { |
| while (isWrapAroundAxis != 0) { |
| final int wrapAroundDimension = Long.numberOfTrailingZeros(isWrapAroundAxis); |
| dimensionBitMask = 1 << wrapAroundDimension; |
| isWrapAroundAxis &= ~dimensionBitMask; // Clear now the bit, for the next iteration. |
| final CoordinateSystemAxis wrapAroundAxis = targetCS.getAxis(wrapAroundDimension); |
| final double min = wrapAroundAxis.getMinimumValue(); |
| final double max = wrapAroundAxis.getMaximumValue(); |
| /* |
| * Iterate over all axes for which a singularity point has been previously found, |
| * excluding the "wrap around axis" currently under consideration. |
| */ |
| for (long am=(includedBoundsValue & ~dimensionBitMask), bm; am != 0; am &= ~bm) { |
| bm = Long.lowestOneBit(am); |
| final int axisIndex = Long.numberOfTrailingZeros(bm); |
| final CoordinateSystemAxis axis = targetCS.getAxis(axisIndex); |
| /* |
| * switch (c) { |
| * case 0: targetPt = (..., singularityMin, ..., wrapAroundMin, ...) |
| * case 1: targetPt = (..., singularityMin, ..., wrapAroundMax, ...) |
| * case 2: targetPt = (..., singularityMax, ..., wrapAroundMin, ...) |
| * case 3: targetPt = (..., singularityMax, ..., wrapAroundMax, ...) |
| * } |
| */ |
| for (int c=0; c<4; c++) { |
| /* |
| * Set the coordinate value along the axis having the singularity point |
| * (cases c=0 and c=2). If the envelope did not included that point, |
| * then skip completely this case and the next one, i.e. skip c={0,1} |
| * or skip c={2,3}. |
| */ |
| double value = max; |
| if ((c & 1) == 0) { // 'true' if we are testing "wrapAroundMin". |
| if (((c == 0 ? includedMinValue : includedMaxValue) & bm) == 0) { |
| c++; // Skip also the case for "wrapAroundMax". |
| continue; |
| } |
| targetPt.setOrdinate(axisIndex, (c == 0) ? axis.getMinimumValue() : axis.getMaximumValue()); |
| value = min; |
| } |
| targetPt.setOrdinate(wrapAroundDimension, value); |
| try { |
| sourcePt = mt.transform(targetPt, sourcePt); |
| } catch (TransformException exception) { |
| if (warning == null) { |
| warning = exception; |
| } else { |
| warning.addSuppressed(exception); |
| } |
| continue; |
| } |
| if (generalEnvelope.contains(sourcePt)) { |
| transformed.add(targetPt); |
| } |
| } |
| targetPt.setOrdinate(axisIndex, centerPt[axisIndex]); |
| } |
| targetPt.setOrdinate(wrapAroundDimension, centerPt[wrapAroundDimension]); |
| } |
| } |
| /* |
| * At this point we finished envelope transformation. Verify if some coordinates 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. We do not wrap |
| * around if the source and target axes use the same range (e.g. the longitude stay [-180 … +180]°) in order |
| * to reduce the risk of discontinuities. If the user really wants unconditional wrap around, (s)he can call |
| * GeneralEnvelope.normalize(). |
| */ |
| final Set<Integer> wrapAroundChanges; |
| if (isOperationComplete && operation instanceof AbstractCoordinateOperation) { |
| wrapAroundChanges = ((AbstractCoordinateOperation) operation).getWrapAroundChanges(); |
| } else { |
| wrapAroundChanges = CoordinateOperations.wrapAroundChanges(sourceCRS, targetCS); |
| } |
| transformed.normalize(targetCS, 0, wrapAroundChanges.size(), wrapAroundChanges.iterator()); |
| if (warning != null) { |
| recoverableException(Envelopes.class, warning); |
| } |
| return transformed; |
| } |
| |
| /** |
| * Returns the bounding box of a geometry defined in <cite>Well Known Text</cite> (WKT) format. |
| * This method does not check the consistency of the provided WKT. For example it does not check |
| * that every points in a {@code LINESTRING} have the same dimension. However this method |
| * ensures that the parenthesis are balanced, in order to catch some malformed WKT. |
| * |
| * <p>Example:</p> |
| * <ul> |
| * <li>{@code BOX(-180 -90, 180 90)} (not really a geometry, but understood by many software products)</li> |
| * <li>{@code POINT(6 10)}</li> |
| * <li>{@code MULTIPOLYGON(((1 1, 5 1, 1 5, 1 1),(2 2, 3 2, 3 3, 2 2)))}</li> |
| * <li>{@code GEOMETRYCOLLECTION(POINT(4 6),LINESTRING(3 8,7 10))}</li> |
| * </ul> |
| * |
| * See {@link GeneralEnvelope#GeneralEnvelope(CharSequence)} for more information about the |
| * parsing rules. |
| * |
| * @param wkt the {@code BOX}, {@code POLYGON} or other kind of element to parse. |
| * @return the envelope of the given geometry. |
| * @throws FactoryException if the given WKT can not be parsed. |
| * |
| * @see #toString(Envelope) |
| * @see CRS#fromWKT(String) |
| * @see org.apache.sis.io.wkt |
| */ |
| public static Envelope fromWKT(final CharSequence wkt) throws FactoryException { |
| ArgumentChecks.ensureNonNull("wkt", wkt); |
| try { |
| return new GeneralEnvelope(wkt); |
| } catch (IllegalArgumentException e) { |
| throw new FactoryException(Errors.format( |
| Errors.Keys.UnparsableStringForClass_2, Envelope.class), e); |
| } |
| } |
| |
| /** |
| * Formats the given envelope as a {@code BOX} element. The output is like below, |
| * where <var>n</var> is the {@linkplain Envelope#getDimension() number of dimensions} |
| * (omitted if equals to 2): |
| * |
| * <blockquote>{@code BOX}<var>n</var>{@code D(}{@linkplain Envelope#getLowerCorner() lower |
| * corner}{@code ,} {@linkplain Envelope#getUpperCorner() upper corner}{@code )}</blockquote> |
| * |
| * <div class="note"><b>Note:</b> |
| * The {@code BOX} element is not part of the standard <cite>Well Known Text</cite> (WKT) format. |
| * However it is understood by many software libraries, for example GDAL and PostGIS.</div> |
| * |
| * The string returned by this method can be {@linkplain GeneralEnvelope#GeneralEnvelope(CharSequence) |
| * parsed} by the {@code GeneralEnvelope} constructor. |
| * |
| * @param envelope the envelope to format. |
| * @return this envelope as a {@code BOX} or {@code BOX3D} (most typical dimensions) element. |
| * |
| * @see #fromWKT(CharSequence) |
| * @see org.apache.sis.io.wkt |
| */ |
| public static String toString(final Envelope envelope) { |
| return AbstractEnvelope.toString(envelope, false); |
| } |
| |
| /** |
| * Formats the given envelope as a {@code POLYGON} element in the <cite>Well Known Text</cite> |
| * (WKT) format. {@code POLYGON} can be used as an alternative to {@code BOX} when the element |
| * needs to be considered as a standard WKT geometry. |
| * |
| * <p>The string returned by this method can be {@linkplain GeneralEnvelope#GeneralEnvelope(CharSequence) |
| * parsed} by the {@code GeneralEnvelope} constructor.</p> |
| * |
| * @param envelope the envelope to format. |
| * @return the envelope as a {@code POLYGON} in WKT format. |
| * @throws IllegalArgumentException if the given envelope can not be formatted. |
| * |
| * @see org.apache.sis.io.wkt |
| */ |
| public static String toPolygonWKT(final Envelope envelope) throws IllegalArgumentException { |
| /* |
| * Get the dimension, ignoring the trailing ones which have infinite values. |
| */ |
| int dimension = envelope.getDimension(); |
| while (dimension != 0) { |
| final double length = envelope.getSpan(dimension - 1); |
| if (Double.isFinite(length)) { |
| break; |
| } |
| dimension--; |
| } |
| if (dimension < 2) { |
| throw new IllegalArgumentException(Errors.format(Errors.Keys.EmptyEnvelope2D)); |
| } |
| final StringBuilder buffer = new StringBuilder("POLYGON("); |
| String separator = "("; |
| for (int corner = 0; corner < CORNERS.length; corner += 2) { |
| for (int i=0; i<dimension; i++) { |
| final double value; |
| switch (i) { |
| case 0: // Fall through |
| case 1: value = CORNERS[corner+i] ? envelope.getMaximum(i) : envelope.getMinimum(i); break; |
| default: value = envelope.getMedian(i); break; |
| } |
| trimFractionalPart(buffer.append(separator).append(value)); |
| separator = " "; |
| } |
| separator = ", "; |
| } |
| return buffer.append("))").toString(); |
| } |
| |
| /** |
| * Enumeration of the 4 corners in an envelope, with repetition of the first point. |
| * The values are (x,y) pairs with {@code false} meaning "minimal value" and {@code true} |
| * meaning "maximal value". This is used by {@link #toPolygonWKT(Envelope)} only. |
| */ |
| private static final boolean[] CORNERS = { |
| false, false, |
| false, true, |
| true, true, |
| true, false, |
| false, false |
| }; |
| } |
