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 a95fbd4ab4177bcf1e5d8d46c1ce4c688758ab2f
Author: Martin Desruisseaux <martin.desruisse...@geomatys.com>
AuthorDate: Thu Dec 19 13:18:43 2024 +0100
Allow GeoTIFF writer to encode "grid to CRS" where the vertical component
is non-linear.
---
.../org/apache/sis/storage/geotiff/Writer.java | 10 +-
.../sis/storage/geotiff/writer/GeoEncoder.java | 144 +++++++++++++--------
.../sis/storage/geotiff/GeoTiffStoreTest.java | 111 ++++++++++++++++
3 files changed, 211 insertions(+), 54 deletions(-)
diff --git
a/endorsed/src/org.apache.sis.storage.geotiff/main/org/apache/sis/storage/geotiff/Writer.java
b/endorsed/src/org.apache.sis.storage.geotiff/main/org/apache/sis/storage/geotiff/Writer.java
index 1ce0d8140c..6801c886f5 100644
---
a/endorsed/src/org.apache.sis.storage.geotiff/main/org/apache/sis/storage/geotiff/Writer.java
+++
b/endorsed/src/org.apache.sis.storage.geotiff/main/org/apache/sis/storage/geotiff/Writer.java
@@ -39,11 +39,14 @@ import javax.measure.IncommensurableException;
import org.opengis.util.FactoryException;
import org.opengis.metadata.Metadata;
import org.opengis.metadata.citation.CitationDate;
+import org.opengis.referencing.operation.TransformException;
import org.apache.sis.image.ImageProcessor;
import org.apache.sis.coverage.grid.GridGeometry;
+import org.apache.sis.coverage.grid.IncompleteGridGeometryException;
import org.apache.sis.coverage.privy.ImageUtilities;
import org.apache.sis.storage.DataStoreException;
import org.apache.sis.storage.DataStoreReferencingException;
+import org.apache.sis.storage.IncompatibleResourceException;
import org.apache.sis.storage.ReadOnlyStorageException;
import org.apache.sis.storage.base.MetadataFetcher;
import org.apache.sis.storage.geotiff.writer.TagValue;
@@ -58,6 +61,9 @@ import org.apache.sis.util.privy.Numerics;
import org.apache.sis.util.resources.Errors;
import org.apache.sis.math.Fraction;
+// Specific to the geoapi-3.1 and geoapi-4.0 branches:
+import org.opengis.coverage.CannotEvaluateException;
+
/**
* An image writer for GeoTIFF files. This writer duplicates the
implementations performed by other libraries,
@@ -403,9 +409,11 @@ final class Writer extends IOBase implements Flushable {
};
mf.accept(metadata);
GeoEncoder geoKeys = null;
- if (grid != null && grid.isDefined(GridGeometry.GRID_TO_CRS)) try {
+ if (grid != null) try {
geoKeys = new GeoEncoder(store.listeners());
geoKeys.write(grid, mf);
+ } catch (IncompleteGridGeometryException | CannotEvaluateException |
TransformException e) {
+ throw new IncompatibleResourceException(e.getMessage(),
e).addAspect("gridGeometry");
} catch (FactoryException | IncommensurableException |
RuntimeException e) {
throw new DataStoreReferencingException(e.getMessage(), e);
}
diff --git
a/endorsed/src/org.apache.sis.storage.geotiff/main/org/apache/sis/storage/geotiff/writer/GeoEncoder.java
b/endorsed/src/org.apache.sis.storage.geotiff/main/org/apache/sis/storage/geotiff/writer/GeoEncoder.java
index 25536bca7b..1ce9787fc4 100644
---
a/endorsed/src/org.apache.sis.storage.geotiff/main/org/apache/sis/storage/geotiff/writer/GeoEncoder.java
+++
b/endorsed/src/org.apache.sis.storage.geotiff/main/org/apache/sis/storage/geotiff/writer/GeoEncoder.java
@@ -16,6 +16,7 @@
*/
package org.apache.sis.storage.geotiff.writer;
+import java.util.Arrays;
import java.util.List;
import java.util.EnumMap;
import java.util.logging.Level;
@@ -45,6 +46,7 @@ import org.opengis.referencing.datum.GeodeticDatum;
import org.opengis.referencing.datum.VerticalDatum;
import org.opengis.referencing.operation.Conversion;
import org.opengis.referencing.operation.Matrix;
+import org.opengis.referencing.operation.TransformException;
import org.opengis.parameter.GeneralParameterValue;
import org.opengis.parameter.ParameterValue;
import org.apache.sis.measure.Units;
@@ -59,6 +61,7 @@ import org.apache.sis.referencing.cs.CoordinateSystems;
import org.apache.sis.referencing.operation.matrix.Matrices;
import org.apache.sis.referencing.operation.transform.MathTransforms;
import org.apache.sis.referencing.factory.UnavailableFactoryException;
+import org.apache.sis.referencing.privy.AxisDirections;
import org.apache.sis.referencing.privy.ReferencingUtilities;
import org.apache.sis.referencing.privy.WKTKeywords;
import org.apache.sis.coverage.grid.GridGeometry;
@@ -84,6 +87,12 @@ import org.apache.sis.pending.jdk.JDK15;
* @author Martin Desruisseaux (Geomatys)
*/
public final class GeoEncoder {
+ /**
+ * Number of dimensions in a rendered image.
+ * Used for identifying codes where a two-dimensional slice is assumed.
+ */
+ private static final int BIDIMENSIONAL = 2;
+
/**
* Size of the model transformation matrix, in number of rows and columns.
* This size is fixed by the GeoTIFF specification.
@@ -102,16 +111,11 @@ public final class GeoEncoder {
private String citation;
/**
- * The coordinate reference system of the grid geometry, or {@code null}
if none.
- * This CRS may contain more dimensions than the 3 dimensions allowed by
GeoTIFF.
+ * The axis directions of the grid geometry, or {@code null} if none.
+ * The array length is 2 or 3, with the vertical axis always last.
* Axis order and axis directions may be different than the (east, north,
up) directions mandated by GeoTIFF.
*/
- private CoordinateReferenceSystem fullCRS;
-
- /**
- * Whether the coordinate system has a vertical component.
- */
- private boolean hasVerticalAxis;
+ private AxisDirection[] axisDirections;
/**
* The conversion from grid coordinates to full CRS, which determines the
model transformation.
@@ -228,31 +232,72 @@ public final class GeoEncoder {
* @throws IncompatibleResourceException if the grid geometry cannot be
encoded.
*/
public void write(final GridGeometry grid, final MetadataFetcher<?>
metadata)
- throws FactoryException, IncommensurableException,
IncompatibleResourceException
+ throws FactoryException, TransformException,
IncommensurableException, IncompatibleResourceException
{
- citation = CollectionsExt.first(metadata.transformationDimension);
- isPoint = CollectionsExt.first(metadata.cellGeometry) ==
CellGeometry.POINT;
- gridToCRS = MathTransforms.getMatrix(grid.getGridToCRS(isPoint ?
PixelInCell.CELL_CENTER : PixelInCell.CELL_CORNER));
- if (gridToCRS == null) {
- String message =
resources().getString(Resources.Keys.CanNotEncodeNonLinearModel);
- throw new
IncompatibleResourceException(message).addAspect("gridToCRS");
- }
- if (grid.isDefined(GridGeometry.CRS)) {
- fullCRS = grid.getCoordinateReferenceSystem();
- final CoordinateReferenceSystem crs =
CRS.getHorizontalComponent(fullCRS);
- if (crs instanceof ProjectedCRS) {
- writeCRS((ProjectedCRS) crs);
- writeCRS(CRS.getVerticalComponent(fullCRS, true));
- } else if (crs instanceof GeodeticCRS) {
- writeCRS((GeodeticCRS) crs, false);
- writeCRS(CRS.getVerticalComponent(fullCRS, true));
- } else if (fullCRS instanceof EngineeringCRS &&
ReferencingUtilities.getDimension(fullCRS) == 2) {
- writeModelType(GeoCodes.userDefined);
- } else {
- throw unsupportedType(fullCRS);
+ citation = CollectionsExt.first(metadata.transformationDimension);
+ isPoint = CollectionsExt.first(metadata.cellGeometry) ==
CellGeometry.POINT;
+ final var anchor = isPoint ? PixelInCell.CELL_CENTER :
PixelInCell.CELL_CORNER;
+ /*
+ * Get the dimension indices of the two-dimensional slice to write.
They should be the first dimensions,
+ * but we allow those dimensions to appear elsewhere.
+ */
+ final int[] dimensions =
grid.getExtent().getSubspaceDimensions(BIDIMENSIONAL);
+ final GridGeometry horizontal = grid.selectDimensions(dimensions);
+ if (grid.isDefined(GridGeometry.GRID_TO_CRS)) {
+ gridToCRS =
MathTransforms.getMatrix(horizontal.getGridToCRS(anchor));
+ if (gridToCRS == null) {
+ String message =
resources().getString(Resources.Keys.CanNotEncodeNonLinearModel);
+ throw new
IncompatibleResourceException(message).addAspect("gridToCRS");
}
- } else {
+ }
+ /*
+ * Write the horizontal component of the CRS. We need to take the CRS
+ * at the same dimensions than the ones selected for the `gridToCRS`.
+ */
+ if (!grid.isDefined(GridGeometry.CRS)) {
writeModelType(GeoCodes.undefined);
+ return;
+ }
+ final CoordinateReferenceSystem crs =
horizontal.getCoordinateReferenceSystem();
+ axisDirections =
CoordinateSystems.getAxisDirections(crs.getCoordinateSystem());
+ if (crs instanceof ProjectedCRS) {
+ writeCRS((ProjectedCRS) crs);
+ } else if (crs instanceof GeodeticCRS) {
+ writeCRS((GeodeticCRS) crs, false);
+ } else if (crs instanceof EngineeringCRS) {
+ writeModelType(GeoCodes.userDefined);
+ return;
+ } else {
+ throw unsupportedType(crs);
+ }
+ /*
+ * Write the vertical component of the CRS, if any. We restrict the
type to `VerticalCRS`,
+ * because this is the semantic of the GeoKeys.
+ */
+ final CoordinateReferenceSystem fullCRS =
grid.getCoordinateReferenceSystem();
+ final VerticalCRS vertical = CRS.getVerticalComponent(fullCRS, true);
+ if (vertical != null) {
+ final CoordinateSystem cs = vertical.getCoordinateSystem();
+ final int vi =
AxisDirections.indexOfColinear(fullCRS.getCoordinateSystem(), cs);
+ if (vi >= 0 && Arrays.binarySearch(dimensions, vi) < 0) {
+ writeCRS(vertical);
+ axisDirections = Arrays.copyOf(axisDirections,
BIDIMENSIONAL+1);
+ axisDirections[BIDIMENSIONAL] = cs.getAxis(0).getDirection();
+ if (gridToCRS != null) {
+ gridToCRS = Matrices.resizeAffine(gridToCRS, MATRIX_SIZE,
MATRIX_SIZE);
+ Matrix more = grid.getLinearGridToCRS(anchor).getMatrix();
+ for (int i=0; i<MATRIX_SIZE; i++) {
+ final int s;
+ switch (i) {
+ default: s = dimensions[i];
break; // Shear from horizontal dimensions.
+ case BIDIMENSIONAL: s = vi;
break; // Scale from vertical dimension.
+ case MATRIX_SIZE-1: s = more.getNumCol() - 1;
break; // Translation.
+ }
+ // Copy the rows of the third dimension.
+ gridToCRS.setElement(BIDIMENSIONAL, i,
more.getElement(BIDIMENSIONAL, s));
+ }
+ }
+ }
}
}
@@ -280,22 +325,19 @@ public final class GeoEncoder {
* @throws IncompatibleResourceException if a unit of measurement cannot
be encoded.
*/
private void writeCRS(final VerticalCRS crs) throws FactoryException,
IncompatibleResourceException {
- if (crs != null) {
- hasVerticalAxis = true;
- if (writeEPSG(GeoKeys.Vertical, crs)) {
- writeName(GeoKeys.VerticalCitation, null, crs);
- addUnits(UnitKey.VERTICAL, crs.getCoordinateSystem());
- final VerticalDatum datum = PseudoDatum.of(crs);
- if (writeEPSG(GeoKeys.VerticalDatum, datum)) {
- /*
- * OGC requirement 25.5 said "VerticalCitationGeoKey SHALL
be populated."
- * But how? Using the same multiple-names convention as
for geodetic CRS?
- *
- * https://github.com/opengeospatial/geotiff/issues/59
- */
- }
- writeUnit(UnitKey.VERTICAL);
+ if (writeEPSG(GeoKeys.Vertical, crs)) {
+ writeName(GeoKeys.VerticalCitation, null, crs);
+ addUnits(UnitKey.VERTICAL, crs.getCoordinateSystem());
+ final VerticalDatum datum = PseudoDatum.of(crs);
+ if (writeEPSG(GeoKeys.VerticalDatum, datum)) {
+ /*
+ * OGC requirement 25.5 said "VerticalCitationGeoKey SHALL be
populated."
+ * But how? Using the same multiple-names convention as for
geodetic CRS?
+ *
+ * https://github.com/opengeospatial/geotiff/issues/59
+ */
}
+ writeUnit(UnitKey.VERTICAL);
}
}
@@ -716,16 +758,12 @@ public final class GeoEncoder {
* If the CRS of the grid geometry has different axis order, we need
to adjust the
* "grid to CRS" transform.
*/
- if (fullCRS != null) {
- final AxisDirection[] source =
CoordinateSystems.getAxisDirections(fullCRS.getCoordinateSystem());
- final AxisDirection[] target = new AxisDirection[hasVerticalAxis ?
3 : 2];
+ if (axisDirections != null) {
+ final var target = axisDirections.clone(); // Preserve vertical
axis direction.
target[0] = AxisDirection.EAST;
target[1] = AxisDirection.NORTH;
- if (hasVerticalAxis) {
- target[2] = AxisDirection.UP;
- }
- gridToCRS = Matrices.createTransform(source,
target).multiply(gridToCRS);
- fullCRS = null; // For avoiding to do the multiplication
again.
+ gridToCRS = Matrices.createTransform(axisDirections,
target).multiply(gridToCRS);
+ axisDirections = null; // For avoiding to do the
multiplication again.
}
/*
* Copy matrix coefficients. This matrix size is always 4×4, no matter
the size of the `gridToCRS` matrix.
diff --git
a/endorsed/src/org.apache.sis.storage.geotiff/test/org/apache/sis/storage/geotiff/GeoTiffStoreTest.java
b/endorsed/src/org.apache.sis.storage.geotiff/test/org/apache/sis/storage/geotiff/GeoTiffStoreTest.java
new file mode 100644
index 0000000000..feb90f7f13
--- /dev/null
+++
b/endorsed/src/org.apache.sis.storage.geotiff/test/org/apache/sis/storage/geotiff/GeoTiffStoreTest.java
@@ -0,0 +1,111 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements. See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License. You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.apache.sis.storage.geotiff;
+
+import java.nio.file.Path;
+import java.nio.file.Files;
+import java.awt.image.BufferedImage;
+import java.awt.image.RenderedImage;
+import org.opengis.referencing.cs.AxisDirection;
+import org.apache.sis.storage.DataStore;
+import org.apache.sis.storage.DataStores;
+import org.apache.sis.storage.StorageConnector;
+import org.apache.sis.storage.GridCoverageResource;
+import org.apache.sis.coverage.grid.GridExtent;
+import org.apache.sis.coverage.grid.GridGeometry;
+import org.apache.sis.coverage.grid.PixelInCell;
+import org.apache.sis.coverage.grid.GridCoverage;
+import org.apache.sis.coverage.grid.GridCoverageBuilder;
+import org.apache.sis.referencing.operation.transform.MathTransforms;
+import org.apache.sis.referencing.operation.matrix.Matrix4;
+import org.apache.sis.referencing.CRS;
+
+// Test dependencies
+import org.junit.jupiter.api.Test;
+import static org.junit.jupiter.api.Assertions.*;
+import org.apache.sis.test.TestCase;
+import org.apache.sis.test.TestUtilities;
+import org.apache.sis.referencing.crs.HardCodedCRS;
+
+// Specific to the geoapi-3.1 and geoapi-4.0 branches:
+import static org.opengis.test.Assertions.assertMatrixEquals;
+import static org.opengis.test.Assertions.assertAxisDirectionsEqual;
+
+
+/**
+ * integration tests for {@link GeoTiffStore}.
+ *
+ * @author Martin Desruisseaux (Geomatys)
+ */
+public final class GeoTiffStoreTest extends TestCase {
+ /**
+ * Creates a new test case.
+ */
+ public GeoTiffStoreTest() {
+ }
+
+ /**
+ * Tests writing an image with a non-linear vertical component in the
"grid to CRS" transform.
+ * This method merely tests that no exception is thrown during the
execution, and that reading
+ * the image back can give back the three-dimensional <abbr>CRS</abbr>.
+ *
+ * @throws Exception if an error occurred while preparing or running the
test.
+ */
+ @Test
+ public void testNonLinearVerticalTransform() throws Exception {
+ final int width = 5;
+ final int height = 4;
+ final var builder = new GridCoverageBuilder();
+ final var extent = new GridExtent(null, null, new long[] {width,
height, 1}, false);
+ final var crs = CRS.compound(HardCodedCRS.WGS84_LATITUDE_FIRST,
HardCodedCRS.DEPTH);
+ final var gridToCRS = MathTransforms.compound(
+ MathTransforms.scale(0.25, -0.25),
+ MathTransforms.interpolate(new double[] {0, 1, 2}, new
double[] {3, 4, 9}));
+ builder.setDomain(new GridGeometry(extent, PixelInCell.CELL_CORNER,
gridToCRS, crs));
+ builder.setValues(new BufferedImage(width, height,
BufferedImage.TYPE_INT_ARGB));
+ final GridCoverage coverage = builder.build();
+ final Path file = Files.createTempFile("sis-test-", ".tiff");
+ try {
+ try (DataStore store = DataStores.openWritable(new
StorageConnector(file), "GeoTIFF")) {
+ assertInstanceOf(GeoTiffStore.class, store).append(coverage,
null);
+ }
+ /*
+ * Read the image that we wrote in above block. This block merely
tests that no exception is thrown,
+ * and that the result has the expected number of dimensions, axis
order and scale factors.
+ */
+ try (DataStore store = DataStores.open(new StorageConnector(file),
"GeoTIFF")) {
+ GridCoverageResource r =
TestUtilities.getSingleton(assertInstanceOf(GeoTiffStore.class,
store).components());
+ GridGeometry gg = r.getGridGeometry();
+ assertEquals(3, gg.getDimension());
+
assertAxisDirectionsEqual(gg.getCoordinateReferenceSystem().getCoordinateSystem(),
+ AxisDirection.EAST, AxisDirection.NORTH,
AxisDirection.DOWN);
+
+ assertMatrixEquals(new Matrix4(0, -0.25, 0, 0,
+ 0.25, 0, 0, 0,
+ 0, 0, 1, 3,
+ 0, 0, 0, 1),
+
MathTransforms.getMatrix(gg.getGridToCRS(PixelInCell.CELL_CORNER)), 0,
"gridToCRS");
+
+ RenderedImage image = r.read(null).render(null);
+ assertEquals(width, image.getWidth(), "width");
+ assertEquals(height, image.getHeight(), "height");
+ }
+ } finally {
+ Files.delete(file);
+ }
+ }
+}
