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 54fd75a Refactor WraparoundAdjustment in a way that makes easier to share code between the method working on Envelope and the method working on DirectPosition. 54fd75a is described below commit 54fd75a91ee549e9474e8431850262122afb3ac5 Author: Martin Desruisseaux AuthorDate: Sun Apr 7 21:05:56 2019 +0200 Refactor WraparoundAdjustment in a way that makes easier to share code between the method working on Envelope and the method working on DirectPosition. --- .../apache/sis/coverage/grid/GridDerivation.java | 69 +-- .../sis/coverage/grid/GridDerivationTest.java | 73 ++- .../internal/referencing/WraparoundAdjustment.java | 568 +++++++++++++-------- .../referencing/WraparoundAdjustmentTest.java | 15 +- 4 files changed, 430 insertions(+), 295 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 a0172e5..1c7ef9c 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 @@ -76,6 +76,7 @@ import org.opengis.coverage.PointOutsideCoverageException; * For dimensionality reduction, see {@link GridGeometry#reduce(int...)}.

* * @author Martin Desruisseaux (Geomatys) + * @author Alexis Manin (Geomatys) * @version 1.0 * * @see GridGeometry#derive() @@ -490,7 +491,7 @@ public class GridDerivation { cornerToCRS = MathTransforms.concatenate(cornerToCRS, baseToAOI.getMathTransform()); } /* - * If the envelope dimensions does not encompass all grid dimensions, the envelope is probably non-invertible. + * If the envelope dimensions do not encompass all grid dimensions, the transform is probably non-invertible. * We need to reduce the number of grid dimensions in the transform for having a one-to-one relationship. */ int dimension = cornerToCRS.getTargetDimensions(); @@ -504,9 +505,8 @@ public class GridDerivation { dimension = baseExtent.getDimension(); // Non-null since 'base.requireGridToCRS()' succeed. GeneralEnvelope indices = null; if (areaOfInterest != null) { - indices = new WraparoundAdjustment(areaOfInterest) - .shiftInto(base.envelope, (baseToAOI != null) ? baseToAOI.getMathTransform() : null) - .result(cornerToCRS.inverse()); + final MathTransform domainToAOI = (baseToAOI != null) ? baseToAOI.getMathTransform() : null; + indices = new WraparoundAdjustment(base.envelope, domainToAOI, cornerToCRS.inverse()).shift(areaOfInterest); clipExtent(indices); } if (indices == null || indices.getDimension() != dimension) { @@ -721,15 +721,17 @@ public class GridDerivation { * before to invoke this method. *

This method does not reduce the number of dimensions of the grid geometry. * For dimensionality reduction, see {@link GridGeometry#reduce(int...)}.

+ *

If the given point is known to be expressed in the same CRS than the grid geometry, + * then the {@linkplain DirectPosition#getCoordinateReferenceSystem() CRS of the point} + * can be left unspecified ({@code null}). It may give a slight performance improvement + * by avoiding the check for coordinate transformation.

* * * @param slicePoint the coordinates where to get a slice. If no coordinate reference system is attached to it, * we consider it's the same as base grid geometry. * @return {@code this} for method call chaining. - * @throws IncompleteGridGeometryException if the base grid geometry has no extent, no "grid to CRS" transform, or - * if its coordinate reference system is unknown while input point one is (in which case we're unable to find a - * transform between base CRS and input point CRS). Note that if you're sure that your point is expressed in base - * CRS, you can give a point without any CRS, to avoid this check. + * @throws IncompleteGridGeometryException if the base grid geometry has no extent, no "grid to CRS" transform, + * or no CRS (unless {@code slicePoint} has no CRS neither, in which case the CRS are assumed the same). * @throws IllegalGridGeometryException if an error occurred while converting the point coordinates to grid coordinates. * @throws PointOutsideCoverageException if the given point is outside the grid extent. */ @@ -741,40 +743,47 @@ public class GridDerivation { if (toBase != null) { gridToCRS = MathTransforms.concatenate(toBase, gridToCRS); } - - /* We'll try to find a link between base coordinate system and input point one. Note that we allow unknown - * CRS on slice point, in which case we consider it to be expressed in base geometry system. However, if the - * point CRS is known, but base geometry one is not, too much ambiguity resides, and we'll throw an error. + /* + * We will try to find a path between grid coordinate reference system (CRS) and given point CRS. Note that we + * allow unknown CRS on the slice point, in which case we consider it to be expressed in grid reference system. + * However, if the point CRS is specified while the base grid CRS is unknown, we are at risk of ambiguity, + * in which case we throw (indirectly) an IncompleteGridGeometryException. */ - final MathTransform baseToSlice; - final CoordinateReferenceSystem targetCRS = slicePoint.getCoordinateReferenceSystem(); - if (targetCRS != null) { - final CoordinateReferenceSystem sourceCRS = base.getCoordinateReferenceSystem(); - baseToSlice = CRS.findOperation(sourceCRS, targetCRS, null) - .getMathTransform(); - gridToCRS = MathTransforms.concatenate(gridToCRS, baseToSlice); + final CoordinateReferenceSystem sliceCRS = slicePoint.getCoordinateReferenceSystem(); + final MathTransform baseToPOI; + if (sliceCRS == null) { + baseToPOI = null; } else { - baseToSlice = null; + final CoordinateReferenceSystem gridCRS = base.getCoordinateReferenceSystem(); // May throw exception. + baseToPOI = CRS.findOperation(gridCRS, sliceCRS, null).getMathTransform(); + gridToCRS = MathTransforms.concatenate(gridToCRS, baseToPOI); } - + /* + * If the point dimensions do not encompass all grid dimensions, the transform is probably non-invertible. + * We need to reduce the number of grid dimensions in the transform for having a one-to-one relationship. + */ final int dimension = gridToCRS.getTargetDimensions(); ArgumentChecks.ensureDimensionMatches("slicePoint", dimension, slicePoint); gridToCRS = dropUnusedDimensions(gridToCRS, dimension); - - GeneralEnvelope sliceEnvelope = new GeneralEnvelope(slicePoint, slicePoint); - sliceEnvelope = new WraparoundAdjustment(sliceEnvelope) - .shiftInto(base.envelope, baseToSlice) - .result(gridToCRS.inverse()); - DirectPosition gridPoint = sliceEnvelope.getMedian(); + /* + * Take in account the case where the point could be inside the grid if we apply a ±360° longitude shift. + * This is the same adjustment than for `subgrid(Envelope)`, but applied on a DirectPosition. Calculation + * is done in units of cells of the GridGeometry to be created by GridDerivation. + */ + DirectPosition gridPoint = new WraparoundAdjustment(base.envelope, baseToPOI, gridToCRS.inverse()).shift(slicePoint); if (scaledExtent != null) { scaledExtent = scaledExtent.slice(gridPoint, modifiedDimensions); } - - // Rebuild the extent without scaling + /* + * Above `scaledExtent` is non-null only if a scale or subsampling has been applied before this `slice` + * method has been invoked. The `baseExtent` below contains same information, but without subsampling. + * The subsampling effect is removed by applying the "scaled to base" transform. Accuracy matter less + * here than for `scaledExtent` since the extent to be returned to user is the later. + */ if (toBase != null) { gridPoint = toBase.transform(gridPoint, gridPoint); } - baseExtent = baseExtent.slice(gridPoint, modifiedDimensions); // Non-null check by 'base.requireGridToCRS()'. + baseExtent = baseExtent.slice(gridPoint, modifiedDimensions); // Non-null check by 'base.requireGridToCRS()'. } catch (FactoryException e) { throw new IllegalGridGeometryException(Resources.format(Resources.Keys.CanNotMapToGridDimensions), e); } catch (TransformException e) { 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 9c57c9d..58a197a 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 @@ -16,18 +16,19 @@ */ package org.apache.sis.coverage.grid; -import org.apache.sis.geometry.Envelope2D; -import org.apache.sis.internal.referencing.j2d.AffineTransform2D; -import org.apache.sis.referencing.CommonCRS; -import org.apache.sis.referencing.crs.DefaultCompoundCRS; -import org.apache.sis.referencing.operation.matrix.Matrices; -import org.junit.Assert; +import java.util.Collections; +import java.util.stream.DoubleStream; +import java.util.stream.IntStream; import org.opengis.geometry.DirectPosition; import org.opengis.geometry.Envelope; import org.opengis.metadata.spatial.DimensionNameType; 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.internal.referencing.j2d.AffineTransform2D; +import org.apache.sis.referencing.crs.DefaultCompoundCRS; +import org.apache.sis.referencing.operation.matrix.Matrices; import org.apache.sis.referencing.operation.matrix.Matrix3; import org.apache.sis.referencing.operation.matrix.Matrix4; import org.apache.sis.referencing.operation.transform.MathTransforms; @@ -40,10 +41,6 @@ import org.apache.sis.test.DependsOn; import org.apache.sis.test.TestCase; import org.junit.Test; -import java.util.Collections; -import java.util.stream.DoubleStream; -import java.util.stream.IntStream; - import static org.apache.sis.test.ReferencingAssert.*; import static org.apache.sis.coverage.grid.GridGeometryTest.assertExtentEquals; @@ -52,6 +49,7 @@ import static org.apache.sis.coverage.grid.GridGeometryTest.assertExtentEquals; * Tests the {@link GridDerivation} implementation. * * @author Martin Desruisseaux (Geomatys) + * @author Alexis Manin (Geomatys) * @version 1.0 * @since 1.0 * @module @@ -264,16 +262,15 @@ public final strictfp class GridDerivationTest extends TestCase { } /** - * Check that wrap-around is well applied when using {@link GridDerivation#slice(DirectPosition)}. + * Checks that wraparound is well applied when using {@link GridDerivation#slice(DirectPosition)}. */ @Test - public void testSliceWrapAround() { + public void testSliceWithWrapAround() { final GridGeometry base = new GridGeometry( PixelInCell.CELL_CORNER, new AffineTransform2D(-0.02, 0, 0, 0.1, 55, 172), - new Envelope2D(CommonCRS.WGS84.geographic(), 42, 172, 13, 51), - GridRoundingMode.NEAREST - ); + new Envelope2D(HardCodedCRS.WGS84_φλ, 42, 172, 13, 51), + GridRoundingMode.NEAREST); final GridGeometry expectedResult = base.derive() .slice(new DirectPosition2D(51, 187)) @@ -283,28 +280,27 @@ public final strictfp class GridDerivationTest extends TestCase { .slice(new DirectPosition2D(51, -173)) .build(); - Assert.assertEquals("Slice with wrap-around", expectedResult, fromWrapAround); - assertBetween( - "Wrapped Y coordinate", - base.envelope.getMinimum(1), - base.envelope.getMaximum(1), - fromWrapAround.envelope.getMedian(1) - ); + assertEquals("Slice with wrap-around", expectedResult, fromWrapAround); + assertBetween("Wrapped Y coordinate", + base.envelope.getMinimum(1), + base.envelope.getMaximum(1), + fromWrapAround.envelope.getMedian(1)); } /** - * Ensure that slicing on a corner point does not fail, but gives back a geometry centered on a pixel corner. - * @throws TransformException If we cannot build our test point. + * Ensures that slicing on a corner point does not fail, but gives back a grid geometry centered on a pixel corner. + * + * @throws TransformException if we cannot build our test point. */ @Test - public void testSliceCorner() throws TransformException { + public void testSliceOnCorner() throws TransformException { GridGeometry base = grid(-132, 327, 986, 597, 2, 3); - - // First of all, we'll try by focusing on the last pixel. + /* + * First of all, we'll try by focusing on the last pixel. + */ final DirectPosition2D gridUpperCorner = new DirectPosition2D( base.extent.getHigh(0), - base.extent.getLow(1) - ); + base.extent.getLow(1)); final DirectPosition geoUpperCorner = base.getGridToCRS(PixelInCell.CELL_CENTER) .transform(gridUpperCorner, null); @@ -315,27 +311,26 @@ public final strictfp class GridDerivationTest extends TestCase { long[] expectedGridPoint = {(long) gridUpperCorner.x, (long) gridUpperCorner.y}; assertExtentEquals(expectedGridPoint, expectedGridPoint, slice.extent); - - /* We will now try to focus on a point near the envelope edge. Note that slicing ensures to return a valid grid - * for any point INSIDE the envelope, it's non-deterministic about points perfectly aligned on the edge. So, - * here we will test a point very near to the envelope edge, but still into it. + /* + * We will now try to focus on a point near the envelope edge. Note that slicing ensures to return a valid grid + * for any point INSIDE the envelope, it's non-determinist about points perfectly aligned on the edge. + * So, here we will test a point very near to the envelope edge, but still into it. */ final GeneralEnvelope grid3d = new GeneralEnvelope(3); grid3d.setEnvelope(0, 0, 0, 1920, 1080, 4); final DefaultCompoundCRS crs3d = new DefaultCompoundCRS( Collections.singletonMap("name", "geo3d"), - CommonCRS.defaultGeographic(), - CommonCRS.Temporal.JULIAN.crs() - ); + HardCodedCRS.WGS84, + HardCodedCRS.TIME); + final GeneralEnvelope geo3d = new GeneralEnvelope(crs3d); geo3d.setEnvelope(-180, -90, 1865.128, 180, 90, 1865.256); base = new GridGeometry( PixelInCell.CELL_CORNER, MathTransforms.linear(Matrices.createTransform(grid3d, geo3d)), geo3d, - GridRoundingMode.NEAREST - ); + GridRoundingMode.NEAREST); final GeneralDirectPosition geo3dUpperCorner = new GeneralDirectPosition(geo3d.getUpperCorner()); IntStream.range(0, geo3dUpperCorner.getDimension()) @@ -348,7 +343,7 @@ public final strictfp class GridDerivationTest extends TestCase { // Build expected grid point focused after slicing. We expect it to be upper corner. expectedGridPoint = DoubleStream.of(grid3d.getUpperCorner().getCoordinate()) .mapToLong(value -> (long) value) - .map(exclusiveValue -> exclusiveValue -1)// Exclusive to inclusive + .map(exclusiveValue -> exclusiveValue - 1) // Exclusive to inclusive. .toArray(); assertExtentEquals(expectedGridPoint, expectedGridPoint, slice.extent); diff --git a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/WraparoundAdjustment.java b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/WraparoundAdjustment.java index fb66213..a18fa73 100644 --- a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/WraparoundAdjustment.java +++ b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/WraparoundAdjustment.java @@ -35,11 +35,13 @@ import org.apache.sis.measure.Units; import org.apache.sis.math.MathFunctions; import org.apache.sis.geometry.Envelopes; import org.apache.sis.geometry.GeneralEnvelope; +import org.apache.sis.geometry.GeneralDirectPosition; /** * Adjustments applied on an envelope for handling wraparound axes. The adjustments consist in shifting - * some axes by an integer amount of periods, typically (not necessarily) 360° of longitude. + * coordinate values on some axes by an integer amount of periods (typically 360° of longitude) in order + * to move an envelope or a position inside a given domain of validity. * * @author Martin Desruisseaux (Geomatys) * @version 1.0 @@ -48,25 +50,70 @@ import org.apache.sis.geometry.GeneralEnvelope; */ public final class WraparoundAdjustment { /** - * The envelope to potentially shift in order to fit in the domain of validity. If a shift is needed, then - * this envelope will be replaced by a new envelope; the user-specified envelope will not be modified. + * The region inside which a given Area Of Interest (AOI) or Point Of Interest (POI) should be located. + * This envelope may be initially in a projected CRS and converted later to geographic CRS in order to + * allow identification of wraparound axes. */ - private Envelope areaOfInterest; + private Envelope domainOfValidity; /** - * If {@link #areaOfInterest} has been converted to a geographic CRS, the transformation back to its original CRS. - * Otherwise {@code null}. + * If the AOI or POI does not use the same CRS than {@link #domainOfValidity}, the transformation from + * {@code domainOfValidity} to the AOI / POI. Otherwise {@code null}. + * + *

Note: + * this class does not check by itself if a coordinate operation is needed; it must be supplied. We do that + * because {@code WraparoundAdjustment} is currently used in contexts where this transform is known anyway, + * so we avoid to compute it twice.

+ */ + private final MathTransform domainToAOI; + + /** + * A transform from the {@link #domainOfValidity} CRS to any user space at caller choice. + * The object returned by {@code shift} will be transformed by this transform after all computations + * have been finished. This is done in order to allow final transforms to be concatenated in a single step. + */ + private final MathTransform domainToAny; + + /** + * If {@code areaOfInterest} or {@code pointOfInterest} has been converted to a geographic CRS, + * the transformation back to its original CRS. Otherwise {@code null}. */ private MathTransform geographicToAOI; /** - * Creates a new instance for adjusting the given envelope. - * The given envelope will not be modified; a copy will be created if needed. + * The coordinate reference system of the Area Of Interest (AOI) or Point of Interest (POI). + * May be replaced by a geographic CRS if the AOI or POI originally used a projected CRS. + */ + private CoordinateReferenceSystem crs; + + /** + * Whether {@link #domainOfValidity} has been transformed to the geographic CRS that is the source + * of {@link #geographicToAOI}. This flag is used for ensuring that {@link #replaceCRS()} performs + * the inverse projection only once. + */ + private boolean isDomainTransformed; + + /** + * Whether the Area Of Interest (AOI) or Point Of Interest (POI) has been transformed in order + * to allow identification of wraparound axes. If {@code true}, then {@link #geographicToAOI} + * needs to be applied in order to restore the AOI or POI to its original projected CRS. + */ + private boolean isResultTransformed; + + /** + * Creates a new instance for adjusting an Area Of Interest (AOI) or Point Of Interest (POI) to the given + * domain of validity. The AOI or POI will be given later, but this method nevertheless requires in advance + * the transform from {@code domainOfValidity} to AOI or POI. * - * @param areaOfInterest the envelope to potentially shift toward the domain of validity. + * @param domainOfValidity the region where a given area or point of interest should be located. + * @param domainToAOI if the AOI or POI are going to use a different CRS than {@code domainOfValidity}, the + * transform from {@code domainOfValidity} to the AOI or POI CRS. Otherwise {@code null}. + * @param domainToAny a transform from the {@code domainOfValidity} CRS to any user space at caller choice. */ - public WraparoundAdjustment(final Envelope areaOfInterest) { - this.areaOfInterest = areaOfInterest; + public WraparoundAdjustment(final Envelope domainOfValidity, final MathTransform domainToAOI, final MathTransform domainToAny) { + this.domainOfValidity = domainOfValidity; + this.domainToAOI = domainToAOI; + this.domainToAny = domainToAny; } /** @@ -74,248 +121,333 @@ public final class WraparoundAdjustment { * 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 cs the coordinate system for which to get wraparound range. * @param dimension dimension of the axis to test. * @return the wraparound range, or {@link Double#NaN} if none. */ static double range(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; + 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; } /** - * Computes an envelope with coordinates equivalent to the {@code areaOfInterest} specified - * at construction time, but potentially shifted for intersecting the given domain of validity. - * The dimensions that may be shifted are the ones having an axis with wraparound meaning. - * The envelope may have been converted to a geographic CRS for performing this operation. - * - *

The coordinate reference system must be specified in the {@code areaOfInterest} - * specified at construction time, or (as a fallback) in the given {@code domainOfValidity}. - * If none of those envelopes have a CRS, then this method does nothing.

+ * Sets {@link #crs} to the given value if it is non-null, or to the {@link #domainOfValidity} CRS otherwise. + * If no non-null CRS is available, returns {@code false} for instructing caller to terminate immediately. * - *

This method does not intersect the area of interest with the domain of validity. - * It is up to the caller to compute that intersection after this method call, if desired.

- * - * @param domainOfValidity the domain of validity, or {@code null} if none. - * @param validToAOI if the envelopes do not use the same CRS, the transformation from {@code domainOfValidity} - * to {@code areaOfInterest}. Otherwise {@code null}. This method does not check by itself if - * a coordinate operation is needed; it must be supplied. - * @return this object, allowing to chain {@link #result(MathTransform)} operation. - * @throws TransformException if an envelope transformation was required but failed. - * - * @see GeneralEnvelope#simplify() + * @return whether a non-null CRS has been set. */ - public WraparoundAdjustment shiftInto(Envelope domainOfValidity, MathTransform validToAOI) throws TransformException { - CoordinateReferenceSystem crs = areaOfInterest.getCoordinateReferenceSystem(); + private boolean setIfNonNull(CoordinateReferenceSystem crs) { if (crs == null) { - crs = domainOfValidity.getCoordinateReferenceSystem(); // Assumed to apply to AOI too. + assert domainToAOI == null || domainToAOI.isIdentity(); + crs = domainOfValidity.getCoordinateReferenceSystem(); // Assumed to apply to AOI or POI too. if (crs == null) { - return this; + return false; } } - /* - * If the coordinate reference system is a projected CRS, it will not have any wraparound axis. - * We need to perform the verification in its base geographic CRS instead, and remember that we - * may need to transform the result later. - */ - final MathTransform projection; - final DirectPosition lowerCorner; - final DirectPosition upperCorner; - GeneralEnvelope shifted = null; // To be initialized to a copy of 'areaOfInterest' when first needed. + this.crs = crs; + return true; + } + + /** + * If the coordinate reference system is a projected CRS, replaces it by another CRS where wraparound axes can + * be identified. The wraparound axes are identifiable in base geographic CRS. If such replacement is applied, + * remember that we may need to transform the result later. + * + * @return whether the replacement has been done. If {@code true}, then {@link #geographicToAOI} is non-null. + */ + private boolean replaceCRS() { if (crs instanceof ProjectedCRS) { final ProjectedCRS p = (ProjectedCRS) crs; - crs = p.getBaseCRS(); - projection = p.getConversionFromBase().getMathTransform(); - shifted = Envelopes.transform(projection.inverse(), areaOfInterest); - lowerCorner = shifted.getLowerCorner(); - upperCorner = shifted.getUpperCorner(); - if (validToAOI == null) { - validToAOI = MathTransforms.identity(projection.getTargetDimensions()); - } + crs = p.getBaseCRS(); // Geographic, so a wraparound axis certainly exists. + geographicToAOI = p.getConversionFromBase().getMathTransform(); + return true; } else { - projection = null; - lowerCorner = areaOfInterest.getLowerCorner(); - upperCorner = areaOfInterest.getUpperCorner(); + return false; } - /* - * We will not read 'areaOfInterest' anymore after we got its two corner points. - * The following loop search for "wraparound" axis. - */ - final CoordinateSystem cs = crs.getCoordinateSystem(); - for (int i=cs.getDimension(); --i >= 0;) { - final double period = range(cs, i); - if (period > 0) { - /* - * Found an axis (typically the longitude axis) with wraparound range meaning. - * We are going to need the domain of validity in the same CRS than the AOI. - * Transform that envelope when first needed. - */ - if (validToAOI != null) { - if (projection != null) { - validToAOI = MathTransforms.concatenate(validToAOI, projection.inverse()); - } - if (!validToAOI.isIdentity()) { - domainOfValidity = Envelopes.transform(validToAOI, domainOfValidity); - } - validToAOI = null; - } - /* - * "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 domain of validity. Compute the distance from current - * lower/upper coordinate to the coordinate of validity domain (the sign tells us whether we - * are before or after). The cases can be: - * - * ┌─────────────┬────────────┬────────────────────────────┬───────────────────────────────┐ - * │lowerIsBefore│upperIsAfter│ Meaning │ Action │ - * ├─────────────┼────────────┼────────────────────────────┼───────────────────────────────┤ - * │ false │ false │ AOI is inside valid area │ Nothing to do │ - * │ true │ true │ AOI encompasses valid area │ Nothing to do │ - * │ true │ false │ AOI on left of valid area │ Add positive amount of period │ - * │ false │ true │ AOI on right of valid area │ Add negative amount of period │ - * └─────────────┴────────────┴────────────────────────────┴───────────────────────────────┘ - * - * 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 under the trajectory of a satellite. - */ - final double validStart = domainOfValidity.getMinimum(i); - final double validEnd = domainOfValidity.getMaximum(i); - final double lowerToValidStart = ((validStart - lower) / period) - lowerCycles; // In number of periods. - final double upperToValidEnd = ((validEnd - upper) / period) - upperCycles; - final boolean lowerIsBefore = (lowerToValidStart > 0); - final boolean upperIsAfter = (upperToValidEnd < 0); - if (lowerIsBefore != upperIsAfter) { - final double upperToValidStart = ((validStart - upper) / period) - upperCycles; - final double lowerToValidEnd = ((validEnd - lower) / period) - lowerCycles; - if (lowerIsBefore) { - /* - * Notation: ⎣x⎦=floor(x) and ⎡x⎤=ceil(x). - * - * We need to add an integer amount of 'period' to both sides in order to move the range - * inside the valid area. We need ⎣lowerToValidStart⎦ for reaching the point where: - * - * (validStart - period) < (new lower) ≦ validStart - * - * But we may add more because there will be no intersection without following condition: - * - * (new upper) ≧ validStart - * - * That second condition is met by ⎡upperToValidStart⎤. However adding more may cause the - * range to move the AOI completely on the right side of the domain of validity. We prevent - * that with a third condition: - * - * (new lower) < validEnd - */ - final double cycles = Math.min(Math.floor(lowerToValidEnd), - Math.max(Math.floor(lowerToValidStart), - Math.ceil (upperToValidStart))); - /* - * If after the shift we see that the following condition hold: - * - * (new lower) + period < validEnd - * - * Then we may have a situation like below: - * ┌────────────────────────────────────────────┐ - * │ Domain of validity │ - * └────────────────────────────────────────────┘ - * ┌────────────────────┐ ┌───── - * │ Area of interest │ │ AOI - * └────────────────────┘ └───── - * ↖……………………………………………………………period……………………………………………………………↗︎ - * - * The user may be requesting two extremums of the domain of validity. We can not express - * that with a single envelope. Instead, we will expand the Area Of Interest to encompass - * the full domain of validity. - */ - if (cycles + 1 < lowerToValidEnd) { - upperCycles += Math.ceil(upperToValidEnd); + } + + /** + * Transforms {@link #domainOfValidity} to the same CRS than the Area Of Interest (AOI) or Point Of Interest (POI). + * This method should be invoked only when the caller detected a wraparound axis. This method transforms the domain + * the first time it is invoked, and does nothing on all subsequent calls. + */ + private void transformDomainToAOI() throws TransformException { + if (!isDomainTransformed) { + isDomainTransformed = true; + MathTransform domainToGeographic = domainToAOI; + if (domainToGeographic == null) { + domainToGeographic = geographicToAOI; + } else if (geographicToAOI != null) { + domainToGeographic = MathTransforms.concatenate(domainToGeographic, geographicToAOI.inverse()); + } + if (domainToGeographic != null && !domainToGeographic.isIdentity()) { + domainOfValidity = Envelopes.transform(domainToGeographic, domainOfValidity); + } + } + } + + /** + * Returns the final transform to apply on the AOI or POI before to return it to the user. + */ + private MathTransform toFinal() throws TransformException { + MathTransform mt = domainToAny; + if (isResultTransformed && geographicToAOI != null) { + mt = MathTransforms.concatenate(geographicToAOI, mt); + } + return mt; + } + + /** + * Computes an envelope with coordinates equivalent to the given {@code areaOfInterest}, but + * potentially shifted for intersecting the domain of validity specified at construction time. + * The dimensions that may be shifted are the ones having an axis with wraparound meaning. + * In order to perform this operation, the envelope may be temporarily converted to a geographic CRS + * and converted back to its original CRS. + * + *

The coordinate reference system should be specified in the {@code areaOfInterest}, + * or (as a fallback) in the {@code domainOfValidity} specified at construction time.

+ * + *

This method does not intersect the area of interest with the domain of validity. + * It is up to the caller to compute that intersection after this method call, if desired.

+ * + * @param areaOfInterest the envelope to potentially shift toward domain of validity. + * If a shift is needed, then given envelope will be replaced by a new envelope; + * the given envelope will not be modified. + * @return envelope potentially expanded or shifted toward the domain of validity. + * @throws TransformException if a coordinate conversion failed. + * + * @see GeneralEnvelope#simplify() + */ + public GeneralEnvelope shift(Envelope areaOfInterest) throws TransformException { + if (setIfNonNull(areaOfInterest.getCoordinateReferenceSystem())) { + /* + * If the coordinate reference system is a projected CRS, it will not have any wraparound axis. + * We need to perform the verification in its base geographic CRS instead, and remember that we + * may need to transform the result later. + */ + final DirectPosition lowerCorner; + final DirectPosition upperCorner; + GeneralEnvelope shifted; // To be initialized to a copy of 'areaOfInterest' when first needed. + if (replaceCRS()) { + shifted = Envelopes.transform(geographicToAOI.inverse(), areaOfInterest); + lowerCorner = shifted.getLowerCorner(); + upperCorner = shifted.getUpperCorner(); + } else { + shifted = null; + lowerCorner = areaOfInterest.getLowerCorner(); + upperCorner = areaOfInterest.getUpperCorner(); + } + /* + * We will not read 'areaOfInterest' anymore after we got its two corner points (except for creating + * a copy if `shifted` is still null). The following loop searches for "wraparound" axes. + */ + final CoordinateSystem cs = crs.getCoordinateSystem(); + for (int i=cs.getDimension(); --i >= 0;) { + final double period = range(cs, i); + if (period > 0) { + /* + * Found an axis (typically the longitude axis) with wraparound range meaning. + * We are going to need the domain of validity in the same CRS than the AOI. + * Transform that envelope when first needed. + */ + transformDomainToAOI(); + /* + * "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; + upperCycles = -cycles; // Will add periods to 'upper'. } - lowerCycles += cycles; - } else { - /* - * Same reasoning than above with sign reverted and lower/upper variables interchanged. - * In this block, 'upperToValidEnd' and 'lowerToValidEnd' are negative, contrarily to - * above block where they were positive. - */ - final double cycles = Math.max(Math.ceil (upperToValidStart), - Math.min(Math.ceil (upperToValidEnd), - Math.floor(lowerToValidEnd))); - if (cycles - 1 > upperToValidStart) { - lowerCycles += Math.floor(lowerToValidStart); - } else { + } + /* + * The range may be before or after the domain of validity. Compute the distance from current + * lower/upper coordinate to the coordinate of validity domain (the sign tells us whether we + * are before or after). The cases can be: + * + * ┌─────────────┬────────────┬────────────────────────────┬───────────────────────────────┐ + * │lowerIsBefore│upperIsAfter│ Meaning │ Action │ + * ├─────────────┼────────────┼────────────────────────────┼───────────────────────────────┤ + * │ false │ false │ AOI is inside valid area │ Nothing to do │ + * │ true │ true │ AOI encompasses valid area │ Nothing to do │ + * │ true │ false │ AOI on left of valid area │ Add positive amount of period │ + * │ false │ true │ AOI on right of valid area │ Add negative amount of period │ + * └─────────────┴────────────┴────────────────────────────┴───────────────────────────────┘ + * + * 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 under the trajectory of a satellite. + */ + final double validStart = domainOfValidity.getMinimum(i); + final double validEnd = domainOfValidity.getMaximum(i); + final double lowerToValidStart = ((validStart - lower) / period) - lowerCycles; // In number of periods. + final double upperToValidEnd = ((validEnd - upper) / period) - upperCycles; + final boolean lowerIsBefore = (lowerToValidStart > 0); + final boolean upperIsAfter = (upperToValidEnd < 0); + if (lowerIsBefore != upperIsAfter) { + final double upperToValidStart = ((validStart - upper) / period) - upperCycles; + final double lowerToValidEnd = ((validEnd - lower) / period) - lowerCycles; + if (lowerIsBefore) { + /* + * Notation: ⎣x⎦=floor(x) and ⎡x⎤=ceil(x). + * + * We need to add an integer amount of 'period' to both sides in order to move the range + * inside the valid area. We need ⎣lowerToValidStart⎦ for reaching the point where: + * + * (validStart - period) < (new lower) ≦ validStart + * + * But we may add more because there will be no intersection without following condition: + * + * (new upper) ≧ validStart + * + * That second condition is met by ⎡upperToValidStart⎤. However adding more may cause the + * range to move the AOI completely on the right side of the domain of validity. We prevent + * that with a third condition: + * + * (new lower) < validEnd + */ + final double cycles = Math.min(Math.floor(lowerToValidEnd), + Math.max(Math.floor(lowerToValidStart), + Math.ceil (upperToValidStart))); + /* + * If after the shift we see that the following condition hold: + * + * (new lower) + period < validEnd + * + * Then we may have a situation like below: + * ┌────────────────────────────────────────────┐ + * │ Domain of validity │ + * └────────────────────────────────────────────┘ + * ┌────────────────────┐ ┌───── + * │ Area of interest │ │ AOI + * └────────────────────┘ └───── + * ↖……………………………………………………………period……………………………………………………………↗︎ + * + * The user may be requesting two extremums of the domain of validity. We can not express + * that with a single envelope. Instead, we will expand the Area Of Interest to encompass + * the full domain of validity. + */ + if (cycles + 1 < lowerToValidEnd) { + upperCycles += Math.ceil(upperToValidEnd); + } else { + upperCycles += cycles; + } lowerCycles += cycles; + } else { + /* + * Same reasoning than above with sign reverted and lower/upper variables interchanged. + * In this block, 'upperToValidEnd' and 'lowerToValidEnd' are negative, contrarily to + * above block where they were positive. + */ + final double cycles = Math.max(Math.ceil (upperToValidStart), + Math.min(Math.ceil (upperToValidEnd), + Math.floor(lowerToValidEnd))); + if (cycles - 1 > upperToValidStart) { + lowerCycles += Math.floor(lowerToValidStart); + } else { + lowerCycles += cycles; + } + upperCycles += cycles; } - upperCycles += cycles; } - } - /* - * If there is change to apply, copy the envelope when first needed and set the fields. - * If we never enter in this block, then 'areaOfInterest' will stay the envelope given - * at construction time. - */ - if (lowerCycles != 0 || upperCycles != 0) { - if (shifted == null) { - shifted = new GeneralEnvelope(areaOfInterest); + /* + * If there is change to apply, copy the envelope when first needed and set the fields. + * If we never enter in this block, then 'areaOfInterest' will stay the envelope given + * at construction time. + */ + if (lowerCycles != 0 || upperCycles != 0) { + isResultTransformed = true; + if (shifted == null) { + shifted = new GeneralEnvelope(areaOfInterest); + } + areaOfInterest = shifted; // 'shifted' may have been set before the loop. + shifted.setRange(i, lower + lowerCycles * period, // TODO: use Math.fma in JDK9. + upper + upperCycles * period); } - areaOfInterest = shifted; // 'shifted' may have been set before the loop. - geographicToAOI = projection; - shifted.setRange(i, lower + lowerCycles * period, // TODO: use Math.fma in JDK9. - upper + upperCycles * period); } } } - return this; + return Envelopes.transform(toFinal(), areaOfInterest); } /** - * Returns the (potentially shifted and/or expanded) area of interest converted by the given transform. + * Computes a position with coordinates equivalent to the given {@code pointOfInterest}, but + * potentially shifted to interior of the domain of validity specified at construction time. + * The dimensions that may be shifted are the ones having an axis with wraparound meaning. + * In order to perform this operation, the position may be temporarily converted to a geographic CRS + * and converted back to its original CRS. + * + *

The coordinate reference system should be specified in the {@code pointOfInterest}, + * or (as a fallback) in the {@code domainOfValidity} specified at construction time.

* - * @param mt a transform from the CRS of the {@code areaOfInterest} given to the constructor. - * @return the area of interest transformed by the given {@code MathTransform}. - * @throws TransformException if the transformation failed. + * @param pointOfInterest the position to potentially shift to domain of validity interior. + * If a shift is needed, then the given position will be replaced by a new position; + * the given position will not be modified. + * @return position potentially shifted to the domain of validity interior. + * @throws TransformException if a coordinate conversion failed. */ - public GeneralEnvelope result(MathTransform mt) throws TransformException { - if (geographicToAOI != null) { - mt = MathTransforms.concatenate(geographicToAOI, mt); + public DirectPosition shift(DirectPosition pointOfInterest) throws TransformException { + if (setIfNonNull(pointOfInterest.getCoordinateReferenceSystem())) { + DirectPosition shifted; + if (replaceCRS()) { + shifted = geographicToAOI.inverse().transform(pointOfInterest, null); + } else { + shifted = pointOfInterest; // To be replaced by a copy of 'pointOfInterest' when first needed. + } + final CoordinateSystem cs = crs.getCoordinateSystem(); + for (int i=cs.getDimension(); --i >= 0;) { + final double period = range(cs, i); + if (period > 0) { + transformDomainToAOI(); + final double x = shifted.getOrdinate(i); + double delta = domainOfValidity.getMinimum(i) - x; + if (delta > 0) { // Test for point before domain of validity. + delta = Math.ceil(delta / period); + } else { + delta = domainOfValidity.getMaximum(i) - x; + if (delta < 0) { // Test for point after domain of validity. + delta = Math.floor(delta / period); + } else { + continue; + } + } + if (delta != 0) { + isResultTransformed = true; + if (shifted == pointOfInterest) { + shifted = new GeneralDirectPosition(pointOfInterest); + } + pointOfInterest = shifted; // 'shifted' may have been set before the loop. + shifted.setOrdinate(i, x + delta * period); // TODO: use Math.fma in JDK9. + } + } + } } - return Envelopes.transform(mt, areaOfInterest); + return toFinal().transform(pointOfInterest, null); } } diff --git a/core/sis-referencing/src/test/java/org/apache/sis/internal/referencing/WraparoundAdjustmentTest.java b/core/sis-referencing/src/test/java/org/apache/sis/internal/referencing/WraparoundAdjustmentTest.java index ba454f5..813ccdb 100644 --- a/core/sis-referencing/src/test/java/org/apache/sis/internal/referencing/WraparoundAdjustmentTest.java +++ b/core/sis-referencing/src/test/java/org/apache/sis/internal/referencing/WraparoundAdjustmentTest.java @@ -58,13 +58,12 @@ public final strictfp class WraparoundAdjustmentTest extends TestCase { private static Envelope adjustWraparoundAxes(Envelope areaOfInterest, Envelope domainOfValidity, MathTransform validToAOI) throws TransformException { - WraparoundAdjustment adj = new WraparoundAdjustment(areaOfInterest); - adj.shiftInto(domainOfValidity, validToAOI); - return adj.result(MathTransforms.identity(2)); + WraparoundAdjustment adj = new WraparoundAdjustment(domainOfValidity, validToAOI, MathTransforms.identity(2)); + return adj.shift(areaOfInterest); } /** - * Tests {@link WraparoundAdjustment#shiftInto(Envelope, MathTransform)} + * Tests {@link WraparoundAdjustment#shift(Envelope)} * with an envelope crossing the anti-meridian. * * @throws TransformException should never happen since this test does not transform coordinates. @@ -88,7 +87,7 @@ public final strictfp class WraparoundAdjustmentTest extends TestCase { } /** - * Tests {@link WraparoundAdjustment#shiftInto(Envelope, MathTransform)} + * Tests {@link WraparoundAdjustment#shift(Envelope)} * with an envelope which is outside the grid valid area. * * @throws TransformException should never happen since this test does not transform coordinates. @@ -108,7 +107,7 @@ public final strictfp class WraparoundAdjustmentTest extends TestCase { } /** - * Tests {@link WraparoundAdjustment#shiftInto(Envelope, MathTransform)} + * Tests {@link WraparoundAdjustment#shift(Envelope)} * with an envelope shifted by 360° before or after the grid valid area. * * @throws TransformException should never happen since this test does not transform coordinates. @@ -144,7 +143,7 @@ public final strictfp class WraparoundAdjustmentTest extends TestCase { } /** - * Tests {@link WraparoundAdjustment#shiftInto(Envelope, MathTransform)} + * Tests {@link WraparoundAdjustment#shift(Envelope)} * with an envelope that cause the method to expand the area of interest. Illustration: * * {@preformat text @@ -178,7 +177,7 @@ public final strictfp class WraparoundAdjustmentTest extends TestCase { } /** - * Tests {@link WraparoundAdjustment#shiftInto(Envelope, MathTransform)} with a projected envelope. + * Tests {@link WraparoundAdjustment#shift(Envelope)} with a projected envelope. * * @throws TransformException if an error occurred while projecting a coordinate. */