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: First attempt to make localization grid robust to cases where the grid crosses the antimeridian.
Date Fri, 08 Feb 2019 12:23:16 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 a44d780  First attempt to make localization grid robust to cases where the grid crosses
the antimeridian.
a44d780 is described below

commit a44d7806ff351fe9d37152204de2526714a5672d
Author: Martin Desruisseaux <martin.desruisseaux@geomatys.com>
AuthorDate: Fri Feb 8 13:22:36 2019 +0100

    First attempt to make localization grid robust to cases where the grid crosses the antimeridian.
---
 .../operation/builder/LinearTransformBuilder.java  | 64 ++++++++++++++++++++--
 .../operation/builder/LocalizationGridBuilder.java | 47 ++++++++++++++++
 .../java/org/apache/sis/internal/netcdf/Axis.java  | 59 ++++++++++++++++++--
 3 files changed, 160 insertions(+), 10 deletions(-)

diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/LinearTransformBuilder.java
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/LinearTransformBuilder.java
index 0a69897..4138d69 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/LinearTransformBuilder.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/LinearTransformBuilder.java
@@ -908,11 +908,14 @@ search:         for (int j=domain(); --j >= 0;) {
         for (int i=0; i<tgtDim; i++) {
             final Vector c = coordinates[i];
             ArgumentChecks.ensureNonNullElement("coordinates", i, c);
-            final int size = c.size();
-            if (size != gridLength) {
-                throw new IllegalArgumentException(Errors.format(Errors.Keys.UnexpectedArrayLength_2,
gridLength, size));
+            int size = c.size();
+            if (size == gridLength) {
+                size = (result[i] = c.doubleValues()).length;
+                if (size == gridLength) {                       // Paranoiac check in case
user overwrite Vector.size().
+                    continue;
+                }
             }
-            result[i] = c.doubleValues();
+            throw new IllegalArgumentException(Errors.format(Errors.Keys.UnexpectedArrayLength_2,
gridLength, size));
         }
         targets     = result;
         transform   = null;
@@ -943,6 +946,59 @@ search:         for (int j=domain(); --j >= 0;) {
     }
 
     /**
+     * Tries to remove discontinuities in coordinates values caused by anti-meridian crossing.
This is the implementation of
+     * {@link LocalizationGridBuilder#resolveWraparoundAxis(int, int, double)} public method.
See that method for javadoc.
+     *
+     * @param  dimension  the dimension to process, from 0 inclusive to {@link #getTargetDimensions()}
exclusive.
+     *                    This is 0 for longitude dimension in a (<var>longitudes</var>,
<var>latitudes</var>) grid.
+     * @param  direction  the direction to walk through: 0 for columns or 1 for rows (higher
dimensions are also possible).
+     *                    Value can be from 0 inclusive to {@link #getSourceDimensions()}
exclusive.
+     *                    The recommended direction is the direction of most stable values,
typically 1 (rows) for longitudes.
+     * @param  period     that wraparound range (typically 360° for longitudes).
+     */
+    final void resolveWraparoundAxis(final int dimension, final int direction, final double
period) {
+        final double[] coordinates = targets[dimension];
+        int stride = 1;
+        for (int i=0; i<direction; i++) {
+            stride *= gridSize[i];                              // Index offset for moving
to next value in the specified direction.
+        }
+        final int page = stride * gridSize[direction];          // Index offset for moving
to next row or whatever is the next dimension.
+        final double threshold = period / 2;
+        double previous = coordinates[0];
+        for (int x=0; x<stride; x++) {                          // For iterating over
dimensions lower than 'dimension'.
+            for (int y=0; y<gridLength; y += page) {            // For iterating over
dimensions greater than 'dimension'.
+                final int stop = y + page;
+                for (int i = x+y; i<stop; i += stride) {
+                    double value = coordinates[i];
+                    double delta = value - previous;
+                    if (Math.abs(delta) > threshold) {
+                        delta = Math.rint(delta / period) * period;
+                        value -= delta;
+                        coordinates[i] = value;
+                    }
+                    previous = value;
+                }
+                /*
+                 * For the next scan, use as a reference the first value of this scan. If
our scan direction is 0
+                 * (each value compared with the value in previous column), then the first
value of next row will
+                 * be compared with the first value of this row. This is illustrated by index
-1 below:
+                 *
+                 *    ┌───┬───┬───┬───┬───┬───┐
+                 *    │-1 │   │   │   │   │   │      coordinates[x]
+                 *    ├───┼───┼───┼───┼───┼───┤
+                 *    │ 0 │ 1 │ 2 │ 3 │ 4 │ 5 │      next row to be scanned
+                 *    └───┴───┴───┴───┴───┴───┘
+                 *
+                 * Since the direction given in argument is the direction of most stable
values, the perpendicular
+                 * direction used for coordinates[x] may have more variation. We assume that
those variations are
+                 * still small enough for taking that nearby value as a reference.
+                 */
+                previous = coordinates[x];
+            }
+        }
+    }
+
+    /**
      * Returns the vector of source coordinates.
      * It is caller responsibility to ensure that this builder is not backed by a grid.
      */
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/LocalizationGridBuilder.java
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/LocalizationGridBuilder.java
index b43a0c2..b3fb15a 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/LocalizationGridBuilder.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/LocalizationGridBuilder.java
@@ -433,6 +433,53 @@ public class LocalizationGridBuilder extends TransformBuilder {
     }
 
     /**
+     * Tries to remove discontinuities in coordinates values caused by anti-meridian crossing.
+     * This method can be invoked when the localization grid may cross the anti-meridian,
+     * where longitude values may suddenly jump from +180° to -180° or conversely.
+     * This method walks through the coordinate values of the given dimension (typically
the longitudes dimension)
+     * in the given direction (grid rows or grid columns).
+     * If a difference greater than {@code period/2} (typically 180°) is found between two
consecutive values,
+     * then a multiple of {@code period} (typically 360°) is added or subtracted in order
to make a value as close
+     * as possible from its previous value.
+     *
+     * <p>This method needs a direction to be specified:</p>
+     * <ul>
+     *   <li>Direction 0 means that each value is compared with the value in the previous
column,
+     *       except the value in the first column which is compared to the value in previous
row.</li>
+     *   <li>Direction 1 means that each value is compared with the value in the previous
row,
+     *       except the value in the first row which is compared to the value in previous
column.</li>
+     * </ul>
+     * The recommended value is the direction of most stable values. Typically, longitude
values increase with column indices
+     * and are almost constant when increasing row indices. In such case, the recommended
direction is 1 for comparing each
+     * value with the value in previous row, since that value should be closer than the value
in previous column.
+     *
+     * <div class="note"><b>Example:</b>
+     * for a grid of (<var>longitude</var>, <var>latitude</var>)
values in decimal degrees where longitude values
+     * vary (increase or decrease) with increasing column indices and latitude values vary
(increase or decrease)
+     * with increasing row indices, the the following method should be invoked for protecting
the grid against
+     * discontinuities on anti-meridian:
+     *
+     * {@preformat java
+     *     grid.resolveWraparoundAxis(0, 1, 360);
+     * }
+     * </div>
+     *
+     * @param  dimension  the dimension to process.
+     *                    This is 0 for longitude dimension in a (<var>longitudes</var>,
<var>latitudes</var>) grid.
+     * @param  direction  the direction to walk through: 0 for columns or 1 for rows.
+     *                    The recommended direction is the direction of most stable values,
typically 1 (rows) for longitudes.
+     * @param  period     that wraparound range (typically 360° for longitudes).
+     *
+     * @since 1.0
+     */
+    public void resolveWraparoundAxis(final int dimension, final int direction, final double
period) {
+        ArgumentChecks.ensureBetween("dimension", 0, linear.getTargetDimensions() - 1, dimension);
+        ArgumentChecks.ensureBetween("direction", 0, linear.getSourceDimensions() - 1, direction);
+        ArgumentChecks.ensureStrictlyPositive("period", period);
+        linear.resolveWraparoundAxis(dimension, direction, period);
+    }
+
+    /**
      * Creates a transform from the source points to the target points.
      * This method assumes that source points are precise and all uncertainty is in the target
points.
      * If this transform is close enough to an affine transform, then an instance of {@link
LinearTransform} is returned.
diff --git a/storage/sis-netcdf/src/main/java/org/apache/sis/internal/netcdf/Axis.java b/storage/sis-netcdf/src/main/java/org/apache/sis/internal/netcdf/Axis.java
index d0ebae9..3398198 100644
--- a/storage/sis-netcdf/src/main/java/org/apache/sis/internal/netcdf/Axis.java
+++ b/storage/sis-netcdf/src/main/java/org/apache/sis/internal/netcdf/Axis.java
@@ -368,6 +368,25 @@ public final class Axis extends NamedElement {
     }
 
     /**
+     * Returns the range of this wraparound axis, or {@link Double#NaN} if this axis is not
a wraparound axis.
+     */
+    @SuppressWarnings("fallthrough")
+    private double wraparoundRange() {
+        if (isWraparound()) {
+            double period = 360;
+            final Unit<?> unit = getUnit();
+            if (unit != null) try {
+                period = unit.getConverterToAny(Units.DEGREE).convert(period);
+            } catch (IncommensurableException e) {
+                warning(e, Errors.Keys.InconsistentUnitsForCS_1, unit);
+                return Double.NaN;
+            }
+            return period;
+        }
+        return Double.NaN;
+    }
+
+    /**
      * Returns {@code true} if coordinates in this axis seem to map cell corner instead than
cell center.
      * A {@code false} value does not necessarily means that the axis maps cell center; it
can be unknown.
      * This method assumes a geographic CRS.
@@ -399,7 +418,7 @@ public final class Axis extends NamedElement {
                 }
                 return uc.convert(data.doubleValue(0)) == min;
             } catch (IncommensurableException e) {
-                coordinates.error(Grid.class, "getGridGeometry", e, Errors.Keys.InconsistentUnitsForCS_1,
unit);
+                warning(e, Errors.Keys.InconsistentUnitsForCS_1, unit);
             }
         }
         return false;
@@ -598,17 +617,32 @@ main:   switch (getDimension()) {
                 final int ro = (xo <= yo) ? 0 : 1;
                 final int width  = getSize(ri ^ 1);     // Fastest varying is right-most
dimension.
                 final int height = getSize(ri    );     // Slowest varying if left-most dimension.
-                if (other.sourceSizes[ro ^ 1] != width ||
-                    other.sourceSizes[ro    ] != height)
+                if (other.sourceSizes[ro ^ 1] == width &&
+                    other.sourceSizes[ro    ] == height)
                 {
-                    coordinates.error(Grid.class, "getGridGeometry", null,
-                            Errors.Keys.MismatchedGridGeometry_2, getName(), other.getName());
-                } else {
                     final LocalizationGridBuilder grid = new LocalizationGridBuilder(width,
height);
                     final Vector vx =  this.read();
                     final Vector vy = other.read();
                     grid.setControlPoints(vx, vy);
+                    /*
+                     * At this point we finished to set values in the localization grid,
but did not computed the transform yet.
+                     * Before to use the grid for calculation, we need to repair discontinuities
sometime found with longitudes.
+                     * If the grid crosses the anti-meridian, some values may suddenly jump
from +180° to -180° or conversely.
+                     * Even when not crossing the anti-meridian, we still observe apparently
random 360° jumps in some files,
+                     * especially close to poles. The methods invoked below try to make the
longitude grid more continuous.
+                     * The "ri" or "ro" argument specifies which dimension varies slowest,
i.e. which dimension have values
+                     * that do not change much when increasing longitudes. This is usually
1 (the rows).
+                     */
+                    double period;
+                    if (!Double.isNaN(period = wraparoundRange())) {
+                        grid.resolveWraparoundAxis(0, ri, period);
+                    }
+                    if (!Double.isNaN(period = other.wraparoundRange())) {
+                        grid.resolveWraparoundAxis(1, ro, period);
+                    }
                     return grid;
+                } else {
+                    warning(null, Errors.Keys.MismatchedGridGeometry_2, getName(), other.getName());
                 }
             }
         }
@@ -616,6 +650,19 @@ main:   switch (getDimension()) {
     }
 
     /**
+     * Reports a non-fatal error that occurred while constructing the grid geometry.
+     * This method pretends that the error come from {@link Grid#getGridGeometry(Decoder)},
+     * which is the caller and is a bit closer to a public API.
+     *
+     * @param  exception  the exception that occurred, or {@code null} if none.
+     * @param  key        one or {@link Errors.Keys} constants.
+     * @param  arguments  values to be formatted in the {@link java.text.MessageFormat} pattern.
+     */
+    private void warning(final Exception exception, final short key, final Object... arguments)
{
+        coordinates.error(Grid.class, "getGridGeometry", exception, key, arguments);
+    }
+
+    /**
      * Returns the coordinates in the localization grid, excluding some trailing NaN values
if any.
      * This method typically returns a cached vector if the coordinates have already been
read.
      *


Mime
View raw message