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 34a9ce8f2aa8fdcc9d92bd6aef233ebef25028f0
Author: Martin Desruisseaux <martin.desruisse...@geomatys.com>
AuthorDate: Sat Dec 29 17:01:35 2018 +0100
Target resolutions in GridGeometry.subgrid(Envelope, double...) should be
in same units than the given envelope.
---
.../org/apache/sis/coverage/grid/GridChange.java | 2 +-
.../org/apache/sis/coverage/grid/GridExtent.java | 3 +-
.../org/apache/sis/coverage/grid/GridGeometry.java | 107 ++---------
.../sis/coverage/grid/SubgridCalculator.java | 206 +++++++++++++++++++++
.../apache/sis/coverage/grid/GridGeometryTest.java | 20 ++
5 files changed, 245 insertions(+), 93 deletions(-)
diff --git
a/core/sis-raster/src/main/java/org/apache/sis/coverage/grid/GridChange.java
b/core/sis-raster/src/main/java/org/apache/sis/coverage/grid/GridChange.java
index 7c5faf1..f044f27 100644
--- a/core/sis-raster/src/main/java/org/apache/sis/coverage/grid/GridChange.java
+++ b/core/sis-raster/src/main/java/org/apache/sis/coverage/grid/GridChange.java
@@ -344,7 +344,7 @@ public class GridChange implements Serializable {
double[] factors = null;
Matrix toGiven = null;
if (strides != null) {
- // Validity of the strides values will be verified by GridExtent
constructor below.
+ // Validity of the strides values will be verified by GridExtent
constructor invoked below.
final GridExtent unscaled = extent;
final int dimension = extent.getDimension();
for (int i = Math.min(dimension, strides.length); --i >= 0;) {
diff --git
a/core/sis-raster/src/main/java/org/apache/sis/coverage/grid/GridExtent.java
b/core/sis-raster/src/main/java/org/apache/sis/coverage/grid/GridExtent.java
index 0e17166..28d1190 100644
--- a/core/sis-raster/src/main/java/org/apache/sis/coverage/grid/GridExtent.java
+++ b/core/sis-raster/src/main/java/org/apache/sis/coverage/grid/GridExtent.java
@@ -213,8 +213,7 @@ public class GridExtent implements Serializable {
/**
* Creates a new grid extent with low and high values adjusted by dividing
the {@linkplain #getSize(int) size}
* by the given strides. The policy of dividing the lower coordinates by
the stride shall be kept consistent
- * with {@link GridGeometry#subgrid(Envelope, double...)} and {@link
GridChange#getTargetGeometry(int...)}
- * computation of grid to CRS.
+ * with {@link GridChange#getTargetGeometry(int...)} computation of grid
to CRS.
*
* <div class="note"><b>Note:</b>
* if a division does not produce an integer, then the size is rounded
toward 0 (or toward negative infinity since
diff --git
a/core/sis-raster/src/main/java/org/apache/sis/coverage/grid/GridGeometry.java
b/core/sis-raster/src/main/java/org/apache/sis/coverage/grid/GridGeometry.java
index 2447b22..853c54e 100644
---
a/core/sis-raster/src/main/java/org/apache/sis/coverage/grid/GridGeometry.java
+++
b/core/sis-raster/src/main/java/org/apache/sis/coverage/grid/GridGeometry.java
@@ -30,11 +30,9 @@ import org.opengis.referencing.datum.PixelInCell;
import org.opengis.referencing.operation.Matrix;
import org.opengis.referencing.operation.MathTransform;
import org.opengis.referencing.operation.TransformException;
-import org.opengis.referencing.operation.CoordinateOperation;
import org.opengis.referencing.crs.CoordinateReferenceSystem;
import org.opengis.referencing.cs.CoordinateSystemAxis;
import org.opengis.referencing.cs.CoordinateSystem;
-import org.opengis.util.FactoryException;
import org.apache.sis.math.MathFunctions;
import org.apache.sis.geometry.Envelopes;
import org.apache.sis.geometry.GeneralEnvelope;
@@ -43,7 +41,6 @@ import org.apache.sis.referencing.IdentifiedObjects;
import org.apache.sis.referencing.operation.matrix.Matrices;
import org.apache.sis.referencing.operation.transform.MathTransforms;
import org.apache.sis.referencing.operation.transform.PassThroughTransform;
-import org.apache.sis.referencing.operation.transform.TransformSeparator;
import org.apache.sis.internal.referencing.DirectPositionView;
import org.apache.sis.internal.system.Modules;
import org.apache.sis.internal.raster.Resources;
@@ -55,7 +52,6 @@ import org.apache.sis.util.resources.Errors;
import org.apache.sis.util.logging.Logging;
import org.apache.sis.util.ArgumentChecks;
import org.apache.sis.util.CharSequences;
-import org.apache.sis.util.ArraysExt;
import org.apache.sis.util.Classes;
import org.apache.sis.util.Debug;
import org.apache.sis.io.TableAppender;
@@ -601,42 +597,10 @@ public class GridGeometry implements Serializable {
if (extent == null) {
throw incomplete(EXTENT, Resources.Keys.UnspecifiedGridExtent);
}
- MathTransform gridToAOI = cornerToCRS;
- if (gridToAOI == null) {
- throw incomplete(GRID_TO_CRS, Resources.Keys.UnspecifiedTransform);
+ if (cornerToCRS == null) {
+ throw incomplete(GridGeometry.GRID_TO_CRS,
Resources.Keys.UnspecifiedTransform);
}
- int[] modifiedDimensions = null;
- try {
- /*
- * If the envelope CRS is different than the expected CRS,
concatenate the envelope transformation
- * to the 'gridToCRS' transform. We should not transform the
envelope here - only concatenate the
- * transforms - because transforming envelopes twice add errors.
- */
- final CoordinateOperation operation =
Envelopes.findOperation(envelope, areaOfInterest);
- if (operation != null) {
- gridToAOI = MathTransforms.concatenate(gridToAOI,
operation.getMathTransform());
- }
- /*
- * If the envelope dimensions does not encompass all grid
dimensions, the envelope 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 modifiedDimensionCount = gridToAOI.getTargetDimensions();
- ArgumentChecks.ensureDimensionMatches("areaOfInterest",
modifiedDimensionCount, areaOfInterest);
- if (modifiedDimensionCount < gridToAOI.getSourceDimensions()) {
- final TransformSeparator sep = new
TransformSeparator(gridToAOI);
- sep.setTrimSourceDimensions(true);
- gridToAOI = sep.separate();
- modifiedDimensions = sep.getSourceDimensions();
- if (modifiedDimensions.length != modifiedDimensionCount) {
- throw new
TransformException(Resources.format(Resources.Keys.CanNotMapToGridDimensions));
- }
- }
- } catch (FactoryException e) {
- throw new
TransformException(Resources.format(Resources.Keys.CanNotMapToGridDimensions),
e);
- }
- final GridExtent sub = new
GridExtent(Envelopes.transform(gridToAOI.inverse(), areaOfInterest),
- GridRoundingMode.NEAREST, null, extent, modifiedDimensions);
- return sub.equals(extent) ? extent : sub;
+ return new SubgridCalculator(this, cornerToCRS, areaOfInterest,
null).extent;
}
/**
@@ -940,69 +904,32 @@ public class GridGeometry implements Serializable {
/**
* Returns a grid geometry over a sub-region of this grid geometry and
optionally with sub-sampling.
* The given envelope can be expressed in any coordinate reference system
(CRS) accepted by {@link #getExtent(Envelope)}.
- * The target resolution, if provided, shall be in the units of CRS axes,
in same order.
+ * The target resolution, if provided, shall be in same units than the
given envelope with axes in the same order.
+ * If the length of {@code targetResolution} array is less than the number
of dimensions of {@code areaOfInterest},
+ * then no sub-sampling will be applied on the missing dimensions.
*
* @param areaOfInterest the desired spatiotemporal region in any CRS
(transformations will be applied as needed),
* or {@code null} for not restricting the
sub-grid to a sub-area.
- * @param targetResolution the desired resolution in units of the CRS
axes, or {@code null} or an empty array
- * if no sub-sampling is desired.
- * @return a grid geometry over the specified sub-region of this grid
geometry.
+ * @param targetResolution the desired resolution in the same units and
order than the axes of the given envelope,
+ * or {@code null} or an empty array if no
sub-sampling is desired.
+ * @return a grid geometry over the specified sub-region of this grid
geometry with the specified resolution.
* @throws IncompleteGridGeometryException if this grid geometry has no
extent or no "grid to CRS" transform.
* @throws TransformException if an error occurred while converting the
envelope coordinates to grid coordinates.
*
* @see #getExtent(Envelope)
*/
public GridGeometry subgrid(final Envelope areaOfInterest, double...
targetResolution) throws TransformException {
- GridExtent domain = extent;
- if (areaOfInterest != null) {
- domain = getExtent(areaOfInterest);
+ if (extent == null) {
+ throw incomplete(EXTENT, Resources.Keys.UnspecifiedGridExtent);
}
- MathTransform mt = null;
- if (targetResolution != null && targetResolution.length != 0) {
- /*
- * Before to compute the strides, make sure we have all required
information.
- * One intent of following calls to getExtent() and
getGridToCRS(…) is to get
- * IncompleteGridGeometryException thrown if an information is
missing.
- */
- if (domain == null) domain = getExtent();
- final MathTransform gridToCRS =
getGridToCRS(PixelInCell.CELL_CENTER);
- /*
- * Convert the target resolutions to sub-samplings as number of
grid cells.
- * The sub-sampling will be scale factors for the new "grid to
CRS" transforms.
- */
- final int dimension = getDimension();
- targetResolution = ArraysExt.resize(targetResolution, dimension);
- Matrix m = gridToCRS.derivative(new
DirectPositionView.Double(domain.getPointOfInterest()));
- final double[] scales =
Matrices.inverse(m).multiply(targetResolution);
- /*
- * Creates the matrix to pre-concatenate to "grid to CRS"
transform. The (low /si) term in translation
- * below must be computed in the same way than the "low"
coordinates in GridExtent(GridExtent, int...)
- * constructor.
- */
- m = Matrices.createIdentity(dimension + 1);
- final int[] strides = new int[dimension];
- Arrays.fill(strides, 1);
- boolean isIdentity = true;
- for (int i=0; i<dimension; i++) {
- final double s = Math.floor(Math.abs(scales[i]));
- if (s > 1) { //
Also for skipping NaN values.
- final int si = (int) s;
- final long low = domain.getLow(i);
- m.setElement(i, i, s);
- m.setElement(i, dimension, low - (low / si) * s); //
(low / si) must be consistent with GridExtent.
- strides[i] = si;
- isIdentity = false;
- }
- }
- if (!isIdentity) {
- mt = MathTransforms.linear(m);
- domain = new GridExtent(domain, strides);
- }
+ if (cornerToCRS == null) {
+ throw incomplete(GRID_TO_CRS, Resources.Keys.UnspecifiedTransform);
}
- if (mt == null && Objects.equals(domain, extent)) {
- return this;
+ final SubgridCalculator sub = new SubgridCalculator(this, cornerToCRS,
areaOfInterest, targetResolution);
+ if (sub.toSubsampled != null || sub.extent != extent) {
+ return new GridGeometry(this, sub.extent, sub.toSubsampled);
}
- return new GridGeometry(this, domain, mt);
+ return this;
}
/**
diff --git
a/core/sis-raster/src/main/java/org/apache/sis/coverage/grid/SubgridCalculator.java
b/core/sis-raster/src/main/java/org/apache/sis/coverage/grid/SubgridCalculator.java
new file mode 100644
index 0000000..60ab1b3
--- /dev/null
+++
b/core/sis-raster/src/main/java/org/apache/sis/coverage/grid/SubgridCalculator.java
@@ -0,0 +1,206 @@
+/*
+ * 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.coverage.grid;
+
+import org.opengis.geometry.Envelope;
+import org.opengis.util.FactoryException;
+import org.opengis.referencing.operation.Matrix;
+import org.opengis.referencing.operation.MathTransform;
+import org.opengis.referencing.operation.TransformException;
+import org.opengis.referencing.operation.CoordinateOperation;
+import org.apache.sis.referencing.operation.transform.MathTransforms;
+import org.apache.sis.referencing.operation.transform.TransformSeparator;
+import org.apache.sis.referencing.operation.matrix.Matrices;
+import org.apache.sis.internal.referencing.DirectPositionView;
+import org.apache.sis.geometry.Envelopes;
+import org.apache.sis.geometry.GeneralEnvelope;
+import org.apache.sis.internal.raster.Resources;
+import org.apache.sis.util.ArgumentChecks;
+import org.apache.sis.util.ArraysExt;
+
+
+/**
+ * Helper class for computing the grid extent of a sub-area of a given grid
geometry.
+ * This class provides the {@link MathTransform} converting source grid
coordinates to target grid coordinates,
+ * where the source is the given {@link GridGeometry} instance and the target
is the sub-sampled grid.
+ *
+ * @author Martin Desruisseaux (Geomatys)
+ * @version 1.0
+ * @since 1.0
+ * @module
+ */
+final class SubgridCalculator {
+ /**
+ * The sub-extent computed by the constructor.
+ */
+ GridExtent extent;
+
+ /**
+ * The conversion from the original grid to the sub-sampled grid, or
{@code null} if no sub-sampling is applied.
+ * This is computed by the constructor.
+ */
+ MathTransform toSubsampled;
+
+ /**
+ * Computes the sub-grid over the given area of interest with the given
resolution.
+ * At least one of {@code areaOfInterest} and {@code resolution} shall be
non-null.
+ * It is caller's responsibility to ensure that {@link
GridGeometry#extent} is non-null.
+ *
+ * @param grid the enclosing grid geometry (mandatory).
+ * @param cornerToCRS the transform from cell corners to grid CRS
(mandatory).
+ * @param areaOfInterest the desired spatiotemporal region in any CRS,
or {@code null} for the whole area.
+ * @param resolution the desired resolution in the same units and
order than the axes of the AOI envelope,
+ * or {@code null} or an empty array if no
sub-sampling is desired.
+ * @throws TransformException if an error occurred while converting the
envelope coordinates to grid coordinates.
+ */
+ SubgridCalculator(final GridGeometry grid, MathTransform cornerToCRS,
final Envelope areaOfInterest, double[] resolution)
+ throws TransformException
+ {
+ /*
+ * List of grid dimensions that are modified by the 'cornerToCRS'
transform, or null for all dimensions.
+ * The length of this array is the number of dimensions of the given
Area Of Interest (AOI). Each value
+ * in this array is between 0 inclusive and 'extent.getDimension()'
exclusive.
+ */
+ int[] modifiedDimensions = null;
+ try {
+ /*
+ * If the envelope CRS is different than the expected CRS,
concatenate the envelope transformation
+ * to the 'gridToCRS' transform. We should not transform the
envelope here - only concatenate the
+ * transforms - because transforming envelopes twice would add
errors.
+ */
+ final CoordinateOperation operation =
Envelopes.findOperation(grid.envelope, areaOfInterest);
+ if (operation != null) {
+ cornerToCRS = MathTransforms.concatenate(cornerToCRS,
operation.getMathTransform());
+ }
+ /*
+ * If the envelope dimensions does not encompass all grid
dimensions, the envelope 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 = cornerToCRS.getTargetDimensions();
+ ArgumentChecks.ensureDimensionMatches("areaOfInterest", dimension,
areaOfInterest);
+ if (dimension < cornerToCRS.getSourceDimensions()) {
+ final TransformSeparator sep = new
TransformSeparator(cornerToCRS);
+ sep.setTrimSourceDimensions(true);
+ cornerToCRS = sep.separate();
+ modifiedDimensions = sep.getSourceDimensions();
+ if (modifiedDimensions.length != dimension) {
+ throw new
TransformException(Resources.format(Resources.Keys.CanNotMapToGridDimensions));
+ }
+ }
+ } catch (FactoryException e) {
+ throw new
TransformException(Resources.format(Resources.Keys.CanNotMapToGridDimensions),
e);
+ }
+ /*
+ * Compute the sub-extent for the given Area Of Interest (AOI),
ignoring for now the sub-sampling.
+ * If no area of interest has been specified, or if the result is
identical to the original extent,
+ * then we will keep the reference to the original GridExtent (i.e. we
share existing instances).
+ */
+ extent = grid.extent;
+ final int dimension = extent.getDimension();
+ GeneralEnvelope indices = null;
+ if (areaOfInterest != null) {
+ indices = Envelopes.transform(cornerToCRS.inverse(),
areaOfInterest);
+ setExtent(indices, extent, modifiedDimensions);
+ }
+ if (indices == null || indices.getDimension() != dimension) {
+ indices = new GeneralEnvelope(dimension);
+ }
+ for (int i=0; i<dimension; i++) {
+ indices.setRange(i, extent.getLow(i), extent.getHigh(i) + 1.0);
+ }
+ /*
+ * Convert the target resolutions to grid cell sub-samplings and
adjust the extent consequently.
+ * We perform this conversion by handling the resolution has a small
translation vector located
+ * at the point of interest, and converting it to a translation vector
in grid coordinates. The
+ * conversion is done by a multiplication with the "CRS to grid"
derivative at that point.
+ *
+ * The sub-sampling will be rounded in such a way that the difference
in grid size is less than
+ * one half of cell. Demonstration:
+ *
+ * e = Math.getExponent(span) → 2^e ≦ span
+ * a = e+1 → 2^a > span → 1/2^a <
1/span
+ * Δs = (s - round(s)) / 2^a
+ * (s - round(s)) ≦ 0.5 → Δs ≦ 0.5/2^a < 0.5/span
+ * Δs < 0.5/span → Δs⋅span < 0.5 cell.
+ */
+ if (resolution != null && resolution.length != 0) {
+ resolution = ArraysExt.resize(resolution,
cornerToCRS.getTargetDimensions());
+ Matrix m = cornerToCRS.derivative(new
DirectPositionView.Double(getPointOfInterest(modifiedDimensions)));
+ resolution = Matrices.inverse(m).multiply(resolution);
+ boolean modified = false;
+ for (int k=0; k<resolution.length; k++) {
+ double s = Math.abs(resolution[k]);
+ if (s > 1) { // Also for
skipping NaN values.
+ final int i = (modifiedDimensions != null) ?
modifiedDimensions[k] : k;
+ final int accuracy = Math.max(0,
Math.getExponent(indices.getSpan(i))) + 1; // Power of 2.
+ s = Math.scalb(Math.rint(Math.scalb(s, accuracy)),
-accuracy);
+ indices.setRange(i, indices.getLower(i) / s,
indices.getUpper(i) / s);
+ modified = true;
+ }
+ resolution[k] = s;
+ }
+ /*
+ * If at least one sub-sampling is effective, build a scale from
the old grid coordinates to the new
+ * grid coordinates. If we had no rounding, the conversion would
be only a scale. But because of rounding,
+ * we need a small translation for the difference between the
"real" coordinate and the integer coordinate.
+ */
+ if (modified) {
+ final GridExtent unscaled = extent;
+ setExtent(indices, null, null);
+ m = Matrices.createIdentity(dimension + 1);
+ for (int k=0; k<resolution.length; k++) {
+ final double s = resolution[k];
+ if (s > 1) { // Also for
skipping NaN values.
+ final int i = (modifiedDimensions != null) ?
modifiedDimensions[k] : k;
+ m.setElement(i, i, s);
+ m.setElement(i, dimension, unscaled.getLow(i) -
extent.getLow(i) * s);
+ }
+ }
+ toSubsampled = MathTransforms.linear(m);
+ }
+ }
+ }
+
+ /**
+ * Sets {@link #extent} to the given envelope, rounded to nearest integers.
+ *
+ * @param indices the envelope to use for setting the grid
extent.
+ * @param enclosing the enclosing grid extent if a sub-sampling
is not yet applied, {@code null} otherwise.
+ * @param modifiedDimensions if {@code enclosing} is non-null, the grid
dimensions to set from the envelope.
+ */
+ private void setExtent(final GeneralEnvelope indices, final GridExtent
enclosing, final int[] modifiedDimensions) {
+ final GridExtent sub = new GridExtent(indices,
GridRoundingMode.NEAREST, null, enclosing, modifiedDimensions);
+ if (!sub.equals(extent)) {
+ extent = sub;
+ }
+ }
+
+ /**
+ * Returns the point of interest of current {@link #extent}.
+ */
+ private double[] getPointOfInterest(final int[] modifiedDimensions) {
+ final double[] pointOfInterest = extent.getPointOfInterest();
+ if (modifiedDimensions == null) {
+ return pointOfInterest;
+ }
+ final double[] filtered = new double[modifiedDimensions.length];
+ for (int i=0; i<filtered.length; i++) {
+ filtered[i] = pointOfInterest[modifiedDimensions[i]];
+ }
+ return filtered;
+ }
+}
diff --git
a/core/sis-raster/src/test/java/org/apache/sis/coverage/grid/GridGeometryTest.java
b/core/sis-raster/src/test/java/org/apache/sis/coverage/grid/GridGeometryTest.java
index bb651e7..fd592b6 100644
---
a/core/sis-raster/src/test/java/org/apache/sis/coverage/grid/GridGeometryTest.java
+++
b/core/sis-raster/src/test/java/org/apache/sis/coverage/grid/GridGeometryTest.java
@@ -25,6 +25,7 @@ import org.apache.sis.referencing.operation.matrix.Matrix4;
import org.apache.sis.referencing.operation.matrix.Matrices;
import org.apache.sis.referencing.operation.transform.MathTransforms;
import org.apache.sis.referencing.crs.HardCodedCRS;
+import org.apache.sis.geometry.DirectPosition2D;
import org.apache.sis.geometry.GeneralEnvelope;
import org.apache.sis.test.DependsOnMethod;
import org.apache.sis.test.DependsOn;
@@ -347,5 +348,24 @@ public final strictfp class GridGeometryTest extends
TestCase {
0, 1, -90,
2, 0, -180,
0, 0, 1),
MathTransforms.getMatrix(grid.getGridToCRS(PixelInCell.CELL_CORNER)), STRICT);
+ /*
+ * A sub-region again but with a requested resolution which is not a
divisor of the actual resolution.
+ * It will force GridGeometry to adjust the translation term to
compensate. We verify that the adustment
+ * is correct by verifying that we still get the same envelope.
+ */
+ grid = grid.subgrid(envelope, 3, 2);
+ assertExtentEquals(new long[] {94, 13}, new long[] {95, 39},
grid.getExtent());
+ assertEnvelopeEquals(envelope, grid.getEnvelope(), STRICT);
+ MathTransform cornerToCRS = grid.getGridToCRS(PixelInCell.CELL_CORNER);
+ assertMatrixEquals("gridToCRS", new Matrix3(
+ 0, 3, -89,
+ 2, 0, -180,
+ 0, 0, 1), MathTransforms.getMatrix(cornerToCRS), STRICT);
+
+ DirectPosition2D src = new DirectPosition2D();
+ DirectPosition2D tgt = new DirectPosition2D();
+ DirectPosition2D exp = new DirectPosition2D();
+ src.x = 94; src.y = 13; exp.x = -50; exp.y = 8; assertEquals("Lower
corner", exp, cornerToCRS.transform(src, tgt));
+ src.x = 96; src.y = 40; exp.x = +31; exp.y = 12; assertEquals("Upper
corner", exp, cornerToCRS.transform(src, tgt));
}
}
