This is an automated email from the ASF dual-hosted git repository. desruisseaux pushed a commit to branch geoapi-4.0 in repository https://gitbox.apache.org/repos/asf/sis.git
commit d1a38eebf9fcff0d2cf50cbc13f6b313208695b6 Author: Martin Desruisseaux <[email protected]> AuthorDate: Mon Jun 16 17:35:03 2025 +0200 Fix a `NonInvertibleTransformException` in calls to `GridGeometry.createTransformTo(…)` when the number of dimensions is not the same in source and target and some dimensions are non-linear. --- .../coverage/grid/CoordinateOperationFinder.java | 116 +++++++++----- .../sis/coverage/grid/GridDerivationTest.java | 28 ++++ .../apache/sis/coverage/grid/GridGeometryTest.java | 61 ++++++-- .../apache/sis/referencing/internal/Resources.java | 5 + .../sis/referencing/internal/Resources.properties | 1 + .../referencing/internal/Resources_fr.properties | 1 + .../apache/sis/referencing/operation/CRSPair.java | 4 +- .../operation/CoordinateOperationFinder.java | 27 ++-- .../referencing/operation/SubOperationInfo.java | 174 +++++++++++++-------- .../operation/transform/MathTransforms.java | 36 ++++- 10 files changed, 328 insertions(+), 125 deletions(-) diff --git a/endorsed/src/org.apache.sis.feature/main/org/apache/sis/coverage/grid/CoordinateOperationFinder.java b/endorsed/src/org.apache.sis.feature/main/org/apache/sis/coverage/grid/CoordinateOperationFinder.java index 52be84d555..ce1971a943 100644 --- a/endorsed/src/org.apache.sis.feature/main/org/apache/sis/coverage/grid/CoordinateOperationFinder.java +++ b/endorsed/src/org.apache.sis.feature/main/org/apache/sis/coverage/grid/CoordinateOperationFinder.java @@ -106,6 +106,14 @@ final class CoordinateOperationFinder implements Supplier<double[]> { */ private CoordinateOperation changeOfCRS; + /** + * The source <abbr>CRS</abbr> of the {@link #changeOfCRS} operation. + * May be non-null when {@link #changeOfCRS} is null if it was determined in an indirect way. + * + * @see #getSourceCRS() + */ + private CoordinateReferenceSystem sourceCRS; + /** * Whether the {@link #changeOfCRS} operation has been determined. * Note that result of determining that operation may be {@code null}, which is why we need this flag. @@ -286,6 +294,16 @@ final class CoordinateOperationFinder implements Supplier<double[]> { isWraparoundDisabled = true; } + /** + * Returns the source <abbr>CRS</abbr> of the {@link #changeOfCRS} transform. + */ + private CoordinateReferenceSystem getSourceCRS() throws FactoryException, TransformException { + if (sourceCRS == null) { + sourceCRS = changeOfCRS().getSourceCRS(); + } + return sourceCRS; + } + /** * 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. @@ -309,31 +327,9 @@ final class CoordinateOperationFinder implements Supplier<double[]> { */ private CoordinateOperation changeOfCRS() throws FactoryException, TransformException { if (!knowChangeOfCRS) { - final Envelope sourceEnvelope = source.envelope; - final Envelope targetEnvelope = target.envelope; try { CoordinateOperations.CONSTANT_COORDINATES.set(this); - if (sourceEnvelope != null && targetEnvelope != null) { - changeOfCRS = Envelopes.findOperation(sourceEnvelope, targetEnvelope); - } - if (changeOfCRS == null && source.isDefined(GridGeometry.CRS) && target.isDefined(GridGeometry.CRS)) { - final CoordinateReferenceSystem sourceCRS = source.getCoordinateReferenceSystem(); - /* - * Unconditionally create operation even if CRS are the same. A non-null operation trig - * the check for wraparound axes, which is necessary even if the transform is identity. - */ - DefaultGeographicBoundingBox areaOfInterest = null; - if (sourceEnvelope != null || targetEnvelope != null) try { - areaOfInterest = new DefaultGeographicBoundingBox(); - areaOfInterest.setBounds(targetEnvelope != null ? targetEnvelope : sourceEnvelope); - } catch (TransformException e) { - areaOfInterest = null; - recoverableException("changeOfCRS", e); - } - changeOfCRS = CRS.findOperation(sourceCRS, target.getCoordinateReferenceSystem(), areaOfInterest); - } - } catch (BackingStoreException e) { // May be thrown by getConstantCoordinates(). - throw e.unwrapOrRethrow(TransformException.class); + changeOfCRS = changeOfCRS(source, target); } finally { CoordinateOperations.CONSTANT_COORDINATES.remove(); } @@ -342,6 +338,47 @@ final class CoordinateOperationFinder implements Supplier<double[]> { return changeOfCRS; } + /** + * Computes the change of <abbr>CRS</abbr> between the given pair of grid geometries. + * + * @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. + * @return coordinate operation from source to target CRS, or {@code null} if none. + * @throws FactoryException if no operation can be found between the source and target CRS. + * @throws TransformException if some coordinates cannot be transformed to the specified target. + */ + private static CoordinateOperation changeOfCRS(final GridGeometry source, final GridGeometry target) + throws FactoryException, TransformException + { + CoordinateOperation changeOfCRS = null; + final Envelope sourceEnvelope = source.envelope; + final Envelope targetEnvelope = target.envelope; + try { + if (sourceEnvelope != null && targetEnvelope != null) { + changeOfCRS = Envelopes.findOperation(sourceEnvelope, targetEnvelope); + } + if (changeOfCRS == null && source.isDefined(GridGeometry.CRS) && target.isDefined(GridGeometry.CRS)) { + final CoordinateReferenceSystem sourceCRS = source.getCoordinateReferenceSystem(); + /* + * Unconditionally create operation even if CRS are the same. A non-null operation trig + * the check for wraparound axes, which is necessary even if the transform is identity. + */ + DefaultGeographicBoundingBox areaOfInterest = null; + if (sourceEnvelope != null || targetEnvelope != null) try { + areaOfInterest = new DefaultGeographicBoundingBox(); + areaOfInterest.setBounds(targetEnvelope != null ? targetEnvelope : sourceEnvelope); + } catch (TransformException e) { + areaOfInterest = null; + recoverableException("changeOfCRS", e); + } + changeOfCRS = CRS.findOperation(sourceCRS, target.getCoordinateReferenceSystem(), areaOfInterest); + } + } catch (BackingStoreException e) { // May be thrown by getConstantCoordinates(). + throw e.unwrapOrRethrow(TransformException.class); + } + return changeOfCRS; + } + /** * 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. @@ -410,8 +447,8 @@ apply: if (forwardChangeOfCRS == null) { targetMedian = sourceMedian; sourceMedian = null; } - final WraparoundApplicator ap = new WraparoundApplicator(sourceMedian, - targetMedian, changeOfCRS.getTargetCRS().getCoordinateSystem()); + final var ap = new WraparoundApplicator(sourceMedian, targetMedian, + changeOfCRS.getTargetCRS().getCoordinateSystem()); forwardChangeOfCRS = ap.forDomainOfUse(forwardChangeOfCRS); } } @@ -432,13 +469,23 @@ apply: if (forwardChangeOfCRS == null) { */ final MathTransform inverse() throws FactoryException, TransformException { final MathTransform sourceCrsToGrid = source.getGridToCRS(anchor).inverse(); - @SuppressWarnings("LocalVariableHidesMemberVariable") - final CoordinateOperation changeOfCRS = changeOfCRS(); - if (changeOfCRS == null) { - return sourceCrsToGrid; - } if (inverseChangeOfCRS == null) { - inverseChangeOfCRS = changeOfCRS.getMathTransform().inverse(); + if (knowChangeOfCRS) { + @SuppressWarnings("LocalVariableHidesMemberVariable") + final CoordinateOperation changeOfCRS = changeOfCRS(); + if (changeOfCRS != null) { + inverseChangeOfCRS = changeOfCRS.getMathTransform().inverse(); + } + } else { + final CoordinateOperation inverse = changeOfCRS(target, source); + if (inverse != null) { + sourceCRS = inverse.getTargetCRS(); + inverseChangeOfCRS = inverse.getMathTransform(); + } + } + if (inverseChangeOfCRS == null) { + return sourceCrsToGrid; + } if (!isWraparoundDisabled) { isWraparoundApplied = false; if (!isWraparoundNeedVerified) { @@ -605,8 +652,7 @@ apply: if (forwardChangeOfCRS == null) { } if (sourceMedian != null) { final MathTransform inverseNoWrap = inverseChangeOfCRS; - final WraparoundApplicator ap = new WraparoundApplicator(targetMedian, - sourceMedian, changeOfCRS().getSourceCRS().getCoordinateSystem()); + final var ap = new WraparoundApplicator(targetMedian, sourceMedian, getSourceCRS().getCoordinateSystem()); inverseChangeOfCRS = ap.forDomainOfUse(inverseNoWrap); if (inverseChangeOfCRS != inverseNoWrap) { crsToGrid = MathTransforms.concatenate(inverseChangeOfCRS, sourceCrsToGrid); @@ -671,8 +717,8 @@ apply: if (forwardChangeOfCRS == null) { * Invoked when the target CRS has some dimensions that the source CRS does not have. * For example, this is invoked during the conversion from (<var>x</var>, <var>y</var>) * coordinates to (<var>x</var>, <var>y</var>, <var>t</var>). If constant values can - * be given to the missing dimensions, than those values are returned. Otherwise this - * method returns {@code null}. + * be given to the missing dimensions, then those values are returned. + * Otherwise this method returns {@code null}. * * <p>The returned array has a length equals to the number of dimensions in the target CRS. * Only coordinates in dimensions without source (<var>t</var> in the above example) will be used. diff --git a/endorsed/src/org.apache.sis.feature/test/org/apache/sis/coverage/grid/GridDerivationTest.java b/endorsed/src/org.apache.sis.feature/test/org/apache/sis/coverage/grid/GridDerivationTest.java index d6055afc6f..86fe67d11e 100644 --- a/endorsed/src/org.apache.sis.feature/test/org/apache/sis/coverage/grid/GridDerivationTest.java +++ b/endorsed/src/org.apache.sis.feature/test/org/apache/sis/coverage/grid/GridDerivationTest.java @@ -357,6 +357,34 @@ public final class GridDerivationTest extends TestCase { assertExtentEquals(new long[] {15, 3, 0}, new long[] {49, 9, 5}, grid.getExtent()); } + /** + * Tests {@link GridDerivation#subgrid(GridGeometry)} using a grid with less dimensions + * than the source grid geometry. + * + * @see <a href="https://issues.apache.org/jira/browse/SIS-610";>SIS-610</a> + */ + @Test + public void testSubgridWithLessDimensions() { + var envelope = new GeneralEnvelope(HardCodedCRS.WGS84_4D); + envelope.setRange(0, 10, 20); + envelope.setRange(1, 30, 40); + envelope.setRange(2, 2, 4); + envelope.setRange(3, 3, 6); + var env2D = new Envelope2D(HardCodedCRS.WGS84, 10, 30, 10, 10); + var extent = new GridExtent(null, null, new long[] {2, 2, 1, 1}, false); + var grid = new GridGeometry(extent, envelope, GridOrientation.DISPLAY); + var aoi = new GridGeometry(new GridExtent(2, 2), env2D, GridOrientation.DISPLAY); + GridGeometry slice = grid.derive().subgrid(aoi).build(); + Matrix gridToCRS = MathTransforms.getMatrix(slice.getGridToCRS(PixelInCell.CELL_CORNER)); + Matrix expected = Matrices.create(5, 5, new double[] { + 5, 0, 0, 0, 10, + 0, -5, 0, 0, 40, + 0, 0, 2, 0, 2, + 0, 0, 0, 3, 3, + 0, 0, 0, 0, 1}); + assertMatrixEquals(expected, gridToCRS, STRICT, "gridToCRS"); + } + /** * Tests {@link GridDerivation#subgrid(Envelope, double...)} on a grid using a polar projection. * The test also uses a geographic envelope with more dimensions than the source grid geometry. diff --git a/endorsed/src/org.apache.sis.feature/test/org/apache/sis/coverage/grid/GridGeometryTest.java b/endorsed/src/org.apache.sis.feature/test/org/apache/sis/coverage/grid/GridGeometryTest.java index 603e8c6259..b423c0b2c2 100644 --- a/endorsed/src/org.apache.sis.feature/test/org/apache/sis/coverage/grid/GridGeometryTest.java +++ b/endorsed/src/org.apache.sis.feature/test/org/apache/sis/coverage/grid/GridGeometryTest.java @@ -277,7 +277,8 @@ public final class GridGeometryTest extends TestCase { 1, 0, 0.5, 0, 1, 0.5, 0, 0, 1)); - final GridGeometry grid = new GridGeometry(extent, PixelInCell.CELL_CENTER, gridToCRS, null); + + final var grid = new GridGeometry(extent, PixelInCell.CELL_CENTER, gridToCRS, null); assertTrue(grid.getGridToCRS(PixelInCell.CELL_CORNER).isIdentity(), "gridToCRS.isIdentity"); verifyGridToCRS(grid); } @@ -324,7 +325,8 @@ public final class GridGeometryTest extends TestCase { 0, 0.5, -90, 0.5, 0, -180, 0, 0, 1)); - final GridGeometry grid = new GridGeometry(PixelInCell.CELL_CORNER, gridToCRS, envelope, GridRoundingMode.NEAREST); + + final var grid = new GridGeometry(PixelInCell.CELL_CORNER, gridToCRS, envelope, GridRoundingMode.NEAREST); assertExtentEquals( new long[] {370, 40}, new long[] {389, 339}, grid.getExtent()); @@ -362,7 +364,7 @@ public final class GridGeometryTest extends TestCase { * Simplest case: no axis flip. * Verification: y = 2 × −25 + 40 = −10 (the minimum value declared in envelope). */ - GridGeometry grid = new GridGeometry(extent, aoi, GridOrientation.HOMOTHETY); + var grid = new GridGeometry(extent, aoi, GridOrientation.HOMOTHETY); Matrix matrix = MathTransforms.getMatrix(grid.getGridToCRS(PixelInCell.CELL_CORNER)); assertMatrixEquals(new Matrix3(0.5, 0, 50, 0, 2, 40, @@ -421,7 +423,7 @@ public final class GridGeometryTest extends TestCase { * Same case as the one tested by `testFromExtentAndEnvelope()`, * but with axis order swapped. */ - GridGeometry grid = new GridGeometry(extent, aoi, GridOrientation.DISPLAY); + var grid = new GridGeometry(extent, aoi, GridOrientation.DISPLAY); Matrix matrix = MathTransforms.getMatrix(grid.getGridToCRS(PixelInCell.CELL_CORNER)); assertMatrixEquals(new Matrix3(0, 0.5, 50, -2, 0, 20, @@ -528,23 +530,23 @@ public final class GridGeometryTest extends TestCase { public void testUpsample() { final GridGeometry grid; { // Source grid - final GridExtent extent = new GridExtent(null, new long[] {10,-8}, new long[] {100, 50}, false); - final Matrix3 mat = new Matrix3( + var extent = new GridExtent(null, new long[] {10,-8}, new long[] {100, 50}, false); + var mat = new Matrix3( 1, 0, 10, 0, -2, 50, 0, 0, 1); - final MathTransform gridToCRS = MathTransforms.linear(mat); + MathTransform gridToCRS = MathTransforms.linear(mat); grid = new GridGeometry(extent, PixelInCell.CELL_CENTER, gridToCRS, HardCodedCRS.CARTESIAN_2D); } final GridGeometry upsampled = grid.upsample(4, 4); final GridGeometry expected; { // Expected grid - GridExtent extent = new GridExtent(null, new long[] {40,-32}, new long[] {400, 200}, false); - final Matrix3 mat = new Matrix3( + var extent = new GridExtent(null, new long[] {40,-32}, new long[] {400, 200}, false); + var mat = new Matrix3( 0.25, 0, 9.625, 0, -0.5, 50.750, 0, 0, 1); - final MathTransform gridToCRS = MathTransforms.linear(mat); + MathTransform gridToCRS = MathTransforms.linear(mat); expected = new GridGeometry(extent, PixelInCell.CELL_CENTER, gridToCRS, HardCodedCRS.CARTESIAN_2D); } assertSame(grid.getEnvelope(), upsampled.getEnvelope(), "envelope"); @@ -559,7 +561,7 @@ public final class GridGeometryTest extends TestCase { */ @Test public void testShiftGrid() { - GridGeometry grid = new GridGeometry( + final var grid = new GridGeometry( new GridExtent(17, 10), PixelInCell.CELL_CENTER, MathTransforms.linear(new Matrix3( @@ -796,6 +798,43 @@ public final class GridGeometryTest extends TestCase { MathTransforms.getMatrix(tr), STRICT, "createTransformTo"); } + /** + * Tests {@link GridGeometry#createTransformTo(GridGeometry, PixelInCell)} with a change + * in the number of dimensions. One of the dimensions is non-linear. + * + * @throws TransformException if the transform cannot be computed. + */ + @Test + public void testCreateTransformWithDimensionChange() throws TransformException { + final var source = new GridGeometry( + new GridExtent(null, null, new long[] {17, 10, 6}, false), + PixelInCell.CELL_CENTER, + MathTransforms.compound( + MathTransforms.linear(new Matrix3( + 1, 0, -7.0, + 0, -1, 50.0, + 0, 0, 1)), + MathTransforms.interpolate(null, new double[] {3, 4, 9})), + HardCodedCRS.WGS84_WITH_TIME); + + final var target = new GridGeometry( + new GridExtent(200, 300), + PixelInCell.CELL_CENTER, + MathTransforms.linear(new Matrix3( + -0.05, 0, 53.0, + 0, 0.1, -8.0, + 0, 0, 1)), + HardCodedCRS.WGS84_LATITUDE_FIRST); + + final Matrix expected = Matrices.create(3, 4, new double[] { + 0, 20, 0, 60, + 10, 0, 0, 10, + 0, 0, 0, 1}); + + final MathTransform tr = source.createTransformTo(target, PixelInCell.CELL_CENTER); + assertMatrixEquals(expected, MathTransforms.getMatrix(tr), STRICT, "createTransformTo"); + } + /** * Tests {@link GridGeometry#contains(GridGeometry)} and {@link GridGeometry#intersects(GridGeometry)}. */ diff --git a/endorsed/src/org.apache.sis.referencing/main/org/apache/sis/referencing/internal/Resources.java b/endorsed/src/org.apache.sis.referencing/main/org/apache/sis/referencing/internal/Resources.java index 96af716fd4..34eeeae90d 100644 --- a/endorsed/src/org.apache.sis.referencing/main/org/apache/sis/referencing/internal/Resources.java +++ b/endorsed/src/org.apache.sis.referencing/main/org/apache/sis/referencing/internal/Resources.java @@ -187,6 +187,11 @@ public class Resources extends IndexedResourceBundle { */ public static final short ConformanceMeansDatumShift = 11; + /** + * A constant value is required for the coordinate at index {0}. + */ + public static final short ConstantCoordinateValueRequired_1 = 109; + /** * This parameter is shown for completeness, but should never have a value different than {0} * for this projection. diff --git a/endorsed/src/org.apache.sis.referencing/main/org/apache/sis/referencing/internal/Resources.properties b/endorsed/src/org.apache.sis.referencing/main/org/apache/sis/referencing/internal/Resources.properties index edb2538662..74149a92a1 100644 --- a/endorsed/src/org.apache.sis.referencing/main/org/apache/sis/referencing/internal/Resources.properties +++ b/endorsed/src/org.apache.sis.referencing/main/org/apache/sis/referencing/internal/Resources.properties @@ -72,6 +72,7 @@ CanNotTransformEnvelopeToGeodetic = Cannot transform envelope to a geodetic refe CanNotTransformGeometry = Cannot transform the given geometry. CanNotUseGeodeticParameters_2 = Cannot use the {0} geodetic parameters. Caused by: {1} ColinearAxisDirections_2 = Axis directions {0} and {1} are colinear. +ConstantCoordinateValueRequired_1 = A constant value is required for the coordinate at index {0}. CoordinateOperationNotFound_2 = Coordinate operation from system \u201c{0}\u201d to \u201c{1}\u201d has not been found. DatumChangesDirectory_1 = Datum shift files are searched in the \u201c{0}\u201d directory. DatumOriginShallBeDate = Origin of temporal datum shall be a date. diff --git a/endorsed/src/org.apache.sis.referencing/main/org/apache/sis/referencing/internal/Resources_fr.properties b/endorsed/src/org.apache.sis.referencing/main/org/apache/sis/referencing/internal/Resources_fr.properties index 94495000ef..e63ca97752 100644 --- a/endorsed/src/org.apache.sis.referencing/main/org/apache/sis/referencing/internal/Resources_fr.properties +++ b/endorsed/src/org.apache.sis.referencing/main/org/apache/sis/referencing/internal/Resources_fr.properties @@ -77,6 +77,7 @@ CanNotTransformEnvelopeToGeodetic = Ne peut pas transformer l\u2019enveloppe ver CanNotTransformGeometry = Ne peut pas transformer la g\u00e9om\u00e9trie donn\u00e9e. CanNotUseGeodeticParameters_2 = Ne peut pas utiliser les param\u00e8tres g\u00e9od\u00e9siques {0}. La cause est\u202f: {1} ColinearAxisDirections_2 = Les directions d\u2019axes {0} et {1} sont colin\u00e9aires. +ConstantCoordinateValueRequired_1 = Une valeur constante est requise pour la coordonn\u00e9e \u00e0 l\u2019index {0}. CoordinateOperationNotFound_2 = L\u2019op\u00e9ration sur les coordonn\u00e9es du syst\u00e8me \u00ab\u202f{0}\u202f\u00bb vers \u00ab\u202f{1}\u202f\u00bb n\u2019a pas \u00e9t\u00e9 trouv\u00e9e. DatumChangesDirectory_1 = Les fichiers de changements de r\u00e9f\u00e9rentiel sont cherch\u00e9s dans le dossier \u00ab\u202f{0}\u202f\u00bb. DatumOriginShallBeDate = L\u2019origine d\u2019un r\u00e9f\u00e9rentiel temporel doit \u00eatre une date. diff --git a/endorsed/src/org.apache.sis.referencing/main/org/apache/sis/referencing/operation/CRSPair.java b/endorsed/src/org.apache.sis.referencing/main/org/apache/sis/referencing/operation/CRSPair.java index 0e46001355..3df3ecb118 100644 --- a/endorsed/src/org.apache.sis.referencing/main/org/apache/sis/referencing/operation/CRSPair.java +++ b/endorsed/src/org.apache.sis.referencing/main/org/apache/sis/referencing/operation/CRSPair.java @@ -35,7 +35,7 @@ import org.apache.sis.util.privy.Strings; * * @author Martin Desruisseaux (Geomatys) */ -final class CRSPair { +final class CRSPair extends org.apache.sis.pending.jdk.Record { /** * The source and target CRS. */ @@ -70,7 +70,7 @@ final class CRSPair { @Override public boolean equals(final Object object) { if (object instanceof CRSPair) { - final CRSPair that = (CRSPair) object; + final var that = (CRSPair) object; return Objects.equals(this.sourceCRS, that.sourceCRS) && Objects.equals(this.targetCRS, that.targetCRS); } diff --git a/endorsed/src/org.apache.sis.referencing/main/org/apache/sis/referencing/operation/CoordinateOperationFinder.java b/endorsed/src/org.apache.sis.referencing/main/org/apache/sis/referencing/operation/CoordinateOperationFinder.java index 2eb7c5a3fa..353eaaa8db 100644 --- a/endorsed/src/org.apache.sis.referencing/main/org/apache/sis/referencing/operation/CoordinateOperationFinder.java +++ b/endorsed/src/org.apache.sis.referencing/main/org/apache/sis/referencing/operation/CoordinateOperationFinder.java @@ -809,6 +809,8 @@ public class CoordinateOperationFinder extends CoordinateOperationRegistry { final SubOperationInfo[] infos; try { infos = SubOperationInfo.createSteps(this, sourceComponents, targetComponents); + } catch (NoninvertibleTransformException e) { + throw new OperationNotFoundException(notFoundMessage(sourceCRS, targetCRS), e); } catch (TransformException e) { throw new FactoryException(notFoundMessage(sourceCRS, targetCRS), e); } @@ -845,10 +847,10 @@ public class CoordinateOperationFinder extends CoordinateOperationRegistry { operation = createFromAffineTransform(AXIS_CHANGES, sourceCRS, stepSourceCRS, null, select); } /* - * For each sub-operation, create a PassThroughOperation for the (stepSourceCRS → stepTargetCRS) operation. - * Each source CRS inside this loop will be for dimensions at indices [startAtDimension … endAtDimension-1]. + * For each sub-operation, create a `PassThroughOperation` for the (`stepSourceCRS` → `stepTargetCRS`) operation. + * Each source CRS in the loop will be for dimensions at indices [`firstAffectedCoordinate` … `endAtDimension`-1]. * Note that those indices are not necessarily the same as the indices in the fields of the same name in - * SubOperationInfo, because those indices are not relative to the same CompoundCRS. + * `SubOperationInfo`, because those indices are not relative to the same `CompoundCRS`. */ int endAtDimension = 0; int remainingSourceDimensions = select.getNumRow() - 1; @@ -874,13 +876,18 @@ public class CoordinateOperationFinder extends CoordinateOperationRegistry { stepTargetCRS = createCompoundCRS(target, stepComponents); } int delta = source.getCoordinateSystem().getDimension(); - final int startAtDimension = endAtDimension; + final int firstAffectedCoordinate = endAtDimension; endAtDimension += delta; final int numTrailingCoordinates = remainingSourceDimensions - endAtDimension; CoordinateOperation subOperation = info.operation; - if ((startAtDimension | numTrailingCoordinates) != 0) { - subOperation = new DefaultPassThroughOperation(IdentifiedObjects.getProperties(subOperation), - stepSourceCRS, stepTargetCRS, subOperation, startAtDimension, numTrailingCoordinates); + if ((firstAffectedCoordinate | numTrailingCoordinates) != 0) { + subOperation = new DefaultPassThroughOperation( + IdentifiedObjects.getProperties(subOperation), + stepSourceCRS, + stepTargetCRS, + subOperation, + firstAffectedCoordinate, + numTrailingCoordinates); } /* * Concatenate the operation with the ones we have found so far, and use the current `stepTargetCRS` @@ -899,11 +906,13 @@ public class CoordinateOperationFinder extends CoordinateOperationRegistry { * source dimensions, add those constants as the last step. It never occur except is some particular * contexts like when computing a transform between two `GridGeometry` instances. */ - if (stepComponents.length < infos.length) { - final Matrix m = SubOperationInfo.createConstantOperation(infos, stepComponents.length, + if (stepComponents.length < infos.length) try { + final Matrix m = SubOperationInfo.createConstantOperation(infos, stepSourceCRS.getCoordinateSystem().getDimension(), targetCRS.getCoordinateSystem().getDimension()); operation = concatenate(operation, createFromAffineTransform(CONSTANTS, stepSourceCRS, targetCRS, null, m)); + } catch (TransformException e) { + throw new FactoryException(notFoundMessage(sourceCRS, targetCRS), e); } return asList(operation); } diff --git a/endorsed/src/org.apache.sis.referencing/main/org/apache/sis/referencing/operation/SubOperationInfo.java b/endorsed/src/org.apache.sis.referencing/main/org/apache/sis/referencing/operation/SubOperationInfo.java index 277276cd67..75f660aebb 100644 --- a/endorsed/src/org.apache.sis.referencing/main/org/apache/sis/referencing/operation/SubOperationInfo.java +++ b/endorsed/src/org.apache.sis.referencing/main/org/apache/sis/referencing/operation/SubOperationInfo.java @@ -18,12 +18,16 @@ package org.apache.sis.referencing.operation; import java.util.List; import org.opengis.referencing.crs.*; +import org.opengis.referencing.operation.Matrix; import org.opengis.referencing.operation.CoordinateOperation; import org.opengis.referencing.operation.OperationNotFoundException; +import org.opengis.referencing.operation.NoninvertibleTransformException; import org.opengis.referencing.operation.TransformException; import org.opengis.util.FactoryException; +import org.apache.sis.referencing.internal.Resources; import org.apache.sis.referencing.operation.matrix.Matrices; import org.apache.sis.referencing.operation.matrix.MatrixSIS; +import org.apache.sis.referencing.operation.transform.MathTransforms; // Specific to the geoapi-4.0 branch: import org.apache.sis.referencing.crs.DefaultImageCRS; @@ -82,32 +86,29 @@ final class SubOperationInfo { /** * The constant values to store in target coordinates, or {@code null} if none. This array is usually null. - * It may be non-null if no source CRS component has been found for a target CRS component. For example, if - * the source CRS has (<var>x</var>, <var>y</var>) axes and the target CRS has (<var>x</var>, <var>y</var>, - * <var>t</var>) axes, then this array may be set to a non-null value for specifying the <var>t</var> value. - * The array length is the number of dimensions in the full (usually compound) target CRS, but only the - * coordinate values for sourceless dimensions are used. Other coordinates are ignored and can be NaN. + * It may be non-null if no source CRS component has been found for a target CRS component. * Exactly one of {@link #operation} or {@link #constants} shall be non-null. * + * <p>The array length is the number of dimensions in the full (usually compound) target CRS, but only the + * coordinate values for sourceless dimensions are used. Other coordinates are ignored and can be NaN.</p> + * + * <h4>Example</h4> + * If the source <abbr>CRS</abbr> has (<var>x</var>, <var>y</var>) axes and the target <abbr>CRS</abbr> + * has (<var>x</var>, <var>y</var>, <var>t</var>) axes, then this array may be set to a non-null value + * for specifying the <var>t</var> value. + * * @see CoordinateOperationContext#getConstantCoordinates() */ private final double[] constants; /** * The first dimension (inclusive) and the last dimension (exclusive) where the {@link SingleCRS} starts/ends - * in the full (usually compound) CRS. It may be an index in source dimensions or in target dimensions, - * depending on the following rules: - * - * <ul> - * <li>If {@link #operation} is non-null, then this is an index in the <em>source</em> dimensions. - * This is the normal case.</li> - * <li>Otherwise the operation is sourceless with coordinate results described by the {@link #constants} array. - * Since there is no source, the index become an index in target dimensions instead of source dimensions.</li> - * </ul> + * in the full (usually compound) CRS. * * @see #sourceToSelected(int, SubOperationInfo[]) */ - private final int startAtDimension, endAtDimension; + private final int sourceLowerDimension, sourceUpperDimension, + targetLowerDimension, targetUpperDimension; /** * Index of this instance in the array of {@code SubOperationInfo} instances, @@ -119,13 +120,17 @@ final class SubOperationInfo { * Creates a new instance wrapping the given coordinate operation or coordinate constants. * Exactly one of {@code operation} or {@code constants} shall be non-null. */ - private SubOperationInfo(final int targetComponentIndex, final CoordinateOperation operation, - final double[] constants, final int startAtDimension, final int endAtDimension) + private SubOperationInfo(final CoordinateOperation operation, final double[] constants, + final int sourceLowerDimension, final int sourceUpperDimension, + final int targetLowerDimension, final int targetUpperDimension, + final int targetComponentIndex) { this.operation = operation; this.constants = constants; - this.startAtDimension = startAtDimension; - this.endAtDimension = endAtDimension; + this.sourceLowerDimension = sourceLowerDimension; + this.sourceUpperDimension = sourceUpperDimension; + this.targetLowerDimension = targetLowerDimension; + this.targetUpperDimension = targetUpperDimension; this.targetComponentIndex = targetComponentIndex; assert (operation == null) != (constants == null); } @@ -140,29 +145,26 @@ final class SubOperationInfo { * @param targets all components of the target CRS. * @return information about each coordinate operation from a source CRS to a target CRS, or {@code null}. * @throws FactoryException if an error occurred while grabbing a coordinate operation. - * @throws TransformException if an error occurred while computing the sourceless coordinate constants. + * @throws TransformException if an error occurred while computing a coordinate value for unmatched dimension. + * @throws NoninvertibleTransformException if the default coordinate value for a unmatched dimension is NaN. */ static SubOperationInfo[] createSteps(final CoordinateOperationFinder caller, final List<? extends SingleCRS> sources, final List<? extends SingleCRS> targets) throws FactoryException, TransformException { - final SubOperationInfo[] infos = new SubOperationInfo[targets.size()]; - final boolean[] sourceComponentIsUsed = new boolean[sources.size()]; + final var infos = new SubOperationInfo[targets.size()]; + final var sourceComponentIsUsed = new boolean[sources.size()]; /* * Iterate over target CRS because all of them must have an operation. * One the other hand, source CRS can be used zero or one (not two) time. - * - * Note: the loops on source CRS and on target CRS follow the same pattern. - * The first 6 lines of code are the same with only `target` replaced - * by `source`. */ - int targetStartAtDimension; - int targetEndAtDimension = 0; + int targetLowerDimension; + int targetUpperDimension = 0; next: for (int targetComponentIndex = 0; targetComponentIndex < infos.length; targetComponentIndex++) { final SingleCRS target = targets.get(targetComponentIndex); - targetStartAtDimension = targetEndAtDimension; - targetEndAtDimension += target.getCoordinateSystem().getDimension(); + targetLowerDimension = targetUpperDimension; + targetUpperDimension += target.getCoordinateSystem().getDimension(); final Class<?> targetType = type(target); OperationNotFoundException failure = null; @@ -171,17 +173,18 @@ next: for (int targetComponentIndex = 0; targetComponentIndex < infos.length; * Some sources may be left unused after this method completion. Check only * sources that may be compatible according the `COMPATIBLE_TYPES` array. */ - int sourceStartAtDimension; - int sourceEndAtDimension = 0; + int sourceLowerDimension; + int sourceUpperDimension = 0; for (int sourceComponentIndex = 0; sourceComponentIndex < sourceComponentIsUsed.length; sourceComponentIndex++) { final SingleCRS source = sources.get(sourceComponentIndex); - sourceStartAtDimension = sourceEndAtDimension; - sourceEndAtDimension += source.getCoordinateSystem().getDimension(); + sourceLowerDimension = sourceUpperDimension; + sourceUpperDimension += source.getCoordinateSystem().getDimension(); if (!sourceComponentIsUsed[sourceComponentIndex]) { - for (final Class<?>[] sourceTypes : COMPATIBLE_TYPES) { - if (sourceTypes[0].isAssignableFrom(targetType)) { - for (final Class<?> sourceType : sourceTypes) { - if (sourceType.isAssignableFrom(type(source))) { + final Class<?> sourceType = type(source); + for (final Class<?>[] compatibleTypes : COMPATIBLE_TYPES) { + if (compatibleTypes[0].isAssignableFrom(targetType)) { + for (final Class<?> compatibleType : compatibleTypes) { + if (compatibleType.isAssignableFrom(sourceType)) { final CoordinateOperation operation; try { operation = caller.createOperation(source, target); @@ -194,8 +197,8 @@ next: for (int targetComponentIndex = 0; targetComponentIndex < infos.length; continue; } /* - * Found an operation. Exclude the source component from the list because each source - * should be used at most once by SourceComponent. Note that the same source may still + * Found an operation. Exclude the source component from the list because each source + * should be used at most once by `SubOperationInfo`. Note that the same source may * be used again in another context if that source is also an interpolation CRS. * * EXAMPLE: consider a coordinate operation from (GeodeticCRS₁, VerticalCRS₁) source @@ -204,13 +207,16 @@ next: for (int targetComponentIndex = 0; targetComponentIndex < infos.length; * to VerticalCRS₂. But the operation on vertical coordinates may need GeodeticCRS₁ * for doing its work, so GeodeticCRS₁ is needed twice. However, when needed for the * vertical coordinate operation, the GeodeticCRS₁ is used as an "interpolation CRS". - * Interpolation CRS are handled in other code paths; it is not the business of this - * SourceComponent class to care about them. From the point of view of this class, + * Interpolation CRS are handled in other code paths, it is not the business of this + * `SubOperationInfo` class to care about them. From the point of view of this class, * GeodeticCRS₁ is used only once. */ sourceComponentIsUsed[sourceComponentIndex] = true; - infos[targetComponentIndex] = new SubOperationInfo(targetComponentIndex, - operation, null, sourceStartAtDimension, sourceEndAtDimension); + infos[targetComponentIndex] = new SubOperationInfo( + operation, null, + sourceLowerDimension, sourceUpperDimension, + targetLowerDimension, targetUpperDimension, + targetComponentIndex); if (failure != null) { CoordinateOperationRegistry.recoverableException("decompose", failure); @@ -227,21 +233,27 @@ next: for (int targetComponentIndex = 0; targetComponentIndex < infos.length; } /* * If we reach this point, we have not been able to find a source CRS that we can map to the target CRS. - * Usually this is fatal; returning null will instruct the caller to throw `OperationNotFoundException`. + * Usually this is fatal, returning null will instruct the caller to throw `OperationNotFoundException`. * However, in some contexts (e.g. when searching for an operation between two `GridGeometry` instances) * it is possible to assign a constant value to the target coordinates. Those values cannot be guessed * by `org.apache.sis.referencing`; they must be provided by caller. If such constants are specified, * then we will try to apply them. */ final double[] constants = CoordinateOperationContext.getConstantCoordinates(); - if (constants == null || constants.length < targetEndAtDimension) { + if (constants == null || constants.length < targetUpperDimension) { return null; } - for (int i = targetStartAtDimension; i < targetEndAtDimension; i++) { - if (Double.isNaN(constants[i])) return null; + for (int i = targetLowerDimension; i < targetUpperDimension; i++) { + if (Double.isNaN(constants[i])) { + throw new NoninvertibleTransformException(caller.resources() + .getString(Resources.Keys.ConstantCoordinateValueRequired_1, i)); + } } - infos[targetComponentIndex] = new SubOperationInfo(targetComponentIndex, - null, constants, targetStartAtDimension, targetEndAtDimension); + infos[targetComponentIndex] = new SubOperationInfo( + null, constants, + sourceUpperDimension, sourceUpperDimension, + targetLowerDimension, targetUpperDimension, + targetComponentIndex); } return infos; } @@ -249,8 +261,8 @@ next: for (int targetComponentIndex = 0; targetComponentIndex < infos.length; /** * Returns the source CRS of given operations. This method modifies the given array in-place by moving all * sourceless operations last. Then an array is returned with the source CRS of only ordinary operations. - * Each CRS at index <var>i</var> in the returned array is the component from {@link #startAtDimension} - * inclusive to {@link #endAtDimension} exclusive in the complete (usually compound) source CRS analyzed + * Each CRS at index <var>i</var> in the returned array is the component from {@link #sourceLowerDimension} + * inclusive to {@link #sourceUpperDimension} exclusive in the complete (usually compound) source CRS analyzed * by {@link CoordinateOperationFinder}. * * @param selected all operations from source to target {@link CompoundCRS}. @@ -267,7 +279,7 @@ next: for (int targetComponentIndex = 0; targetComponentIndex < infos.length; n--; } } - final CoordinateReferenceSystem[] stepComponents = new CoordinateReferenceSystem[n]; + final var stepComponents = new CoordinateReferenceSystem[n]; for (int i=0; i<n; i++) { stepComponents[i] = selected[i].operation.getSourceCRS(); } @@ -300,7 +312,7 @@ next: for (int targetComponentIndex = 0; targetComponentIndex < infos.length; * the source CRS are not necessarily picked in the same order as they appear in the list. * * <h4>Example</h4> - * if the source CRS has (<var>x</var>, <var>y</var>, <var>t</var>) coordinates and the target CRS has + * If the source CRS has (<var>x</var>, <var>y</var>, <var>t</var>) coordinates and the target CRS has * (<var>t</var>, <var>x</var>, <var>y</var>) coordinates with some operation applied on <var>x</var> * and <var>y</var>, then the operations will be applied in that order: * @@ -332,13 +344,14 @@ next: for (int targetComponentIndex = 0; targetComponentIndex < infos.length; int selectedDimensions = 0; for (final SubOperationInfo component : selected) { if (component.operation == null) break; - selectedDimensions += component.endAtDimension - component.startAtDimension; + selectedDimensions += component.sourceUpperDimension - component.sourceLowerDimension; } final MatrixSIS select = Matrices.createZero(selectedDimensions + 1, sourceDimensions + 1); select.setElement(selectedDimensions, sourceDimensions, 1); int j = 0; for (final SubOperationInfo component : selected) { - for (int i=component.startAtDimension; i<component.endAtDimension; i++) { + if (component.operation == null) break; + for (int i = component.sourceLowerDimension; i < component.sourceUpperDimension; i++) { select.setElement(j++, i, 1); } } @@ -353,26 +366,55 @@ next: for (int targetComponentIndex = 0; targetComponentIndex < infos.length; } /** - * Returns the matrix of an operation setting some coordinates to a constant values. + * Returns the matrix of an operation setting some coordinates to constant values. * * @param selected all operations from source to target {@link CompoundCRS}. - * @param n index of the first selected operation which describe a constant value. * @param srcDim number of dimensions in the target CRS of previous operation step. * @param tgtDim number of dimensions in the full (usually compound) target CRS. */ - static MatrixSIS createConstantOperation(final SubOperationInfo[] selected, int n, - final int srcDim, final int tgtDim) + static MatrixSIS createConstantOperation(final SubOperationInfo[] selected, final int srcDim, final int tgtDim) + throws TransformException { final boolean[] targetDimensionIsUsed = new boolean[tgtDim]; final MatrixSIS m = Matrices.createZero(tgtDim + 1, srcDim + 1); m.setElement(tgtDim, srcDim, 1); - do { - final SubOperationInfo component = selected[n]; - for (int j = component.startAtDimension; j < component.endAtDimension; j++) { - m.setElement(j, srcDim, component.constants[j]); - targetDimensionIsUsed[j] = true; + for (final SubOperationInfo component : selected) { + if (component.constants != null) { + /* + * Component for which no operation has been found. + * Set all coordinate values of that component to the specified constants + */ + for (int j = component.targetLowerDimension; j < component.targetUpperDimension; j++) { + m.setElement(j, srcDim, component.constants[j]); + targetDimensionIsUsed[j] = true; + } + } else { + /* + * Component for which an operation has been found, but maybe this is the inverse of the + * "Geographic 3D to 2D conversion" (EPSG:9659) which sets the height to `DEFAULT_HEIGHT`. + * If all values in the row are zero, then the coordinate value is unconditionally zero. + * Replace that value (usually the default ellipsoidal height) by the specified constant. + */ + final Matrix last = MathTransforms.getMatrix(MathTransforms.getLastStep(component.operation.getMathTransform())); + if (last != null) { +otherRow: for (int j = last.getNumRow() - 1; --j >= 0;) { // Ignore the last row. + for (int i = last.getNumCol(); --i >= 0;) { + if (last.getElement(j, i) != 0) { + continue otherRow; + } + } + final double[] constants = CoordinateOperationContext.getConstantCoordinates(); + if (constants == null) break; + final int k = component.targetLowerDimension + j; + m.setElement(j, srcDim, constants[k]); + targetDimensionIsUsed[k] = true; + } + } } - } while (++n < selected.length); + } + /* + * All coordinates that have not been set to a constant shall be propagated unchanged (scale factor of 1). + */ for (int i=0,j=0; j<tgtDim; j++) { if (!targetDimensionIsUsed[j]) { m.setElement(j, i++, 1); diff --git a/endorsed/src/org.apache.sis.referencing/main/org/apache/sis/referencing/operation/transform/MathTransforms.java b/endorsed/src/org.apache.sis.referencing/main/org/apache/sis/referencing/operation/transform/MathTransforms.java index 7be5444319..1c5c2429c1 100644 --- a/endorsed/src/org.apache.sis.referencing/main/org/apache/sis/referencing/operation/transform/MathTransforms.java +++ b/endorsed/src/org.apache.sis.referencing/main/org/apache/sis/referencing/operation/transform/MathTransforms.java @@ -635,8 +635,8 @@ public final class MathTransforms extends Static { * <li>Otherwise returns the given transform in a list of size 1.</li> * </ul> * - * @param transform the transform for which to get the components, or {@code null}. - * @return all single math transforms performed by this concatenated transform. + * @param transform the transform for which to get the steps, or {@code null}. + * @return all steps performed by the given transform. */ public static List<MathTransform> getSteps(final MathTransform transform) { if (transform != null) { @@ -650,6 +650,38 @@ public final class MathTransforms extends Static { } } + /** + * Returns the first step of the given (potentially concatenated) transform. + * Invoking this method is equivalent to invoking {@link #getSteps(MathTransform)} + * and retaining only the first element of the list. + * + * @param transform the transform for which to get the first step, or {@code null}. + * @return the first step performed by the given transform, or {@code null} if the argument was null. + * @since 1.5 + */ + public static MathTransform getFirstStep(MathTransform transform) { + while (transform instanceof ConcatenatedTransform) { + transform = ((ConcatenatedTransform) transform).transform1; + } + return transform; + } + + /** + * Returns the last step of the given (potentially concatenated) transform. + * Invoking this method is equivalent to invoking {@link #getSteps(MathTransform)} + * and retaining only the last element of the list. + * + * @param transform the transform for which to get the last step, or {@code null}. + * @return the last step performed by the given transform, or {@code null} if the argument was null. + * @since 1.5 + */ + public static MathTransform getLastStep(MathTransform transform) { + while (transform instanceof ConcatenatedTransform) { + transform = ((ConcatenatedTransform) transform).transform2; + } + return transform; + } + /** * Returns whether the given object is a linear transform. * This method is defined here for keeping it consistent with {@link #getMatrix(MathTransform)}.
