sis-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From desruisse...@apache.org
Subject [sis] 02/02: Initial support of Area Of Interest (AOI) crossing anti-meridian in call to GridDerivation.subgrid(…).
Date Mon, 25 Feb 2019 20:17:46 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

commit 0d35997a981ba75463065f10741e5daa6addabfd
Author: Martin Desruisseaux <martin.desruisseaux@geomatys.com>
AuthorDate: Mon Feb 25 21:17:12 2019 +0100

    Initial support of Area Of Interest (AOI) crossing anti-meridian in call to GridDerivation.subgrid(…).
---
 .../apache/sis/coverage/grid/GridDerivation.java   | 144 +++++++++++++++++++--
 .../org/apache/sis/coverage/grid/GridGeometry.java |   1 +
 .../sis/coverage/grid/GridDerivationTest.java      |   9 +-
 .../org/apache/sis/geometry/GeneralEnvelope.java   |   4 +-
 .../sis/internal/referencing/ExtentSelector.java   |   2 +-
 .../internal/referencing/ReferencingUtilities.java |  34 +++++
 .../referencing/ReferencingUtilitiesTest.java      |  13 +-
 .../src/main/java/org/apache/sis/math/Vector.java  |   2 +-
 8 files changed, 188 insertions(+), 21 deletions(-)

diff --git a/core/sis-raster/src/main/java/org/apache/sis/coverage/grid/GridDerivation.java
b/core/sis-raster/src/main/java/org/apache/sis/coverage/grid/GridDerivation.java
index c79bef1..4e811a3 100644
--- a/core/sis-raster/src/main/java/org/apache/sis/coverage/grid/GridDerivation.java
+++ b/core/sis-raster/src/main/java/org/apache/sis/coverage/grid/GridDerivation.java
@@ -28,13 +28,16 @@ import org.opengis.referencing.operation.TransformException;
 import org.opengis.referencing.operation.CoordinateOperation;
 import org.opengis.referencing.operation.NoninvertibleTransformException;
 import org.opengis.referencing.crs.CoordinateReferenceSystem;
+import org.opengis.referencing.cs.CoordinateSystem;
 import org.apache.sis.referencing.operation.transform.MathTransforms;
 import org.apache.sis.referencing.operation.transform.TransformSeparator;
 import org.apache.sis.referencing.operation.matrix.Matrices;
 import org.apache.sis.referencing.CRS;
+import org.apache.sis.internal.referencing.ReferencingUtilities;
 import org.apache.sis.internal.referencing.DirectPositionView;
 import org.apache.sis.geometry.GeneralDirectPosition;
 import org.apache.sis.geometry.GeneralEnvelope;
+import org.apache.sis.geometry.AbstractEnvelope;
 import org.apache.sis.geometry.Envelopes;
 import org.apache.sis.internal.raster.Resources;
 import org.apache.sis.util.resources.Vocabulary;
@@ -47,6 +50,7 @@ import org.apache.sis.util.Debug;
 import org.apache.sis.util.collection.DefaultTreeTable;
 import org.apache.sis.util.collection.TableColumn;
 import org.apache.sis.util.collection.TreeTable;
+import org.apache.sis.math.MathFunctions;
 
 // Branch-dependent imports
 import org.opengis.coverage.PointOutsideCoverageException;
@@ -431,6 +435,123 @@ public class GridDerivation {
     }
 
     /**
+     * Returns an envelope with coordinates equivalent to the given coordinates,
+     * but potentially shifted for intersecting the envelope of the grid geometry.
+     * This method assumes that {@code base.envelope} is non-null, which should be
+     * the case if {@link GridGeometry#requireGridToCRS()} has been invoked first.
+     *
+     * @see GeneralEnvelope#simplify()
+     */
+    private Envelope adjustWraparoundAxes(final Envelope areaOfInterest, final CoordinateOperation
baseToAOI) throws TransformException {
+        CoordinateReferenceSystem crs = areaOfInterest.getCoordinateReferenceSystem();
+        if (crs == null) {
+            crs = base.envelope.getCoordinateReferenceSystem();
+        }
+        if (crs != null) {
+            GeneralEnvelope  shifted  = null;
+            AbstractEnvelope validity = null;
+            final DirectPosition lowerCorner = areaOfInterest.getLowerCorner();
+            final DirectPosition upperCorner = areaOfInterest.getUpperCorner();
+            final CoordinateSystem cs = crs.getCoordinateSystem();
+            for (int i=cs.getDimension(); --i >= 0;) {
+                final double period = ReferencingUtilities.getWraparoundRange(cs, i);
+                if (period > 0) {
+                    /*
+                     * Found an axis (typically the longitude axis) with wraparound range
meaning.
+                     * We are going to need grid geometry envelope in the same CRS than the
AOI.
+                     * Compute that envelope when first needed.
+                     */
+                    if (validity == null) {
+                        validity = base.envelope;
+                        if (baseToAOI != null) {
+                            final MathTransform mt = baseToAOI.getMathTransform();
+                            if (!mt.isIdentity()) {
+                                validity = Envelopes.transform(mt, validity);
+                            }
+                        }
+                    }
+                    /*
+                     * "Unroll" the range. For example if we have [+160 … -170]° of longitude,
we can replace by [160 … 190]°.
+                     * We do not change the 'lower' or 'upper' value now in order to avoid
rounding error. Instead we compute
+                     * how many periods we need to add to those values. We adjust the side
which results in the value closest
+                     * to zero, in order to reduce rounding error if no more adjustment is
done in the next block.
+                     */
+                    final double lower = lowerCorner.getOrdinate(i);
+                    final double upper = upperCorner.getOrdinate(i);
+                    double lowerCycles = 0;                             // In number of periods.
+                    double upperCycles = 0;
+                    double delta = upper - lower;
+                    if (MathFunctions.isNegative(delta)) {              // Use 'isNegative'
for catching [+0 … -0] range.
+                        final double cycles = (delta == 0) ? -1 : Math.floor(delta / period);
        // Always negative.
+                        delta = cycles * period;
+                        if (Math.abs(lower + delta) < Math.abs(upper - delta)) {
+                            lowerCycles = cycles;                                    // Will
subtract periods to 'lower'.
+                        } else {
+                            upperCycles = -cycles;                                   // Will
add periods to 'upper'.
+                        }
+                    }
+                    /*
+                     * The range may be before or after the range of the grid geometry. Compute
the distance from current
+                     * lower/upper coordinate to the coordinate of valid area. The sign tells
us whether we are inside or
+                     * outside valid area. If both points are inside, there is nothing to
do. If both points are outside,
+                     * we can do nothing for now and clip to the valid area later. If one
point is inside while the other
+                     * is outside, try to shift the range for moving that point inside (or
at least make it closer).
+                     *
+                     * We try to compute multiples of 'periods' instead than just adding
or subtracting 'periods' once in
+                     * order to support images that cover more than one period, for example
images over 720° of longitude.
+                     * It may happen for example if an image shows data below the trajectory
of a satellite.
+                     */
+                    final double  lowerToValid   = Math.ceil ((validity.getLower(i) - lower)
/ period) - lowerCycles;
+                    final double  upperToValid   = Math.floor((validity.getUpper(i) - upper)
/ period) - upperCycles;
+                    final boolean lowerIsOutside = (lowerToValid > 0);
+                    final boolean upperIsOutside = (upperToValid < 0);
+                    if (lowerIsOutside != upperIsOutside) {
+                        final double cycles;
+                        if (lowerIsOutside) {
+                            /*
+                             * We need to add 'lowerToValid' to both sides in order to move
the range inside the valid area.
+                             * The new distances (in numbers of period) to valid area after
the shift will be:
+                             *
+                             *     upperToValid -= lowerToValid;
+                             *     lowerToValid -= lowerToValid;        // == 0
+                             *
+                             * If the shift causes the upper value to go outside the valid
area (new upperToValid < 0), it may
+                             * change which side of the valid area is intersected by the
AOI. For example before the shift the
+                             * AOI may intersect the left part of the valid area, but after
the shift AOI would intersect the
+                             * right part of the valid area. In order to avoid that change,
we have to stop the shift before
+                             * the upper value goes outside valid area.
+                             */
+                            cycles = Math.min(lowerToValid, upperToValid);              
   // All values here are positive.
+                        } else {
+                            /*
+                             * We need to subtract 'upperToValid' to both sides in order
to move the range inside the valid area.
+                             * The same argument than above apply, with sign reverted and
lower/upper variables interchanged.
+                             */
+                            cycles = Math.max(lowerToValid, upperToValid);              
   // All values here are negative.
+                        }
+                        lowerCycles += cycles;
+                        upperCycles += cycles;
+                    }
+                    /*
+                     * If there is change to apply, copy the envelope when first needed.
+                     */
+                    if (lowerCycles != 0 || upperCycles != 0) {
+                        if (shifted == null) {
+                            shifted = new GeneralEnvelope(areaOfInterest);
+                        }
+                        shifted.setRange(i, lower + lowerCycles * period,
+                                            upper + upperCycles * period);
+                    }
+                }
+            }
+            if (shifted != null) {
+                return shifted;
+            }
+        }
+        return areaOfInterest;
+    }
+
+    /**
      * Requests a grid geometry over a sub-region of the base grid geometry and optionally
with subsampling.
      * The given envelope does not need to be expressed in the same coordinate reference
system (CRS)
      * than {@linkplain GridGeometry#getCoordinateReferenceSystem() the CRS of the base grid
geometry};
@@ -465,7 +586,7 @@ public class GridDerivation {
      *
      * @see GridExtent#subsample(int[])
      */
-    public GridDerivation subgrid(final Envelope areaOfInterest, double... resolution) {
+    public GridDerivation subgrid(Envelope areaOfInterest, double... resolution) {
         ensureSubgridNotSet();
         MathTransform cornerToCRS = base.requireGridToCRS();
         subGridSetter = "subgrid";
@@ -475,9 +596,9 @@ public class GridDerivation {
              * to the 'gridToCRS' transform.  We should not transform the envelope here -
only concatenate the
              * transforms - because transforming envelopes twice would add errors.
              */
-            final CoordinateOperation operation = Envelopes.findOperation(base.envelope,
areaOfInterest);
-            if (operation != null) {
-                cornerToCRS = MathTransforms.concatenate(cornerToCRS, operation.getMathTransform());
+            final CoordinateOperation baseToAOI = Envelopes.findOperation(base.envelope,
areaOfInterest);
+            if (baseToAOI != null) {
+                cornerToCRS = MathTransforms.concatenate(cornerToCRS, baseToAOI.getMathTransform());
             }
             /*
              * If the envelope dimensions does not encompass all grid dimensions, the envelope
is probably non-invertible.
@@ -494,6 +615,7 @@ public class GridDerivation {
             dimension = baseExtent.getDimension();      // Non-null since 'base.requireGridToCRS()'
succeed.
             GeneralEnvelope indices = null;
             if (areaOfInterest != null) {
+                areaOfInterest = adjustWraparoundAxes(areaOfInterest, baseToAOI);
                 indices = Envelopes.transform(cornerToCRS.inverse(), areaOfInterest);
                 clipExtent(indices);
             }
@@ -521,11 +643,11 @@ public class GridDerivation {
             if (resolution != null && resolution.length != 0) {
                 resolution = ArraysExt.resize(resolution, cornerToCRS.getTargetDimensions());
                 Matrix m = cornerToCRS.derivative(new DirectPositionView.Double(getPointOfInterest()));
-                resolution = Matrices.inverse(m).multiply(resolution);
+                final double[] subsampling = Matrices.inverse(m).multiply(resolution);
                 final int[] modifiedDimensions = this.modifiedDimensions;               
     // Will not change anymore.
                 boolean modified = false;
-                for (int k=0; k<resolution.length; k++) {
-                    double s = Math.abs(resolution[k]);
+                for (int k=0; k<subsampling.length; k++) {
+                    double s = Math.abs(subsampling[k]);
                     if (s > 1) {                                // Also for skipping NaN
values.
                         final int i = (modifiedDimensions != null) ? modifiedDimensions[k]
: k;
                         final int accuracy = Math.max(0, Math.getExponent(indices.getSpan(i)))
+ 1;         // Power of 2.
@@ -534,7 +656,7 @@ public class GridDerivation {
                                             indices.getUpper(i) / s);
                         modified = true;
                     }
-                    resolution[k] = s;
+                    subsampling[k] = s;
                 }
                 /*
                  * If at least one subsampling is effective, build a scale from the old grid
coordinates to the new
@@ -547,8 +669,8 @@ public class GridDerivation {
                     scaledExtent = new GridExtent(indices, rounding, null, null, modifiedDimensions);
                     if (baseExtent.equals(scaledExtent)) scaledExtent = baseExtent;
                     m = Matrices.createIdentity(dimension + 1);
-                    for (int k=0; k<resolution.length; k++) {
-                        final double s = resolution[k];
+                    for (int k=0; k<subsampling.length; k++) {
+                        final double s = subsampling[k];
                         if (s > 1) {                            // Also for skipping NaN
values.
                             final int i = (modifiedDimensions != null) ? modifiedDimensions[k]
: k;
                             m.setElement(i, i, s);
@@ -556,7 +678,7 @@ public class GridDerivation {
                         }
                     }
                     toBase = MathTransforms.linear(m);
-                    scales = resolution;                        // For information purpose
only.
+                    scales = subsampling;                       // For information purpose
only.
                 }
             }
         } catch (FactoryException | TransformException e) {
diff --git a/core/sis-raster/src/main/java/org/apache/sis/coverage/grid/GridGeometry.java
b/core/sis-raster/src/main/java/org/apache/sis/coverage/grid/GridGeometry.java
index 150499e..2c48f0e 100644
--- a/core/sis-raster/src/main/java/org/apache/sis/coverage/grid/GridGeometry.java
+++ b/core/sis-raster/src/main/java/org/apache/sis/coverage/grid/GridGeometry.java
@@ -935,6 +935,7 @@ public class GridGeometry implements Serializable {
     /**
      * Verifies that this grid geometry defines an {@linkplain #extent} and a {@link #cornerToCRS}
transform.
      * They are the information required for mapping the grid to a spatiotemporal envelope.
+     * Note that this implies that {@link #envelope} is non-null.
      *
      * @return {@link #cornerToCRS}.
      */
diff --git a/core/sis-raster/src/test/java/org/apache/sis/coverage/grid/GridDerivationTest.java
b/core/sis-raster/src/test/java/org/apache/sis/coverage/grid/GridDerivationTest.java
index 7752630..13588d5 100644
--- a/core/sis-raster/src/test/java/org/apache/sis/coverage/grid/GridDerivationTest.java
+++ b/core/sis-raster/src/test/java/org/apache/sis/coverage/grid/GridDerivationTest.java
@@ -32,7 +32,6 @@ import org.apache.sis.test.DependsOnMethod;
 import org.apache.sis.test.DependsOn;
 import org.apache.sis.test.TestCase;
 import org.junit.Test;
-import org.junit.Ignore;
 
 import static org.apache.sis.test.ReferencingAssert.*;
 import static org.apache.sis.coverage.grid.GridGeometryTest.assertExtentEquals;
@@ -255,7 +254,6 @@ public final strictfp class GridDerivationTest extends TestCase {
      * Tests deriving a grid geometry with an envelope crossing the antimeridian.
      */
     @Test
-    @Ignore("TODO: not yet fixed.")
     public void testSubgridCrossingAntiMeridian() {
         final GridGeometry grid = new GridGeometry(
                 new GridExtent(200, 180), PixelInCell.CELL_CORNER,
@@ -265,11 +263,12 @@ public final strictfp class GridDerivationTest extends TestCase {
                         0,  0,  1)), HardCodedCRS.WGS84);
 
         final GeneralEnvelope aoi = new GeneralEnvelope(HardCodedCRS.WGS84);
-        aoi.setRange(0, 140, -179);
+        aoi.setRange(0, 140, -179);                 // Cross anti-meridian.
         aoi.setRange(1, -90,   90);
 
         final GridGeometry subgrid = grid.derive().subgrid(aoi).build();
-        Envelope subEnv = subgrid.getEnvelope();
-        assertEquals(aoi, subEnv);
+        final Envelope subEnv = subgrid.getEnvelope();
+        aoi.setRange(0, 140, 181);
+        assertEnvelopeEquals(aoi, subEnv);
     }
 }
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/geometry/GeneralEnvelope.java
b/core/sis-referencing/src/main/java/org/apache/sis/geometry/GeneralEnvelope.java
index 195fc27..f4d7a9b 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/geometry/GeneralEnvelope.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/geometry/GeneralEnvelope.java
@@ -973,7 +973,7 @@ public class GeneralEnvelope extends ArrayEnvelope implements Cloneable,
Seriali
     }
 
     /**
-     * Ensures that <var>lower</var> &lt;= <var>upper</var> for
every dimensions.
+     * Ensures that <var>lower</var> ≦ <var>upper</var> for every
dimensions.
      * If a {@linkplain #getUpper(int) upper ordinate value} is less than a
      * {@linkplain #getLower(int) lower ordinate value}, then there is a choice:
      *
@@ -1005,7 +1005,7 @@ public class GeneralEnvelope extends ArrayEnvelope implements Cloneable,
Seriali
             final int iUpper = iLower + d;
             final double lower = ordinates[iLower];
             final double upper = ordinates[iUpper];
-            if (isNegative(upper - lower)) {
+            if (isNegative(upper - lower)) {                            // Use 'isNegative'
for catching [+0 … -0] range.
                 final CoordinateSystemAxis axis = getAxis(crs, i);
                 if (isWrapAround(axis)) {
                     ordinates[iLower] = axis.getMinimumValue();
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/ExtentSelector.java
b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/ExtentSelector.java
index c6b0e9c..6201282 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/ExtentSelector.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/ExtentSelector.java
@@ -31,7 +31,7 @@ import org.apache.sis.metadata.iso.extent.Extents;
  *
  * @param <T>  the type of object to be selected.
  *
- * @since   0.4
+ * @since 0.4
  * @module
  */
 public final class ExtentSelector<T> {
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/ReferencingUtilities.java
b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/ReferencingUtilities.java
index a5441de..1663be8 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/ReferencingUtilities.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/ReferencingUtilities.java
@@ -40,6 +40,9 @@ import org.opengis.referencing.datum.VerticalDatumType;
 import org.opengis.referencing.operation.CoordinateOperationFactory;
 import org.opengis.util.FactoryException;
 import org.apache.sis.internal.system.DefaultFactories;
+import org.apache.sis.internal.metadata.AxisDirections;
+import org.apache.sis.measure.Longitude;
+import org.apache.sis.measure.Units;
 import org.apache.sis.util.Static;
 import org.apache.sis.util.Utilities;
 import org.apache.sis.util.CharSequences;
@@ -96,6 +99,37 @@ public final class ReferencingUtilities extends Static {
     }
 
     /**
+     * Returns the range (maximum - minimum) of the given axis if it has wraparound meaning,
+     * or {@link Double#NaN} otherwise. This method implements a fallback for the longitude
+     * axis if it does not declare the minimum and maximum values as expected.
+     *
+     * @param  cs         the coordinate system for which to get wraparound range, or {@code
null}.
+     * @param  dimension  dimension of the axis to test.
+     * @return the wraparound range, or {@link Double#NaN} if none.
+     */
+    public static double getWraparoundRange(final CoordinateSystem cs, final int dimension)
{
+        if (cs != null) {
+            final CoordinateSystemAxis axis = cs.getAxis(dimension);
+            if (axis != null && RangeMeaning.WRAPAROUND.equals(axis.getRangeMeaning()))
{
+                double period = axis.getMaximumValue() - axis.getMinimumValue();
+                if (period > 0 && period != Double.POSITIVE_INFINITY) {
+                    return period;
+                }
+                final AxisDirection dir = AxisDirections.absolute(axis.getDirection());
+                if (AxisDirection.EAST.equals(dir) && cs instanceof EllipsoidalCS)
{
+                    period = Longitude.MAX_VALUE - Longitude.MIN_VALUE;
+                    final Unit<?> unit = axis.getUnit();
+                    if (unit != null) {
+                        period = Units.DEGREE.getConverterTo(Units.ensureAngular(unit)).convert(period);
+                    }
+                    return period;
+                }
+            }
+        }
+        return Double.NaN;
+    }
+
+    /**
      * Returns the unit used for all axes in the given coordinate system.
      * If not all axes use the same unit, then this method returns {@code null}.
      *
diff --git a/core/sis-referencing/src/test/java/org/apache/sis/internal/referencing/ReferencingUtilitiesTest.java
b/core/sis-referencing/src/test/java/org/apache/sis/internal/referencing/ReferencingUtilitiesTest.java
index 2893598..c8be532 100644
--- a/core/sis-referencing/src/test/java/org/apache/sis/internal/referencing/ReferencingUtilitiesTest.java
+++ b/core/sis-referencing/src/test/java/org/apache/sis/internal/referencing/ReferencingUtilitiesTest.java
@@ -25,6 +25,7 @@ import org.opengis.referencing.datum.VerticalDatum;
 import org.opengis.referencing.IdentifiedObject;
 import org.apache.sis.referencing.datum.HardCodedDatum;
 import org.apache.sis.referencing.crs.HardCodedCRS;
+import org.apache.sis.referencing.cs.HardCodedCS;
 import org.apache.sis.util.Utilities;
 import org.apache.sis.measure.Units;
 import org.apache.sis.test.TestCase;
@@ -38,7 +39,7 @@ import static org.apache.sis.internal.referencing.ReferencingUtilities.*;
  * Tests {@link ReferencingUtilities}.
  *
  * @author  Martin Desruisseaux (Geomatys)
- * @version 0.8
+ * @version 1.0
  * @since   0.4
  * @module
  */
@@ -57,6 +58,16 @@ public final strictfp class ReferencingUtilitiesTest extends TestCase {
     }
 
     /**
+     * Tests {@link ReferencingUtilities#getWraparoundRange(CoordinateSystem, int)}.
+     */
+    @Test
+    public void testGetWraparoundRange() {
+        assertTrue  (Double.isNaN(getWraparoundRange(HardCodedCS.GEODETIC_φλ, 0)));
+        assertEquals(360, getWraparoundRange(HardCodedCS.GEODETIC_φλ, 1), STRICT);
+        assertEquals(400, getWraparoundRange(HardCodedCS.ELLIPSOIDAL_gon, 0), STRICT);
+    }
+
+    /**
      * Tests {@link ReferencingUtilities#isEllipsoidalHeight(VerticalDatum)}.
      */
     @Test
diff --git a/core/sis-utility/src/main/java/org/apache/sis/math/Vector.java b/core/sis-utility/src/main/java/org/apache/sis/math/Vector.java
index 7331c30..652b3e7 100644
--- a/core/sis-utility/src/main/java/org/apache/sis/math/Vector.java
+++ b/core/sis-utility/src/main/java/org/apache/sis/math/Vector.java
@@ -691,7 +691,7 @@ search:     for (;;) {
      * the [0 … {@link #size()} - 1] range:
      *
      * <blockquote><code>{@linkplain Math#abs(double) abs}({@linkplain #doubleValue(int)
doubleValue}(<var>i</var>)
-     * - (doubleValue(0) + increment*<var>i</var>)) &lt;= tolerance</code></blockquote>
+     * - (doubleValue(0) + increment*<var>i</var>)) ≦ tolerance</code></blockquote>
      *
      * The tolerance threshold can be zero if exact matches are desired.
      * The return value (if non-null) is always a signed value,


Mime
View raw message