in endorsed/src/org.apache.sis.feature/main/org/apache/sis/coverage/grid/ResampledGridCoverage.java [296:481]
static GridCoverage create(final GridCoverage source, final GridGeometry target, final ImageProcessor processor,
final boolean allowOperationReplacement)
throws FactoryException, TransformException
{
final GridGeometry sourceGG = source.getGridGeometry();
final CoordinateOperationFinder changeOfCRS = new CoordinateOperationFinder(sourceGG, target);
changeOfCRS.verifyPresenceOfCRS(true);
/*
* Compute the transform from source pixels to target CRS (to be completed to target pixels later).
* The following lines may throw IncompleteGridGeometryException, which is desired because if that
* transform is missing, we cannot continue (we have no way to guess it).
*/
// Finder is initialized to PixelInCell.CELL_CORNER.
final MathTransform sourceCornerToCRS = changeOfCRS.gridToCRS();
final MathTransform crsToSourceCorner = changeOfCRS.inverse();
changeOfCRS.setAnchor(PixelInCell.CELL_CENTER);
final MathTransform sourceCenterToCRS = changeOfCRS.gridToCRS();
final MathTransform crsToSourceCenter = changeOfCRS.inverse();
/*
* Compute the transform from target grid to target CRS. This transform may be unspecified,
* in which case we need to compute a default transform trying to preserve resolution at the
* point of interest.
*/
boolean isGeometryExplicit = target.isDefined(GridGeometry.EXTENT);
GridExtent targetExtent = isGeometryExplicit ? target.getExtent() : null;
final MathTransform targetCenterToCRS;
if (target.isDefined(GridGeometry.GRID_TO_CRS)) {
isGeometryExplicit = true;
targetCenterToCRS = target.getGridToCRS(PixelInCell.CELL_CENTER);
if (targetExtent == null) {
targetExtent = targetExtent(sourceGG.getExtent(), sourceCornerToCRS,
target.getGridToCRS(PixelInCell.CELL_CORNER).inverse(), false);
}
} else {
/*
* We will try to preserve resolution at the point of interest, which is typically in the center.
* We will also try to align the grids in such a way that integer coordinates close to the point
* of interest are integers in both grids. This correction is given by the fractional digits of
* `originToPOI` vector (the integer digits do not matter; they will be cancelled later).
*/
final GridExtent sourceExtent = sourceGG.getExtent();
final double[] sourcePOI = sourceExtent.getPointOfInterest(PixelInCell.CELL_CENTER);
final double[] targetPOI = new double[sourceCenterToCRS.getTargetDimensions()];
final MatrixSIS vectors = MatrixSIS.castOrCopy(MathTransforms.derivativeAndTransform(
sourceCenterToCRS, sourcePOI, 0, targetPOI, 0));
final double[] originToPOI = vectors.multiply(sourcePOI);
/*
* The first `vectors` column gives the displacement in target CRS when moving in source grid by one cell
* toward right, and the second column gives the displacement when moving one cell toward up (positive y).
* More columns may exist in 3D, 4D, etc. cases. We retain only the magnitudes of those vectors, in order
* to build new vectors with directions parallel with target grid axes. There is one magnitude value for
* each target CRS dimension. If there is more target grid dimensions than the number of magnitude values
* (unusual, but not forbidden), some grid dimensions will be ignored provided that their size is 1
* (otherwise a SubspaceNotSpecifiedException is thrown).
*/
final MatrixSIS magnitudes = vectors.normalizeColumns(); // Length in dimension of source grid.
final int crsDim = vectors.getNumRow(); // Number of dimensions of target CRS.
final int gridDim = target.getDimension(); // Number of dimensions of target grid.
final int mappedDim = Math.min(magnitudes.getNumCol(), Math.min(crsDim, gridDim));
final MatrixSIS crsToGrid = Matrices.create(gridDim + 1, crsDim + 1, ExtendedPrecisionMatrix.CREATE_ZERO);
final int[] dimSelect = (gridDim > crsDim && targetExtent != null) ?
targetExtent.getSubspaceDimensions(crsDim) : null;
/*
* The goal below is to build a target "gridToCRS" which perform the same axis swapping and flipping than
* the source "gridToCRS". For example if the source "gridToCRS" was flipping y axis, then we want target
* "gridToCRS" to also flips that axis, unless the transformation from "source CRS" to "target CRS" flips
* that axis, in which case the result in target "gridToCRS" would be to not flip again:
*
* (source gridToCRS) → (source CRS to target CRS) → (target gridToCRS)⁻¹
* flip y axis no axis direction change flip y axis
* or flip y axis flip y axis no additional flip
*
* For each column, the row index of the greatest absolute value is taken as the target dimension where
* to set vector magnitude in target "crsToGrid" matrix. That way, if transformation from source CRS to
* target CRS does not flip or swap axis, target `gridToCRS` matrix should looks like source `gridToCRS`
* matrix (i.e. zero and non-zero coefficients at the same places). If two vectors have their greatest
* value on the same row, the largest value win (since the `vectors` matrix has been normalized to unit
* vectors, values in different columns are comparable).
*/
for (;;) {
double max = -1;
int sign = 1;
int tgDim = -1; // Grid dimension of maximal value.
int tcDim = -1; // CRS dimension of maximal value.
for (int i=0; i<mappedDim; i++) {
// `ci` differs from `i` only if the source grid has "too much" dimensions.
final int ci = (dimSelect != null) ? dimSelect[i] : i;
for (int j=0; j<crsDim; j++) {
final double v = vectors.getElement(j, ci);
final double m = Math.abs(v);
if (m > max) {
max = m;
sign = (v < 0) ? -1 : 1; // Like `Math.signum(…)` but without 0.
tcDim = j;
tgDim = ci;
}
}
}
if (tgDim < 0) break; // No more non-zero value in the `vectors` matrix.
for (int j=0; j<crsDim; j++) {
vectors.setElement(j, tgDim, Double.NaN); // For preventing this column to be selected again.
}
for (int i=0; i<gridDim; i++) {
vectors.setElement(tcDim, i, Double.NaN); // For preventing this row to be selected again.
}
DoubleDouble m = DoubleDouble.of(sign);
m = m.divide(magnitudes.getNumber(0, tgDim), false);
crsToGrid.setNumber(tgDim, tcDim, m); // Scale factor from CRS coordinates to grid coordinates.
/*
* Move the point of interest in a place where conversion to source grid coordinates
* will be close to integer. The exact location does not matter; an additional shift
* will be applied later for translating to target grid extent.
*/
m = m.multiply(DoubleDouble.sum(originToPOI[tcDim], -targetPOI[tcDim]));
crsToGrid.setNumber(tgDim, crsDim, m);
}
crsToGrid.setElement(gridDim, crsDim, 1);
/*
* At this point we got a first estimation of "target CRS to grid" transform, without translation terms.
* Apply the complete transform chain on source extent; this will give us a tentative target extent.
* This tentative extent will be compared with desired target.
*/
final GridExtent tentative = targetExtent(sourceExtent, sourceCornerToCRS,
MathTransforms.linear(crsToGrid), true);
if (targetExtent == null) {
// Create an extent of same size but with lower coordinates set to 0.
if (tentative.startsAtZero()) {
targetExtent = tentative;
} else {
final long[] coordinates = new long[gridDim * 2];
for (int i=0; i<gridDim; i++) {
coordinates[i + gridDim] = tentative.getSize(i) - 1;
}
targetExtent = new GridExtent(tentative, coordinates);
}
}
/*
* At this point we have the desired target extent and the extent that we actually got by applying
* full "source to target" transform. Compute the scale and offset differences between target and
* actual extents, then adjust matrix coefficients for compensating those differences.
*/
for (int j=0; j<gridDim; j++) {
DoubleDouble span = DoubleDouble.of(targetExtent.getSize(j));
DoubleDouble scale = DoubleDouble.of(tentative.getSize(j));
scale = span.divide(scale);
DoubleDouble offset = DoubleDouble.of(targetExtent.getLow(j));
offset = offset.subtract(scale.multiply(tentative.getLow(j)));
crsToGrid.convertAfter(j, scale, offset);
}
targetCenterToCRS = MathTransforms.linear(crsToGrid.inverse());
}
/*
* At this point all target grid geometry components are non-null.
* Build the final target GridGeometry if any components were missing.
* If an envelope is defined, resample only that sub-region.
*/
GridGeometry complete = target;
ComparisonMode mode = ComparisonMode.IGNORE_METADATA;
if (!target.isDefined(GridGeometry.EXTENT | GridGeometry.GRID_TO_CRS | GridGeometry.CRS)) {
final CoordinateReferenceSystem targetCRS = changeOfCRS.getTargetCRS();
complete = new GridGeometry(targetExtent, PixelInCell.CELL_CENTER, targetCenterToCRS, targetCRS);
mode = ComparisonMode.APPROXIMATE;
if (target.isDefined(GridGeometry.ENVELOPE)) {
final MathTransform targetCornerToCRS = complete.getGridToCRS(PixelInCell.CELL_CORNER);
GeneralEnvelope bounds = new GeneralEnvelope(complete.getEnvelope());
bounds.intersect(target.getEnvelope());
bounds = Envelopes.transform(targetCornerToCRS.inverse(), bounds);
targetExtent = new GridExtent(bounds, false, GridRoundingMode.NEAREST, GridClippingMode.STRICT, null, null, targetExtent, null);
complete = new GridGeometry(targetExtent, PixelInCell.CELL_CENTER, targetCenterToCRS, targetCRS);
isGeometryExplicit = true;
}
}
if (sourceGG.equals(complete, mode)) {
return source;
}
/*
* Complete the "target to source" transform.
*/
final MathTransform targetCornerToCRS = complete.getGridToCRS(PixelInCell.CELL_CORNER);
final ResampledGridCoverage resampled = new ResampledGridCoverage(source, complete,
MathTransforms.concatenate(targetCornerToCRS, crsToSourceCorner),
MathTransforms.concatenate(targetCenterToCRS, crsToSourceCenter),
changeOfCRS, processor);
return resampled.specialize(!isGeometryExplicit, allowOperationReplacement);
}