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));
}
}
|