sis-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From desruisse...@apache.org
Subject [sis] branch geoapi-4.0 updated: Deeper analysis of situations where a wraparound may occur during the transformation of a `GridGeometry` envelope. We must be careful to not apply too large translations on envelope corners: a lower-left corner should not become greater than upper-right corner for example.
Date Mon, 28 Sep 2020 18:48:14 GMT
This is an automated email from the ASF dual-hosted git repository.

desruisseaux pushed a commit to branch geoapi-4.0
in repository https://gitbox.apache.org/repos/asf/sis.git


The following commit(s) were added to refs/heads/geoapi-4.0 by this push:
     new ee73474  Deeper analysis of situations where a wraparound may occur during the transformation of a `GridGeometry` envelope. We must be careful to not apply too large translations on envelope corners: a lower-left corner should not become greater than upper-right corner for example.
ee73474 is described below

commit ee73474bb6d514e8567d7dab9a5617bbe95f2e43
Author: Martin Desruisseaux <martin.desruisseaux@geomatys.com>
AuthorDate: Mon Sep 28 20:40:45 2020 +0200

    Deeper analysis of situations where a wraparound may occur during the transformation of a `GridGeometry` envelope.
    We must be careful to not apply too large translations on envelope corners: a lower-left corner should not become greater than upper-right corner for example.
---
 .../coverage/grid/CoordinateOperationFinder.java   |  81 +++++----
 .../apache/sis/coverage/grid/GridDerivation.java   |  11 +-
 .../org/apache/sis/coverage/grid/GridExtent.java   |   6 +-
 .../org/apache/sis/coverage/grid/GridGeometry.java |   2 +-
 .../sis/coverage/grid/ResampledGridCoverage.java   |   2 +-
 .../sis/coverage/grid/GridDerivationTest.java      |  62 ++++++-
 .../internal/referencing/WraparoundInEnvelope.java | 200 +++++++++++++++++++++
 .../internal/referencing/WraparoundTransform.java  | 135 ++++++++++----
 .../internal/referencing/provider/Wraparound.java  |   2 +-
 .../referencing/WraparoundTransformTest.java       |   8 +-
 .../java/org/apache/sis/measure/Longitude.java     |   5 +-
 11 files changed, 427 insertions(+), 87 deletions(-)

diff --git a/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/CoordinateOperationFinder.java b/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/CoordinateOperationFinder.java
index f7afcb0..3c24f7b 100644
--- a/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/CoordinateOperationFinder.java
+++ b/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/CoordinateOperationFinder.java
@@ -66,7 +66,7 @@ import org.apache.sis.internal.util.Numerics;
  *   <li>Use the area of interest and grid resolution for refining the coordinate operation between two CRS.</li>
  * </ul>
  *
- * <b>Note:</b> except for the {@link #gridToGrid(PixelInCell)} convenience method,
+ * <b>Note:</b> except for the {@link #gridToGrid()} convenience method,
  * this class does not provide the complete chain of operations from grid to grid.
  * It provides only the operation from <em>cell indices</em> in source grid to <em>coordinates in the CRS</em>
  * of destination grid. Callers must add the last step (conversion from target CRS to cell indices) themselves.
@@ -180,7 +180,7 @@ final class CoordinateOperationFinder implements Supplier<double[]> {
     private boolean isWraparoundApplied;
 
     /**
-     * Creates a new finder.
+     * Creates a new finder initialized to {@link PixelInCell#CELL_CORNER} anchor.
      *
      * @param  source  the grid geometry which is the source of the coordinate operation to find.
      * @param  target  the grid geometry which is the target of the coordinate operation to find.
@@ -188,31 +188,7 @@ final class CoordinateOperationFinder implements Supplier<double[]> {
     CoordinateOperationFinder(final GridGeometry source, final GridGeometry target) {
         this.source = source;
         this.target = target;
-    }
-
-    /**
-     * Returns the CRS of the source grid geometry. If neither the source and target grid geometry
-     * define a CRS, then this method returns {@code null}.
-     *
-     * @throws IncompleteGridGeometryException if the target grid geometry has a CRS but the source
-     *         grid geometry has none. Note that the converse is allowed, in which case the target
-     *         CRS is assumed the same than the source.
-     */
-    private CoordinateReferenceSystem getSourceCRS() {
-        return source.isDefined(GridGeometry.CRS) ||
-               target.isDefined(GridGeometry.CRS) ? source.getCoordinateReferenceSystem() : null;
-    }
-
-    /**
-     * Returns the target of the "corner to CRS" transform.
-     * May be {@code null} if the neither the source and target grid geometry define a CRS.
-     *
-     * @throws IncompleteGridGeometryException if the target grid geometry has a CRS but the source
-     *         grid geometry has none. Note that the converse is allowed, in which case the target
-     *         CRS is assumed the same than the source.
-     */
-    final CoordinateReferenceSystem getTargetCRS() {
-        return (changeOfCRS != null) ? changeOfCRS.getTargetCRS() : getSourceCRS();
+        this.anchor = PixelInCell.CELL_CORNER;
     }
 
     /**
@@ -243,6 +219,31 @@ final class CoordinateOperationFinder implements Supplier<double[]> {
     }
 
     /**
+     * Returns the CRS of the source grid geometry. If neither the source and target grid geometry
+     * define a CRS, then this method returns {@code null}.
+     *
+     * @throws IncompleteGridGeometryException if the target grid geometry has a CRS but the source
+     *         grid geometry has none. Note that the converse is allowed, in which case the target
+     *         CRS is assumed the same than the source.
+     */
+    private CoordinateReferenceSystem getSourceCRS() {
+        return source.isDefined(GridGeometry.CRS) ||
+               target.isDefined(GridGeometry.CRS) ? source.getCoordinateReferenceSystem() : null;
+    }
+
+    /**
+     * Returns the target of the "corner to CRS" transform.
+     * May be {@code null} if the neither the source and target grid geometry define a CRS.
+     *
+     * @throws IncompleteGridGeometryException if the target grid geometry has a CRS but the source
+     *         grid geometry has none. Note that the converse is allowed, in which case the target
+     *         CRS is assumed the same than the source.
+     */
+    final CoordinateReferenceSystem getTargetCRS() {
+        return (changeOfCRS != null) ? changeOfCRS.getTargetCRS() : getSourceCRS();
+    }
+
+    /**
      * Returns the coordinate operation from source CRS to target CRS. It may be the identity operation.
      * We try to take envelopes in account because the operation choice may depend on the geographic area.
      *
@@ -295,14 +296,16 @@ final class CoordinateOperationFinder implements Supplier<double[]> {
      * Computes the transform from “grid coordinates of the source” to “grid coordinates of the target”.
      * This is a concatenation of {@link #gridToCRS()} with target "CRS to grid" transform.
      *
-     * @param  anchor  whether the operation is between cell centers or cell corners.
+     * <p><b>WARNING:</b> this method may return a mutable transform.
+     * That transform should be only short lived (e.g. just the time to transform an envelope).
+     * See {@link WraparoundTransform#transform(MathTransform, Envelope)}.</p>
+     *
      * @return operation from source grid indices to target grid indices.
      * @throws FactoryException if no operation can be found between the source and target CRS.
      * @throws TransformException if some coordinates can not be transformed to the specified target.
      * @throws IncompleteGridGeometryException if required CRS or a "grid to CRS" information is missing.
      */
-    final MathTransform gridToGrid(final PixelInCell anchor) throws FactoryException, TransformException {
-        setAnchor(anchor);
+    final MathTransform gridToGrid() throws FactoryException, TransformException {
         final MathTransform step1 = gridToCRS();
         final MathTransform step2 = target.getGridToCRS(anchor);
         if (step1.equals(step2)) {                                          // Optimization for a common case.
@@ -329,6 +332,10 @@ final class CoordinateOperationFinder implements Supplier<double[]> {
      *   <li>{@link #forwardChangeOfCRS} — cached for next invocation of this {@code gridToCRS()} method.</li>
      * </ul>
      *
+     * <p><b>WARNING:</b> this method may return a mutable transform.
+     * That transform should be only short lived (e.g. just the time to transform an envelope).
+     * See {@link WraparoundTransform#transform(MathTransform, Envelope)}.</p>
+     *
      * @return operation from source grid indices to target geospatial coordinates.
      * @throws FactoryException if no operation can be found between the source and target CRS.
      * @throws TransformException if some coordinates can not be transformed to the specified target.
@@ -350,13 +357,15 @@ final class CoordinateOperationFinder implements Supplier<double[]> {
                  */
 apply:          if (forwardChangeOfCRS == null) {
                     forwardChangeOfCRS = changeOfCRS.getMathTransform();
-                    DirectPosition median = median(target, null);
-                    if (median == null) {
-                        median = median(source, forwardChangeOfCRS);
-                        if (median == null) break apply;
+                    DirectPosition sourceMedian = median(source, forwardChangeOfCRS);
+                    DirectPosition targetMedian = median(target, null);
+                    if (targetMedian == null) {
+                        if (sourceMedian == null) break apply;
+                        targetMedian = sourceMedian;
+                        sourceMedian = null;
                     }
                     forwardChangeOfCRS = WraparoundTransform.forDomainOfUse(forwardChangeOfCRS,
-                                            changeOfCRS.getTargetCRS().getCoordinateSystem(), median);
+                            sourceMedian, targetMedian, changeOfCRS.getTargetCRS().getCoordinateSystem());
                 }
                 gridToCRS = MathTransforms.concatenate(gridToCRS, forwardChangeOfCRS);
             }
@@ -542,7 +551,7 @@ apply:          if (forwardChangeOfCRS == null) {
             if (median != null) {
                 final MathTransform inverseNoWrap = inverseChangeOfCRS;
                 inverseChangeOfCRS = WraparoundTransform.forDomainOfUse(inverseNoWrap,
-                        changeOfCRS().getSourceCRS().getCoordinateSystem(), median);
+                        null, median, changeOfCRS().getSourceCRS().getCoordinateSystem());
 
                 if (inverseChangeOfCRS != inverseNoWrap) {
                     crsToGrid = MathTransforms.concatenate(inverseChangeOfCRS, sourceCrsToGrid);
diff --git a/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/GridDerivation.java b/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/GridDerivation.java
index 9287864..8e6a3c4 100644
--- a/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/GridDerivation.java
+++ b/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/GridDerivation.java
@@ -422,18 +422,19 @@ public class GridDerivation {
         ensureSubgridNotSet();
         subGridSetter = "subgrid";
         if (!base.equals(gridOfInterest)) {
-            final MathTransform mapCorners, mapCenters;
-            final GridExtent domain = gridOfInterest.getExtent();                  // May throw IncompleteGridGeometryException.
+            final MathTransform mapCenters;
+            final GridExtent domain = gridOfInterest.getExtent();       // May throw IncompleteGridGeometryException.
             try {
                 final CoordinateOperationFinder finder = new CoordinateOperationFinder(gridOfInterest, base);
-                mapCorners = finder.gridToGrid(PixelInCell.CELL_CORNER);
-                mapCenters = finder.gridToGrid(PixelInCell.CELL_CENTER);
+                final MathTransform mapCorners = finder.gridToGrid();
+                finder.setAnchor(PixelInCell.CELL_CENTER);
+                mapCenters = finder.gridToGrid();
                 clipExtent(domain.toCRS(mapCorners, mapCenters, null));
             } catch (FactoryException | TransformException e) {
                 throw new IllegalGridGeometryException(e, "gridOfInterest");
             }
             if (baseExtent != base.extent && baseExtent.equals(gridOfInterest.extent)) {
-                baseExtent = gridOfInterest.extent;                                                 // Share common instance.
+                baseExtent = gridOfInterest.extent;                                         // Share common instance.
             }
             /*
              * The subsampling will be determined by scale factors of the transform from the given desired grid geometry to
diff --git a/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/GridExtent.java b/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/GridExtent.java
index 0e97c1b..c05ac7d 100644
--- a/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/GridExtent.java
+++ b/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/GridExtent.java
@@ -38,15 +38,15 @@ import org.apache.sis.util.ArgumentChecks;
 import org.apache.sis.util.resources.Errors;
 import org.apache.sis.util.resources.Vocabulary;
 import org.apache.sis.util.collection.WeakValueHashMap;
-import org.apache.sis.internal.referencing.ExtendedPrecisionMatrix;
 import org.apache.sis.internal.referencing.AxisDirections;
+import org.apache.sis.internal.referencing.ExtendedPrecisionMatrix;
+import org.apache.sis.internal.referencing.WraparoundTransform;
 import org.apache.sis.internal.feature.Resources;
 import org.apache.sis.internal.util.Numerics;
 import org.apache.sis.internal.util.Strings;
 import org.apache.sis.internal.util.DoubleDouble;
 import org.apache.sis.geometry.AbstractEnvelope;
 import org.apache.sis.geometry.GeneralEnvelope;
-import org.apache.sis.geometry.Envelopes;
 import org.apache.sis.coverage.SubspaceNotSpecifiedException;
 import org.apache.sis.referencing.operation.matrix.Matrices;
 import org.apache.sis.referencing.operation.matrix.MatrixSIS;
@@ -914,7 +914,7 @@ public class GridExtent implements GridEnvelope, LenientComparable, Serializable
             if (high != Long.MAX_VALUE) high++;             // Make the coordinate exclusive before cast.
             envelope.setRange(i, coordinates[i], high);     // Possible loss of precision in cast to `double` type.
         }
-        envelope = Envelopes.transform(cornerToCRS, envelope);
+        envelope = WraparoundTransform.transform(cornerToCRS, envelope);
         if (envelope.isEmpty()) try {
             /*
              * If the envelope contains some NaN values, try to replace them by constant values inferred from the math transform.
diff --git a/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/GridGeometry.java b/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/GridGeometry.java
index 28b2b8b..586bde7 100644
--- a/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/GridGeometry.java
+++ b/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/GridGeometry.java
@@ -1310,9 +1310,9 @@ public class GridGeometry implements LenientComparable, Serializable {
          * effort for using `WraparoundTransform` only if needed (contrarily to `gridToCRS(…)`).
          */
         final CoordinateOperationFinder finder = new CoordinateOperationFinder(target, this);
+        finder.setAnchor(anchor);
         final MathTransform tr;
         try {
-            finder.setAnchor(anchor);
             tr = finder.inverse(true);
         } catch (FactoryException e) {
             throw new TransformException(e);
diff --git a/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/ResampledGridCoverage.java b/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/ResampledGridCoverage.java
index 791f6c7..6fbe4c6 100644
--- a/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/ResampledGridCoverage.java
+++ b/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/ResampledGridCoverage.java
@@ -231,7 +231,7 @@ final class ResampledGridCoverage extends GridCoverage {
          * The following lines may throw IncompleteGridGeometryException, which is desired because if that
          * transform is missing, we can not continue (we have no way to guess it).
          */
-        changeOfCRS.setAnchor(PixelInCell.CELL_CORNER);
+        // Finder is initialized to PixelInCell.CELL_CORNER.
         final MathTransform sourceCornerToCRS = changeOfCRS.gridToCRS();
         final MathTransform crsToSourceCorner = changeOfCRS.inverse(true);
         changeOfCRS.setAnchor(PixelInCell.CELL_CENTER);
diff --git a/core/sis-feature/src/test/java/org/apache/sis/coverage/grid/GridDerivationTest.java b/core/sis-feature/src/test/java/org/apache/sis/coverage/grid/GridDerivationTest.java
index e1849f7..3e200fa 100644
--- a/core/sis-feature/src/test/java/org/apache/sis/coverage/grid/GridDerivationTest.java
+++ b/core/sis-feature/src/test/java/org/apache/sis/coverage/grid/GridDerivationTest.java
@@ -26,6 +26,9 @@ import org.opengis.referencing.datum.PixelInCell;
 import org.opengis.referencing.operation.MathTransform;
 import org.opengis.referencing.operation.TransformException;
 import org.apache.sis.geometry.Envelope2D;
+import org.apache.sis.geometry.GeneralEnvelope;
+import org.apache.sis.geometry.DirectPosition2D;
+import org.apache.sis.geometry.GeneralDirectPosition;
 import org.apache.sis.internal.referencing.Formulas;
 import org.apache.sis.internal.referencing.j2d.AffineTransform2D;
 import org.apache.sis.referencing.cs.AxesConvention;
@@ -36,9 +39,6 @@ import org.apache.sis.referencing.operation.matrix.Matrix4;
 import org.apache.sis.referencing.operation.transform.MathTransforms;
 import org.apache.sis.referencing.operation.HardCodedConversions;
 import org.apache.sis.referencing.crs.HardCodedCRS;
-import org.apache.sis.geometry.DirectPosition2D;
-import org.apache.sis.geometry.GeneralDirectPosition;
-import org.apache.sis.geometry.GeneralEnvelope;
 import org.apache.sis.test.DependsOnMethod;
 import org.apache.sis.test.DependsOn;
 import org.apache.sis.test.TestCase;
@@ -464,6 +464,62 @@ public final strictfp class GridDerivationTest extends TestCase {
     }
 
     /**
+     * Tests deriving a grid geometry when all involved grid geometries cross the anti-meridian.
+     * Illustration:
+     *
+     * {@preformat text
+     *   ──────────────┐                    ┌──────────────────
+     *        Grid     │                    │       Grid
+     *   ──────────────┘                    └──────────────────
+     *             120°W                    100°E
+     *   ─────────────────┐       ┌────────────────────────────
+     *          AOI       │       │      Area Of Interest
+     *   ─────────────────┘       └────────────────────────────
+     *                102°W       22°W
+     *   ↖…………………………………………………………………………………………………360° period…………↗︎
+     * }
+     */
+    @Test
+    public void testAntiMeridianCrossingInBothGrids() {
+        /*
+         * Longitudes from 100°E to 240°E (in WGS84 geographic CRS), which is equivalent to 100°E to 120°W.
+         * Latitude range is from 21°S to 60°N, but this is secondary for this test.
+         */
+        final GridGeometry grid = new GridGeometry(
+                new GridExtent(null, null, new long[] {8400, 4860}, true), PixelInCell.CELL_CENTER,
+                MathTransforms.linear(new Matrix3(
+                        0.016664682775860015,  0,  100.00833234138793,   0,
+                        0.016663238016868958,      -20.991668380991566,  0, 0, 1)),
+                HardCodedCRS.WGS84.forConvention(AxesConvention.POSITIVE_RANGE));
+        /*
+         * 22°27′34″W to 102°27′35″W in Mercator projection.
+         * Latitude range is about 66°S to 80°N, but it is not important for this test.
+         */
+        final GridGeometry areaOfInterest = new GridGeometry(
+                new GridExtent(null, null, new long[] {865, 725}, true), PixelInCell.CELL_CENTER,
+                MathTransforms.linear(new Matrix3(
+                         35992.44506018084,  0,  -2482193.2646860380,  0,
+                        -35992.44506018084,      16123175.875372937,   0, 0, 1)),
+                HardCodedConversions.mercator());
+        /*
+         * Verify that the area of interest contains the full grid.
+         */
+        final GeneralEnvelope e1 = new GeneralEnvelope(grid.getGeographicExtent().get());
+        final GeneralEnvelope e2 = new GeneralEnvelope(areaOfInterest.getGeographicExtent().get());
+        assertTrue(e2.contains(e1));
+        /*
+         * Expect the same longitude range than `grid` since it is fully included in `areaOfInterest`.
+         */
+        final GridGeometry result = grid.derive().subgrid(areaOfInterest).build();
+        assertEnvelopeEquals(new Envelope2D(null, 100, -21, 240 - 100, 60 - -21),
+                             result.getEnvelope(), Formulas.ANGULAR_TOLERANCE);
+        /*
+         * Following is empirical values provided as an anti-regression test.
+         */
+        assertExtentEquals(new long[2], new long[] {442, 285}, result.getExtent());
+    }
+
+    /**
      * Tests deriving a grid geometry from an area of interest shifted by 360° before or after the grid valid area.
      */
     @Test
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/WraparoundInEnvelope.java b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/WraparoundInEnvelope.java
new file mode 100644
index 0000000..ad3957c
--- /dev/null
+++ b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/WraparoundInEnvelope.java
@@ -0,0 +1,200 @@
+/*
+ * 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.internal.referencing;
+
+import org.apache.sis.internal.util.Numerics;
+import org.apache.sis.util.ComparisonMode;
+
+
+/**
+ * A {@link WraparoundTransform} where the number of cycles added or removed does not exceed a given limit.
+ * The bound is determined by whether the coordinate to transform is before or after a median point.
+ * If the coordinate is before the median, this class puts a limit on positive number of cycles added.
+ * If the coordinate is after  the median, this class puts a limit on positive number of cycles removed.
+ * The intent is to avoid that the lower bound of an envelope is shifted by a greater number of cycles
+ * than the upper bound, which may result in lower bound becoming greater than upper bound.
+ *
+ * The final result is that envelopes transformed using {@code WraparoundInEnvelope} may be larger
+ * than envelopes transformed using {@link WraparoundTransform} but should never be smaller.
+ *
+ * <h2>Mutability</h2>
+ * <b>This class is mutable.</b> This class records the translations that {@link #shift(double)} wanted to apply
+ * but could not because of the {@linkplain #limit} documented in above paragraph. When such missed translations
+ * are detected, caller should {@linkplain #translate() translate the limit} and transform same envelope again.
+ * We do that because each envelope transformation should use a consistent {@linkplain #limit} for all corners.
+ * Since this strategy breaks usual {@link org.apache.sis.referencing.operation.transform.AbstractMathTransform}
+ * contract about immutability, this class should be used only for temporary transforms to apply on an envelope.
+ *
+ * @author  Martin Desruisseaux (Geomatys)
+ * @version 1.1
+ * @since   1.1
+ * @module
+ */
+final class WraparoundInEnvelope extends WraparoundTransform {
+    /**
+     * For cross-version compatibility.
+     */
+    private static final long serialVersionUID = -1590996680159048327L;
+
+    /**
+     * The median source coordinate where the minimum or maximum number of cycles is determined.
+     * This <em>source</em> coordinate is not necessarily the same as the median <em>target</em> coordinate
+     * (which is set to zero by application of normalization matrix) because those medians were determined
+     * from source and target envelopes, which do not need to be the same.
+     */
+    private final double sourceMedian;
+
+    /**
+     * Number of cycles at the {@linkplain #sourceMedian} position. This is the minimum or maximum number
+     * of cycles to remove to a coordinate, depending if that coordinate is before or after the median.
+     */
+    private double limit;
+
+    /**
+     * The minimum and maximum number of {@linkplain #period period}s that {@link #shift(double)} wanted
+     * to add to the coordinate before to be constrained to the {@link #limit}.
+     */
+    private double minCycles, maxCycles;
+
+    /**
+     * Whether {@link #minCycles} or {@link #maxCycles} has been updated.
+     */
+    private boolean minChanged, maxChanged;
+
+    /**
+     * Creates a new transform with a wraparound behavior in the given dimension.
+     * Input and output values in the wraparound dimension shall be normalized in
+     * the [−p/2 … +p/2] range where <var>p</var> is the period (e.g. 360°).
+     */
+    WraparoundInEnvelope(final int dimension, final int wraparoundDimension,
+                         final double period, final double sourceMedian)
+    {
+        super(dimension, wraparoundDimension, period);
+        this.sourceMedian = sourceMedian;
+        reset();
+    }
+
+    /**
+     * Copy constructor for {@link #redim(int)}.
+     */
+    private WraparoundInEnvelope(final WraparoundInEnvelope source, final int dimension) {
+        super(dimension, source.wraparoundDimension, source.period);
+        sourceMedian = source.sourceMedian;
+        synchronized (source) {
+            minCycles = maxCycles = limit = source.limit;
+        }
+    }
+
+    /**
+     * Returns an instance with the specified number of dimensions while keeping
+     * {@link #wraparoundDimension} and {@link #period} unchanged.
+     */
+    @Override
+    WraparoundTransform redim(final int n) {
+        return new WraparoundInEnvelope(this, n);
+    }
+
+    /**
+     * Applies the wraparound on the given value. This implementation ensures that coordinates smaller than
+     * {@link #sourceMedian} before wraparound are still smaller than the (possibly shifted) median after wraparound,
+     * and conversely for coordinates greater than the median.
+     *
+     * The final result is that envelopes transformed using {@code WraparoundInEnvelope} may be larger
+     * than envelopes transformed using {@link WraparoundTransform} but should never be smaller.
+     *
+     * @param  x  the value on which to apply wraparound.
+     * @return the value after wraparound.
+     */
+    @Override
+    final synchronized double shift(final double x) {
+        double n = Math.rint(x / period);
+        if (x < sourceMedian) {
+            if (n < limit) {
+                if (n < minCycles) {
+                    minCycles = n;
+                    minChanged = true;
+                }
+                n = limit;
+            }
+        } else {
+            if (n > limit) {
+                if (n > maxCycles) {
+                    maxCycles = n;
+                    maxChanged = true;
+                }
+                n = limit;
+            }
+        }
+        return x - n * period;
+    }
+
+    /**
+     * Modifies this transform with a translation for enabling the wraparound that could not be applied in previous
+     * {@link #shift(double)} executions. If this method returns {@code true}, then this transform computes different
+     * output coordinates for the same input coordinates.
+     *
+     * <h4>Usage</h4>
+     * This method can be invoked after transforming an envelope. If this method returns {@code true}, then
+     * the same envelope should be transformed again and the new result added to the previous result (union).
+     *
+     * @return {@code true} if this transform has been modified.
+     */
+    final synchronized boolean translate() {
+        if (minChanged) {
+            minChanged = false;
+            limit = minCycles;
+            return true;
+        }
+        if (maxChanged) {
+            maxChanged = false;
+            limit = maxCycles;
+            return true;
+        }
+        return false;
+    }
+
+    /**
+     * Resets this transform to its initial state.
+     */
+    final synchronized void reset() {
+        minCycles = maxCycles = limit = Math.rint(sourceMedian / period);
+    }
+
+    /**
+     * Compares this transform with the given object for equality.
+     *
+     * @param  object  the object to compare with this transform.
+     * @param  mode    ignored, can be {@code null}.
+     */
+    @Override
+    public boolean equals(final Object object, final ComparisonMode mode) {
+        if (super.equals(object, mode)) {
+            final WraparoundInEnvelope other = (WraparoundInEnvelope) object;
+            return Numerics.equals(sourceMedian, other.sourceMedian) &&
+                   Numerics.equals(limit, other.limit);
+        }
+        return false;
+    }
+
+    /**
+     * Computes a hash code value for this transform.
+     */
+    @Override
+    protected int computeHashCode() {
+        return super.computeHashCode() + 7*Double.hashCode(limit);
+    }
+}
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/WraparoundTransform.java b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/WraparoundTransform.java
index 9b0573d..b58dedd 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/WraparoundTransform.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/WraparoundTransform.java
@@ -16,9 +16,10 @@
  */
 package org.apache.sis.internal.referencing;
 
-import java.io.Serializable;
 import java.util.List;
+import java.io.Serializable;
 import org.opengis.util.FactoryException;
+import org.opengis.geometry.Envelope;
 import org.opengis.geometry.DirectPosition;
 import org.opengis.parameter.ParameterValueGroup;
 import org.opengis.parameter.ParameterDescriptorGroup;
@@ -34,10 +35,13 @@ import org.apache.sis.referencing.operation.matrix.Matrices;
 import org.apache.sis.referencing.operation.transform.AbstractMathTransform;
 import org.apache.sis.referencing.operation.transform.MathTransforms;
 import org.apache.sis.internal.referencing.provider.Wraparound;
+import org.apache.sis.geometry.GeneralEnvelope;
+import org.apache.sis.geometry.Envelopes;
 import org.apache.sis.internal.system.Modules;
 import org.apache.sis.internal.util.Numerics;
 import org.apache.sis.measure.Longitude;
 import org.apache.sis.parameter.Parameters;
+import org.apache.sis.util.ArraysExt;
 import org.apache.sis.util.ArgumentChecks;
 import org.apache.sis.util.ComparisonMode;
 import org.apache.sis.util.logging.Logging;
@@ -46,9 +50,7 @@ import org.apache.sis.util.collection.BackingStoreException;
 
 /**
  * Enforces coordinate values in the range of a wraparound axis (typically longitude).
- * This transform is usually not needed for the [-180 … +180]° range since it is the
- * range of trigonometric functions. However this transform is useful for shifting
- * transformation results in the [0 … 360]° range.
+ * For example this transform can shift longitudes in the [-180 … +180]° range to the [0 … 360]° range.
  *
  * <p>{@code WraparoundTransform}s are not created automatically by {@link org.apache.sis.referencing.CRS#findOperation
  * CRS.findOperation(…)} because they introduce a discontinuity in coordinate transformations. Such discontinuities are
@@ -70,7 +72,7 @@ import org.apache.sis.util.collection.BackingStoreException;
  * @since   1.1
  * @module
  */
-public final class WraparoundTransform extends AbstractMathTransform implements Serializable {
+public class WraparoundTransform extends AbstractMathTransform implements Serializable {
     /**
      * For cross-version compatibility.
      */
@@ -104,12 +106,19 @@ public final class WraparoundTransform extends AbstractMathTransform implements
      * @param  dimension            number of dimensions of the transform to create.
      * @param  wraparoundDimension  dimension where wraparound happens.
      * @param  period               period on wraparound axis.
+     * @param  sourceMedian         coordinate in the center of source envelope, or 0 if none.
      * @return the wraparound transform (may be a cached instance).
      */
-    public static WraparoundTransform create(final int dimension, final int wraparoundDimension, final double period) {
+    public static WraparoundTransform create(final int dimension, final int wraparoundDimension,
+                                             final double period, final double sourceMedian)
+    {
         ArgumentChecks.ensureStrictlyPositive("dimension", dimension);
         ArgumentChecks.ensureBetween("wraparoundDimension", 0, dimension - 1, wraparoundDimension);
         ArgumentChecks.ensureStrictlyPositive("period", period);
+        if (sourceMedian != 0) {
+            ArgumentChecks.ensureFinite("sourceMedian", sourceMedian);
+            return new WraparoundInEnvelope(dimension, wraparoundDimension, period, sourceMedian);
+        }
         if (period == (Longitude.MAX_VALUE - Longitude.MIN_VALUE) && (wraparoundDimension & ~1) == 0) {
             final int i = ((dimension - FIRST_CACHED_DIMENSION) << 1) | wraparoundDimension;
             if (i >= 0 && i < CACHE.length) {
@@ -122,7 +131,7 @@ public final class WraparoundTransform extends AbstractMathTransform implements
     /**
      * The dimension of source and target coordinates.
      */
-    private final int dimension;
+    final int dimension;
 
     /**
      * The dimension where to apply wraparound.
@@ -141,20 +150,28 @@ public final class WraparoundTransform extends AbstractMathTransform implements
      * affine transforms which can efficiently be combined with other affine transforms, but it caused
      * more rounding errors.
      */
-    private final double period;
+    final double period;
 
     /**
      * Creates a new transform with a wraparound behavior in the given dimension.
      * Input and output values in the wraparound dimension shall be normalized in
      * the [−p/2 … +p/2] range where <var>p</var> is the period (e.g. 360°).
      */
-    private WraparoundTransform(final int dimension, final int wraparoundDimension, final double period) {
+    WraparoundTransform(final int dimension, final int wraparoundDimension, final double period) {
         this.dimension = dimension;
         this.wraparoundDimension = wraparoundDimension;
         this.period = period;
     }
 
     /**
+     * Returns an instance with the specified number of dimensions while keeping
+     * {@link #wraparoundDimension} and {@link #period} unchanged.
+     */
+    WraparoundTransform redim(final int n) {
+        return create(n, wraparoundDimension, period, 0);
+    }
+
+    /**
      * Returns an instance with the number of dimensions compatible with the given matrix, while keeping
      * {@link #wraparoundDimension} and {@link #period} unchanged.
      * If no instance can be created for the given number of dimensions, then this method returns {@code null}.
@@ -166,7 +183,7 @@ public final class WraparoundTransform extends AbstractMathTransform implements
         final int n = (applyOtherFirst ? other.getNumRow() : other.getNumCol()) - 1;
         if (n == dimension) return this;
         if (n >= wraparoundDimension) return null;
-        return create(n, wraparoundDimension, period);
+        return redim(n);
     }
 
     /**
@@ -186,7 +203,7 @@ public final class WraparoundTransform extends AbstractMathTransform implements
         MathTransform tr = op.getMathTransform();
         final CoordinateSystem targetCS = op.getTargetCRS().getCoordinateSystem();
         for (final int wraparoundDimension : CoordinateOperations.wrapAroundChanges(op)) {
-            tr = concatenate(tr, wraparoundDimension, targetCS, null);
+            tr = concatenate(tr, wraparoundDimension, null, null, targetCS);
         }
         return tr;
     }
@@ -202,18 +219,19 @@ public final class WraparoundTransform extends AbstractMathTransform implements
      * (e.g. when transforming a raster) and we want that output coordinates to be continuous
      * in that domain, independently of axis ranges.</p>
      *
-     * @param  tr        the transform to augment with "wrap around" behavior.
-     * @param  targetCS  the target coordinate system.
-     * @param  median    the coordinates to put at the center of new ranges.
+     * @param  tr            the transform to augment with "wrap around" behavior.
+     * @param  sourceMedian  the coordinates at the center of source ranges, or {@code null} if none.
+     * @param  targetMedian  the coordinates to keep at the center of new ranges.
+     * @param  targetCS      the target coordinate system.
      * @return the math transform with wraparound if needed.
      * @throws TransformException if a coordinate can not be computed.
      */
-    public static MathTransform forDomainOfUse(MathTransform tr, final CoordinateSystem targetCS,
-            final DirectPosition median) throws TransformException
+    public static MathTransform forDomainOfUse(MathTransform tr, final DirectPosition sourceMedian,
+            final DirectPosition targetMedian, final CoordinateSystem targetCS) throws TransformException
     {
         final int dimension = targetCS.getDimension();
         for (int i=0; i<dimension; i++) {
-            tr = concatenate(tr, i, targetCS, median);
+            tr = concatenate(tr, i, sourceMedian, targetMedian, targetCS);
         }
         return tr;
     }
@@ -226,25 +244,27 @@ public final class WraparoundTransform extends AbstractMathTransform implements
      *
      * @param  tr                   the transform to concatenate with a wraparound transform.
      * @param  wraparoundDimension  the dimension where "wrap around" behavior may apply.
-     * @param  targetCS             the target coordinate system.
-     * @param  median               the coordinate to put at the center of new range,
+     * @param  sourceMedian         the coordinate at the center of source envelope, or {@code null} if none.
+     * @param  targetMedian         the coordinate to put at the center of new range,
      *                              or {@code null} for standard axis center.
+     * @param  targetCS             the target coordinate system.
      * @return the math transform with "wrap around" behavior in the specified dimension.
      * @throws TransformException if a coordinate can not be computed.
      */
     private static MathTransform concatenate(final MathTransform tr, final int wraparoundDimension,
-            final CoordinateSystem targetCS, final DirectPosition median) throws TransformException
+            final DirectPosition sourceMedian, final DirectPosition targetMedian, final CoordinateSystem targetCS)
+            throws TransformException
     {
         final double period = WraparoundAdjustment.range(targetCS, wraparoundDimension);
         if (!(period > 0 && period != Double.POSITIVE_INFINITY)) {
             return tr;
         }
         double m;
-        if (median == null) {
+        if (targetMedian == null) {
             final CoordinateSystemAxis axis = targetCS.getAxis(wraparoundDimension);
             m = (axis.getMinimumValue() + axis.getMaximumValue()) / 2;
         } else try {
-            m = median.getOrdinate(wraparoundDimension);
+            m = targetMedian.getOrdinate(wraparoundDimension);
         } catch (BackingStoreException e) {
             // Some implementations compute coordinates only when first needed.
             throw e.unwrapOrRethrow(TransformException.class);
@@ -256,7 +276,7 @@ public final class WraparoundTransform extends AbstractMathTransform implements
              */
             final int power = 10 - Math.getExponent(m);
             m = Math.scalb(Math.rint(Math.scalb(m, power)), -power);
-        } else if (median == null) {
+        } else if (targetMedian == null) {
             /*
              * May happen if `WraparoundAdjustment.range(…)` recognized a longitude axis
              * despite the `CoordinateSystemAxis` not declarining minimum/maximum values.
@@ -267,8 +287,9 @@ public final class WraparoundTransform extends AbstractMathTransform implements
             // Invalid median value. Assume caller means "no wrap".
             return tr;
         }
+        final double sm = (sourceMedian != null) ? sourceMedian.getOrdinate(wraparoundDimension) - m : 0;
         final int dimension = tr.getTargetDimensions();
-        MathTransform wraparound = create(dimension, wraparoundDimension, period);
+        MathTransform wraparound = create(dimension, wraparoundDimension, period, sm);
         if (m != 0) {
             final double[] t = new double[dimension];
             t[wraparoundDimension] = m;
@@ -308,7 +329,7 @@ public final class WraparoundTransform extends AbstractMathTransform implements
      * @return the dimension of input points.
      */
     @Override
-    public int getSourceDimensions() {
+    public final int getSourceDimensions() {
         return dimension;
     }
 
@@ -318,7 +339,7 @@ public final class WraparoundTransform extends AbstractMathTransform implements
      * @return the dimension of output points.
      */
     @Override
-    public int getTargetDimensions() {
+    public final int getTargetDimensions() {
         return dimension;
     }
 
@@ -339,6 +360,16 @@ public final class WraparoundTransform extends AbstractMathTransform implements
     }
 
     /**
+     * Applies the wraparound on the given value.
+     *
+     * @param  x  the value on which to apply wraparound.
+     * @return the value after wraparound.
+     */
+    double shift(final double x) {
+        return Math.IEEEremainder(x, period);
+    }
+
+    /**
      * Wraparounds a single coordinate point in an array,
      * and optionally computes the transform derivative at that location.
      */
@@ -349,7 +380,7 @@ public final class WraparoundTransform extends AbstractMathTransform implements
         if (dstPts != null) {
             System.arraycopy(srcPts, srcOff, dstPts, dstOff, dimension);
             dstOff += wraparoundDimension;
-            dstPts[dstOff] = Math.IEEEremainder(dstPts[dstOff], period);
+            dstPts[dstOff] = shift(dstPts[dstOff]);
         }
         return derivate ? Matrices.createIdentity(dimension) : null;
     }
@@ -364,7 +395,7 @@ public final class WraparoundTransform extends AbstractMathTransform implements
         System.arraycopy(srcPts, srcOff, dstPts, dstOff, numPts * dimension);
         dstOff += wraparoundDimension;
         while (--numPts >= 0) {
-            dstPts[dstOff] = Math.IEEEremainder(dstPts[dstOff], period);
+            dstPts[dstOff] = shift(dstPts[dstOff]);
             dstOff += dimension;
         }
     }
@@ -379,7 +410,7 @@ public final class WraparoundTransform extends AbstractMathTransform implements
         System.arraycopy(srcPts, srcOff, dstPts, dstOff, numPts * dimension);
         dstOff += wraparoundDimension;
         while (--numPts >= 0) {
-            dstPts[dstOff] = (float) Math.IEEEremainder(dstPts[dstOff], period);
+            dstPts[dstOff] = (float) shift(dstPts[dstOff]);
             dstOff += dimension;
         }
     }
@@ -422,7 +453,7 @@ public final class WraparoundTransform extends AbstractMathTransform implements
          * If the other transform is also a `WraparoundTransform` for the same dimension,
          * then there is no need to concatenate those two consecutive redudant transforms.
          */
-        if (equals(other, null)) {
+        if (equals(other)) {
             return this;
         }
         final List<MathTransform> steps = MathTransforms.getSteps(other);
@@ -613,7 +644,7 @@ public final class WraparoundTransform extends AbstractMathTransform implements
      */
     @Override
     public boolean equals(final Object object, final ComparisonMode mode) {
-        if (object instanceof WraparoundTransform) {
+        if (super.equals(object, mode)) {
             final WraparoundTransform other = (WraparoundTransform) object;
             return other.dimension == dimension && other.wraparoundDimension == wraparoundDimension
                     && Numerics.equals(period, other.period);
@@ -628,4 +659,46 @@ public final class WraparoundTransform extends AbstractMathTransform implements
     protected int computeHashCode() {
         return dimension * 31 + wraparoundDimension + Double.hashCode(period);
     }
+
+    /**
+     * Transforms an envelope using the given math transform with special checks for wraparounds.
+     * The transformation is only approximated: the returned envelope may be bigger than necessary.
+     *
+     * <p>This may method modifies the given transform with translations for enabling wraparounds
+     * that could not be applied in previous {@link #shift(double)} executions.
+     * If the {@link WraparoundInEnvelope#translate()} method returns {@code true}, then the given
+     * transform will compute different output coordinates for the same input coordinates.</p>
+     *
+     * @param  transform  the transform to use.
+     * @param  envelope   envelope to transform. This envelope will not be modified.
+     * @return the transformed envelope.
+     * @throws TransformException if a transform failed.
+     */
+    public static GeneralEnvelope transform(final MathTransform transform, final Envelope envelope) throws TransformException {
+        final GeneralEnvelope result = Envelopes.transform(transform, envelope);
+        WraparoundInEnvelope[] wraparounds = null;
+        for (final MathTransform t : MathTransforms.getSteps(transform)) {
+            if (t instanceof WraparoundInEnvelope) {
+                if (wraparounds == null) {
+                    wraparounds = new WraparoundInEnvelope[] {(WraparoundInEnvelope) t};
+                } else {
+                    wraparounds = ArraysExt.append(wraparounds, (WraparoundInEnvelope) t);
+                }
+            }
+        }
+        if (wraparounds != null) {
+            for (;;) {
+                boolean done = false;
+                for (final WraparoundInEnvelope tr : wraparounds) {
+                    done |= tr.translate();
+                }
+                if (!done) break;
+                result.add(Envelopes.transform(transform, envelope));
+            }
+            for (final WraparoundInEnvelope tr : wraparounds) {
+                tr.reset();
+            }
+        }
+        return result;
+    }
 }
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/Wraparound.java b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/Wraparound.java
index 26f36a3..34c95fc 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/Wraparound.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/Wraparound.java
@@ -115,6 +115,6 @@ public final class Wraparound extends AbstractProvider {
         return WraparoundTransform.create(
                 pg.intValue(DIMENSION),
                 pg.intValue(WRAPAROUND_DIMENSION),
-                pg.doubleValue(PERIOD));
+                pg.doubleValue(PERIOD), 0);
     }
 }
diff --git a/core/sis-referencing/src/test/java/org/apache/sis/internal/referencing/WraparoundTransformTest.java b/core/sis-referencing/src/test/java/org/apache/sis/internal/referencing/WraparoundTransformTest.java
index dcdf148..59fb0e2 100644
--- a/core/sis-referencing/src/test/java/org/apache/sis/internal/referencing/WraparoundTransformTest.java
+++ b/core/sis-referencing/src/test/java/org/apache/sis/internal/referencing/WraparoundTransformTest.java
@@ -49,10 +49,10 @@ public final strictfp class WraparoundTransformTest extends TestCase {
     public void testCache() {
         final double period = 360;
         final WraparoundTransform t1, t2, t3, t4;
-        assertSame   (WraparoundTransform.create(3, 0, period), t1 = WraparoundTransform.create(3, 0, period));
-        assertNotSame(WraparoundTransform.create(3, 0, period), t2 = WraparoundTransform.create(3, 1, period));
-        assertNotSame(WraparoundTransform.create(3, 0, period), t3 = WraparoundTransform.create(2, 0, period));
-        assertNotSame(WraparoundTransform.create(3, 0, period), t4 = WraparoundTransform.create(3, 2, period));
+        assertSame   (WraparoundTransform.create(3, 0, period, 0), t1 = WraparoundTransform.create(3, 0, period, 0));
+        assertNotSame(WraparoundTransform.create(3, 0, period, 0), t2 = WraparoundTransform.create(3, 1, period, 0));
+        assertNotSame(WraparoundTransform.create(3, 0, period, 0), t3 = WraparoundTransform.create(2, 0, period, 0));
+        assertNotSame(WraparoundTransform.create(3, 0, period, 0), t4 = WraparoundTransform.create(3, 2, period, 0));
         assertEquals(3, t1.getSourceDimensions());
         assertEquals(3, t2.getSourceDimensions());
         assertEquals(2, t3.getSourceDimensions());
diff --git a/core/sis-utility/src/main/java/org/apache/sis/measure/Longitude.java b/core/sis-utility/src/main/java/org/apache/sis/measure/Longitude.java
index 96a0d35..2e03d63 100644
--- a/core/sis-utility/src/main/java/org/apache/sis/measure/Longitude.java
+++ b/core/sis-utility/src/main/java/org/apache/sis/measure/Longitude.java
@@ -18,6 +18,7 @@ package org.apache.sis.measure;
 
 import org.opengis.geometry.DirectPosition;
 import org.opengis.referencing.cs.AxisDirection;
+import org.apache.sis.math.MathFunctions;
 
 
 /**
@@ -176,7 +177,7 @@ public final class Longitude extends Angle {
      * @since 1.1
      */
     public static boolean isWraparound(final double west, final double east) {
-        return (west > east) || (Double.doubleToRawLongBits(west) == 0 &&
-                                 Double.doubleToRawLongBits(east) == Long.MIN_VALUE);
+        return (west > east) || (MathFunctions.isPositiveZero(west) &&
+                                 MathFunctions.isNegativeZero(east));
     }
 }


Mime
View raw message