sis-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From desruisse...@apache.org
Subject [sis] branch master updated: Handle the case where the longitude value of the point to transform using a datum shift is on the other side of the ±180° meridian compared to the region where the grid is defined. Example: NADCON file for Alaska which is defined in the [−194° … −127.875°] range of longitude value, so the +170° longitude needs to be converted to −190° before it can be used.
Date Wed, 19 Feb 2020 13:22:55 GMT
This is an automated email from the ASF dual-hosted git repository.

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


The following commit(s) were added to refs/heads/master by this push:
     new a4a9730  Handle the case where the longitude value of the point to transform using
a datum shift is on the other side of the ±180° meridian compared to the region where the
grid is defined. Example: NADCON file for Alaska which is defined in the [−194° … −127.875°]
range of longitude value, so the +170° longitude needs to be converted to −190° before
it can be used.
a4a9730 is described below

commit a4a9730d363df5775d860ee4b52d04c4d885a297
Author: Martin Desruisseaux <martin.desruisseaux@geomatys.com>
AuthorDate: Wed Feb 19 14:17:30 2020 +0100

    Handle the case where the longitude value of the point to transform using a datum shift
is on the other side of the ±180° meridian compared to the region where the grid is defined.
Example: NADCON file for Alaska which is defined in the [−194° … −127.875°] range
of longitude value, so the +170° longitude needs to be converted to −190° before it can
be used.
    
    https://issues.apache.org/jira/browse/SIS-488
---
 .../referencing/provider/DatumShiftGridFile.java   |  64 +++++++++--
 .../referencing/provider/DatumShiftGridGroup.java  |   2 +-
 .../sis/referencing/datum/DatumShiftGrid.java      | 124 +++++++++++++++++----
 3 files changed, 153 insertions(+), 37 deletions(-)

diff --git a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/DatumShiftGridFile.java
b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/DatumShiftGridFile.java
index f77bba7..7ebb2d7 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/DatumShiftGridFile.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/DatumShiftGridFile.java
@@ -37,6 +37,7 @@ import org.opengis.referencing.operation.MathTransformFactory;
 import org.opengis.referencing.operation.NoninvertibleTransformException;
 import org.opengis.referencing.operation.TransformException;
 import org.apache.sis.measure.Units;
+import org.apache.sis.measure.Longitude;
 import org.apache.sis.math.DecimalFunctions;
 import org.apache.sis.util.collection.Cache;
 import org.apache.sis.util.collection.TreeTable;
@@ -47,6 +48,7 @@ import org.apache.sis.util.resources.Errors;
 import org.apache.sis.parameter.Parameters;
 import org.apache.sis.referencing.datum.DatumShiftGrid;
 import org.apache.sis.internal.referencing.j2d.AffineTransform2D;
+import org.apache.sis.referencing.operation.matrix.AffineTransforms2D;
 import org.apache.sis.referencing.operation.transform.MathTransforms;
 import org.apache.sis.referencing.operation.transform.InterpolatedTransform;
 
@@ -62,9 +64,9 @@ import org.apache.sis.referencing.operation.transform.InterpolatedTransform;
  * @author  Martin Desruisseaux (Geomatys)
  * @version 1.1
  *
- * @param <C>  dimension of the coordinate unit (usually {@link javax.measure.quantity.Angle}).
- * @param <T>  dimension of the translation unit (usually {@link javax.measure.quantity.Angle}
- *             or {@link javax.measure.quantity.Length}).
+ * @param <C>  dimension of the coordinate unit (usually {@link Angle}).
+ * @param <T>  dimension of the translation unit. Usually {@link Angle},
+ *             but can also be {@link javax.measure.quantity.Length}.
  *
  * @see org.apache.sis.referencing.operation.transform.InterpolatedTransform
  *
@@ -115,6 +117,16 @@ abstract class DatumShiftGridFile<C extends Quantity<C>, T extends
Quantity<T>>
     protected final int nx;
 
     /**
+     * Number of cells that the grid would have if it was spanning 360° of longitude, or
0 if no wraparound
+     * should be applied. Current implementation rounds to nearest integer on the assumption
that we expect
+     * an integer number of cells in 360°. This value is used for longitude values that
are on the other side
+     * of the ±180° meridian compared to the region where the grid is defined.
+     *
+     * @see #replaceOutsideGridCoordinate(int, double)
+     */
+    private final double cycle;
+
+    /**
      * The best translation accuracy that we can expect from this file.
      * The unit of measurement depends on {@link #isCellValueRatio()}.
      *
@@ -150,8 +162,8 @@ abstract class DatumShiftGridFile<C extends Quantity<C>, T extends
Quantity<T>>
      * @param coordinateUnit    the unit of measurement of input values, before conversion
to grid indices by {@code coordinateToGrid}.
      * @param translationUnit   the unit of measurement of output values.
      * @param isCellValueRatio  {@code true} if results of {@link #interpolateInCell interpolateInCell(…)}
are divided by grid cell size.
-     * @param x0                longitude in degrees of the center of the cell at grid index
(0,0).
-     * @param y0                latitude in degrees of the center of the cell at grid index
(0,0).
+     * @param x0                longitude in degrees of the center of the cell at grid index
(0,0), positive east.
+     * @param y0                latitude in degrees of the center of the cell at grid index
(0,0), positive north.
      * @param Δx                increment in <var>x</var> value between cells
at index <var>gridX</var> and <var>gridX</var> + 1.
      * @param Δy                increment in <var>y</var> value between cells
at index <var>gridY</var> and <var>gridY</var> + 1.
      * @param nx                number of cells along the <var>x</var> axis in
the grid.
@@ -173,6 +185,16 @@ abstract class DatumShiftGridFile<C extends Quantity<C>, T extends
Quantity<T>>
         this.descriptor = descriptor;
         this.files      = files;
         this.nx         = nx;
+        if (Units.isAngular(coordinateUnit)) {
+            cycle = Math.rint((Longitude.MAX_VALUE - Longitude.MIN_VALUE) / Math.abs(Δx));
+        } else {
+            cycle = 0;
+            /*
+             * Note: non-angular source coordinates are currently never used in this package.
+             * If it continue to be like that in the future, we should remove the check for
+             * Units.isAngular(…) and replace the C parameterized type by Angle directly.
+             */
+        }
     }
 
     /**
@@ -188,26 +210,31 @@ abstract class DatumShiftGridFile<C extends Quantity<C>, T
extends Quantity<T>>
         nx         = other.nx;
         accuracy   = other.accuracy;
         subgrids   = other.subgrids;
+        cycle      = other.cycle;
     }
 
     /**
      * Creates a new datum shift grid with the same configuration than the given grid,
      * except the size and transform which are set to the given values.
+     * This is used for creating a {@link DatumShiftGridGroup} containing many grids,
+     * using one grid as a template for setting parameter values.
      * The {@link #accuracy} is initialized to zero and should be updated by the caller.
      *
-     * @param other             the other datum shift grid from which to copy parameters.
-     * @param coordinateToGrid  conversion from the "real world" coordinates to grid indices
including fractional parts.
-     * @param nx                number of cells along the <var>x</var> axis in
the grid.
-     * @param ny                number of cells along the <var>y</var> axis in
the grid.
+     * @param  other      the other datum shift grid from which to copy parameters.
+     * @param  gridToCRS  conversion from grid indices to "real world" coordinates.
+     * @param  nx         number of cells along the <var>x</var> axis in the
grid.
+     * @param  ny         number of cells along the <var>y</var> axis in the
grid.
      */
-    DatumShiftGridFile(final DatumShiftGridFile<C,T> other,
-                       final AffineTransform2D coordinateToGrid, final int nx, final int
ny)
+    DatumShiftGridFile(final DatumShiftGridFile<C,T> other, final AffineTransform2D
gridToCRS, final int nx, final int ny)
+            throws NoninvertibleTransformException
     {
-        super(other.getCoordinateUnit(), coordinateToGrid, new int[] {nx, ny},
+        super(other.getCoordinateUnit(), gridToCRS.inverse(), new int[] {nx, ny},
               other.isCellValueRatio(), other.getTranslationUnit());
         descriptor = other.descriptor;
         files      = other.files;
         this.nx    = nx;
+        cycle      = (other.cycle == 0) ? 0 :
+                Math.rint((Longitude.MAX_VALUE - Longitude.MIN_VALUE) / AffineTransforms2D.getScaleX0(gridToCRS));
         // Accuracy to be set by caller. Initial value needs to be zero.
     }
 
@@ -363,6 +390,19 @@ abstract class DatumShiftGridFile<C extends Quantity<C>, T extends
Quantity<T>>
     }
 
     /**
+     * Invoked when a {@code gridX} or {@code gridY} coordinate is outside the range of valid
grid coordinates.
+     * If the coordinate outside the range is a longitude value and if we handle those values
as cyclic, brings
+     * that coordinate inside the range.
+     */
+    @Override
+    protected double replaceOutsideGridCoordinate(final int dimension, final double gridCoordinate)
{
+        if (dimension == 0 && cycle != 0) {
+            return Math.IEEEremainder(gridCoordinate, cycle);
+        }
+        return super.replaceOutsideGridCoordinate(dimension, gridCoordinate);
+    }
+
+    /**
      * Returns the descriptor specified at construction time.
      *
      * @return a description of the values in this grid.
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/DatumShiftGridGroup.java
b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/DatumShiftGridGroup.java
index df2cefa..269c9c1 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/DatumShiftGridGroup.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/DatumShiftGridGroup.java
@@ -128,7 +128,7 @@ final class DatumShiftGridGroup<C extends Quantity<C>, T extends
Quantity<T>> ex
                                 final Dimension gridSize)
             throws IOException, NoninvertibleTransformException
     {
-        super(grids.get(tiles[0]), gridToCRS.inverse(), gridSize.width, gridSize.height);
+        super(grids.get(tiles[0]), gridToCRS, gridSize.width, gridSize.height);
         final int n = grids.size();
         regions  = new Region[n];
         subgrids = new DatumShiftGridFile[n];
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/datum/DatumShiftGrid.java
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/datum/DatumShiftGrid.java
index 3313b97..158540a 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/datum/DatumShiftGrid.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/datum/DatumShiftGrid.java
@@ -88,9 +88,9 @@ import org.apache.sis.measure.Units;
  *
  *   <li><b>Localization grid of raster data</b><br>
  *   Some remote sensing raster data are provided with a <cite>localization grid</cite>
giving pixel coordinates
- *   (e.g. latitude and longitude). This can been seen as a change from {@linkplain DefaultImageDatum
image datum}
- *   to {@linkplain DefaultGeodeticDatum geodetic datum}. The coordinate transformation process
can sometime be
- *   performed by a mathematical conversion (for example an affine transform) applied as
a
+ *   (e.g. latitude and longitude). This can been seen as a change from {@linkplain DefaultEngineeringDatum
+ *   image datum} to {@linkplain DefaultGeodeticDatum geodetic datum}. The coordinate transformation
process
+ *   can sometime be performed by a mathematical conversion (for example an affine transform)
applied as a
  *   {@linkplain org.apache.sis.referencing.operation.builder.LinearTransformBuilder first
approximation},
  *   followed by small corrections for the residual part.
  *   {@code DatumShiftGrid} can describe the small corrections part.
@@ -118,6 +118,14 @@ import org.apache.sis.measure.Units;
  * in more than two dimensions. See the above <cite>datum shift by geocentric translations</cite>
use case for
  * an example.
  *
+ * <h2>Longitude wraparound</h2>
+ * Some grids are defined over an area beyond the [−180° … +180°] range of longitudes.
+ * For example NADCON grid for Alaska is defined in a [−194° … −127.875°] range,
+ * in which case a longitude of 170° needs to be replaced by −190° before it can be processed
by the grid.
+ * The default {@code DatumShiftGrid} class does not apply longitude wraparound automatically
+ * (it does not even know which axis, if any, is longitude),
+ * but subclasses can add this support by overriding the {@link #replaceOutsideGridCoordinate(int,
double)} method.
+ *
  * <h2>Sub-grids</h2>
  * Some datum shift grid files provide a grid valid on a wide region, refined with denser
sub-grids in smaller regions.
  * For each point to transform, the {@link org.opengis.referencing.operation.MathTransform}
should search and use the
@@ -277,7 +285,9 @@ public abstract class DatumShiftGrid<C extends Quantity<C>, T
extends Quantity<T
     /**
      * Returns the domain of validity of input coordinates that can be specified to the
      * {@link #interpolateAt interpolateAt(…)} method. Coordinates outside that domain
-     * will still be accepted, but results will be extrapolations possibly far from reality.
+     * will still be accepted, but results may be extrapolations far from reality.
+     * This method does not take in account longitude wraparound
+     * (i.e. the returned envelope may cross the ±180° meridian).
      *
      * <p>The envelope coordinates are computed at cell centers; the envelope does
not contain
      * the margin of 0.5 cell between cell center and cell border at the edges of the envelope.
@@ -493,24 +503,32 @@ public abstract class DatumShiftGrid<C extends Quantity<C>,
T extends Quantity<T
      * @see #isCellInGrid(double, double)
      */
     public void interpolateInCell(double gridX, double gridY, final double[] vector) {
-        boolean skipX = false;
-        boolean skipY = false;                          // Whether to skip derivative calculation
for X or Y.
-        if (gridX < 0) {gridX = 0; skipX = true;}
-        if (gridY < 0) {gridY = 0; skipY = true;}
-        int ix = (int) gridX;  gridX -= ix;
-        int iy = (int) gridY;  gridY -= iy;
-        int n;
-        if (ix > (n = gridSize[0] - 2)) {
-            skipX |= (ix != n+1 || gridX != 0);         // Keep value 'false' if gridX ==
gridSize[0] - 1.
-            ix     = n;
-            gridX  = 1;
+        final int xmax = gridSize[0] - 2;
+        final int ymax = gridSize[1] - 2;
+        int ix = (int) gridX;                               // Really want rounding toward
zero (not floor).
+        if (ix < 0 || ix > xmax) {
+            gridX = replaceOutsideGridCoordinate(0, gridX);
+            ix = Math.max(0, Math.min(xmax, (int) gridX));
         }
-        if (iy > (n = gridSize[1] - 2)) {
-            skipY |= (iy != n+1 || gridY != 0);         // Keep value 'false' if gridY ==
gridSize[1] - 1.
-            iy     = n;
-            gridY  = 1;
+        int iy = (int) gridY;
+        if (iy < 0 || iy > ymax) {
+            gridY = replaceOutsideGridCoordinate(1, gridY);
+            iy = Math.max(0, Math.min(ymax, (int) gridY));
         }
-        n = getTranslationDimensions();
+        gridX -= ix;                                        // If was negative, will continue
to be negative.
+        gridY -= iy;
+        /*
+         * The `skipX` and `skipY` flags tell us whether a coordinate is outside the grid,
+         * in which case we will need to skip derivative calculation for component outside.
+         * More specifically, the Jacobian in dimensions outside the grid must be identity.
+         * Note that we want the `false` value if gridX == (gridSize[0] - 1) == (xmax + 1)
+         * if which case we have gridX = 1 (in all other cases, gridX < 1). Same for y.
+         */
+        boolean skipX = (gridX < 0); if (skipX) gridX = 0;
+        boolean skipY = (gridY < 0); if (skipY) gridY = 0;
+        if (gridX > 1)  {gridX = 1; skipX = true;}
+        if (gridY > 1)  {gridY = 1; skipY = true;}
+        final int n = getTranslationDimensions();
         boolean derivative = (vector.length >= n + INTERPOLATED_DIMENSIONS * INTERPOLATED_DIMENSIONS);
         for (int dim = 0; dim < n; dim++) {
             double dx, dy;
@@ -587,8 +605,18 @@ public abstract class DatumShiftGrid<C extends Quantity<C>,
T extends Quantity<T
      * @see #interpolateInCell(double, double, double[])
      */
     public Matrix derivativeInCell(double gridX, double gridY) {
-        final int ix = Math.max(0, Math.min(gridSize[0] - 2, (int) gridX));
-        final int iy = Math.max(0, Math.min(gridSize[1] - 2, (int) gridY));
+        final int xmax = gridSize[0] - 2;
+        final int ymax = gridSize[1] - 2;
+        int ix = (int) gridX;
+        if (ix < 0 || ix > xmax) {
+            gridX = replaceOutsideGridCoordinate(0, gridX);
+            ix = Math.max(0, Math.min(xmax, (int) gridX));
+        }
+        int iy = (int) gridY;
+        if (iy < 0 || iy > ymax) {
+            gridY = replaceOutsideGridCoordinate(1, gridY);
+            iy = Math.max(0, Math.min(ymax, (int) gridY));
+        }
         gridX -= ix;
         gridY -= iy;
         final boolean skipX = (gridX < 0 || gridX > 1);
@@ -716,8 +744,56 @@ public abstract class DatumShiftGrid<C extends Quantity<C>,
T extends Quantity<T
      *
      * @since 1.0
      */
-    public boolean isCellInGrid(final double gridX, final double gridY) {
-        return gridX >= 0 && gridY >= 0 && gridX <= gridSize[0]
- 1 && gridY <= gridSize[1] - 1;
+    public boolean isCellInGrid(double gridX, double gridY) {
+        double max = gridSize[1] - 1;
+        if (!(gridY >= 0 && gridY <= max)) {
+            gridY = replaceOutsideGridCoordinate(1, gridY);
+            if (!(gridY >= 0 && gridY <= max)) return false;
+        }
+        max = gridSize[0] - 1;
+        if (!(gridX >= 0 && gridX <= max)) {
+            gridX = replaceOutsideGridCoordinate(0, gridX);
+            if (!(gridX >= 0 && gridX <= max)) return false;
+        }
+        return true;
+    }
+
+    /**
+     * Invoked when a {@code gridX} or {@code gridY} coordinate is outside the range of valid
grid coordinates.
+     * This method can replace the invalid coordinate by a valid one. The main purpose is
to handle datum shift
+     * grids crossing the anti-meridian. For example NADCON grid for Alaska is defined in
a [−194° … −127.875°]
+     * longitude range, so a longitude of 170° needs to be converted to a longitude of −190°
before it can be
+     * processed by that grid.
+     *
+     * <p>The default implementation returns {@code gridCoordinate} unchanged. Subclasses
need to override this
+     * method if they want to handle longitude cycles. Note that the coordinate value is
a grid index, not a
+     * longitude value. So the cycle to add or remove is the number of cells that the grid
would have if it was
+     * spanning 360° of longitude.</p>
+     *
+     * <div class="note"><b>Example:</b>
+     * this method may be implemented as below:
+     *
+     * {@preformat java
+     *     &#64;Override
+     *     protected double replaceOutsideGridCoordinate(int dimension, double gridCoordinate)
{
+     *         if (dimension == 0) {
+     *             return Math.IEEEremainder(gridCoordinate, cycle);
+     *         } else {
+     *             return super.replaceOutsideGridCoordinate(dimension, gridCoordinate);
+     *         }
+     *     }
+     * }</div>
+     *
+     * @param  dimension       0 if the given coordinate is {@code gridX}, or 1 if it is
{@code gridY}.
+     * @param  gridCoordinate  the grid coordinate which is outside the range of valid values.
+     * @return the grid coordinate to use. May still be outside the range of valid grid values.
+     *
+     * @see Math#IEEEremainder(double, double)
+     *
+     * @since 1.1
+     */
+    protected double replaceOutsideGridCoordinate(final int dimension, final double gridCoordinate)
{
+        return gridCoordinate;
     }
 
     /**


Mime
View raw message