in endorsed/src/org.apache.sis.feature/main/org/apache/sis/coverage/grid/GridDerivation.java [672:807]
public GridDerivation subgrid(Envelope areaOfInterest, double... resolution) {
ensureSubgridNotSet();
final boolean isEnvelopeOnly = base.isEnvelopeOnly() && (resolution == null || resolution.length == 0);
MathTransform cornerToCRS = isEnvelopeOnly ? MathTransforms.identity(base.envelope.getDimension())
: base.requireGridToCRS(false); // Normal case.
subGridSetter = "subgrid";
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.
*/
MathTransform baseToAOI = null;
if (areaOfInterest != null) {
final CoordinateReferenceSystem crs = areaOfInterest.getCoordinateReferenceSystem();
if (crs != null) {
areaOfInterest = new DimensionReducer(base, crs).apply(areaOfInterest);
baseToAOI = findBaseToAOI(areaOfInterest.getCoordinateReferenceSystem());
cornerToCRS = MathTransforms.concatenate(cornerToCRS, baseToAOI);
}
}
/*
* If the grid geometry contains only an envelope, and if user asked nothing more than intersecting
* envelopes, then we will return a new GridGeometry with that intersection and nothing else.
*/
if (isEnvelopeOnly) {
if (areaOfInterest != null) {
intersection = new GeneralEnvelope(base.envelope);
if (baseToAOI != null && !baseToAOI.isIdentity()) {
areaOfInterest = Envelopes.transform(baseToAOI.inverse(), areaOfInterest);
}
intersection.intersect(areaOfInterest);
}
return this;
}
/*
* If the envelope dimensions do not encompass all grid dimensions, the transform is probably non-invertible.
* We need to reduce the number of grid dimensions in the transform for having a one-to-one relationship.
*/
int dimension = cornerToCRS.getTargetDimensions();
ArgumentChecks.ensureDimensionMatches("areaOfInterest", dimension, areaOfInterest);
cornerToCRS = dropUnusedDimensions(cornerToCRS, dimension);
/*
* Compute the sub-extent for the given Area Of Interest (AOI), ignoring for now the subsampling.
* 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).
*/
dimension = baseExtent.getDimension(); // Non-null since `base.requireGridToCRS()` succeed.
GeneralEnvelope indices = null;
if (areaOfInterest != null) {
indices = wraparound(baseToAOI, cornerToCRS).shift(areaOfInterest);
setBaseExtentClipped(false, indices);
}
if (indices == null || indices.getDimension() != dimension) {
indices = new GeneralEnvelope(dimension);
}
final GridExtent extent = getBaseExtentExpanded(true);
for (int i=0; i<dimension; i++) {
long high = extent.getHigh(i);
if (high != Long.MAX_VALUE) high++; // Increment before conversion to `double`.
indices.setRange(i, extent.getLow(i), high);
}
/*
* Convert the target resolutions to grid cell subsampling and adjust the extent consequently.
* We perform this conversion by handling the resolutions as 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 subsampling 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 affine = cornerToCRS.derivative(new DirectPositionView.Double(getPointOfInterest()));
final double[] subsampling = Matrices.inverse(affine).multiply(resolution);
@SuppressWarnings("LocalVariableHidesMemberVariable")
final int[] modifiedDimensions = this.modifiedDimensions; // Will not change anymore.
boolean scaled = false;
for (int k=0; k < subsampling.length; k++) {
double s = Math.abs(subsampling[k]);
if (s > 1) { // Also for skipping NaN values.
scaled = true;
final int i = (modifiedDimensions != null) ? modifiedDimensions[k] : k;
if (chunkSize != null && i < chunkSize.length && chunkSize[i] != 1) {
s = roundSubsampling(s, i); // Include clamp to `maximumSubsampling`.
} else {
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);
if (maximumSubsampling != null && i < maximumSubsampling.length) {
final double max = maximumSubsampling[i];
if (s > max) s = max;
}
}
indices.setRange(i, indices.getLower(i) / s,
indices.getUpper(i) / s);
}
subsampling[k] = s;
}
/*
* If at least one subsampling 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.
*
* TODO: need to clip to baseExtent, taking in account the difference in resolution.
*/
if (scaled) {
/*
* The `margin` and `chunkSize` arguments must be null because `scaledExtent` uses different units.
* The margin and chunk size were applied during `baseExtent` computation and copied in `indices`.
*/
scaledExtent = new GridExtent(indices, false, rounding, clipping, null, null, null, modifiedDimensions);
if (extent.equals(scaledExtent)) scaledExtent = extent; // Share common instance.
affine = Matrices.createIdentity(dimension + 1);
for (int k=0; k<subsampling.length; k++) {
final double s = subsampling[k];
if (s > 1) { // Also for skipping NaN values.
final int i = (modifiedDimensions != null) ? modifiedDimensions[k] : k;
affine.setElement(i, i, s);
affine.setElement(i, dimension, Math.fma(-s, scaledExtent.getLow(i), extent.getLow(i)));
}
}
toBase = MathTransforms.linear(affine);
}
}
} catch (FactoryException | TransformException e) {
throw new IllegalGridGeometryException(e, "areaOfInterest");
}
modifiedDimensions = null; // Not needed anymore.
return this;
}