Google Git
Sign in
apache / sis / 76bddb705ed107711db475c6f348b79748ab850d / . / endorsed / src / org.apache.sis.feature / main / org / apache / sis / coverage / grid / GridGeometry.java
blob: 89880ba706fc537000dc261dd913f2dd3fc3183e [file] [log] [blame]
/*
* 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 java.util.Locale;
import java.util.Objects;
import java.util.Optional;
import java.time.Instant;
import java.text.NumberFormat;
import java.math.RoundingMode;
import java.io.Serializable;
import java.io.IOException;
import java.io.UncheckedIOException;
import java.awt.image.RenderedImage; // For javadoc only.
import org.opengis.util.FactoryException;
import org.opengis.metadata.Identifier;
import org.opengis.metadata.extent.GeographicBoundingBox;
import org.opengis.geometry.Envelope;
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.operation.NoninvertibleTransformException;
import org.opengis.referencing.crs.CoordinateReferenceSystem;
import org.opengis.referencing.crs.DerivedCRS;
import org.opengis.referencing.cs.CoordinateSystemAxis;
import org.opengis.referencing.cs.CoordinateSystem;
import org.apache.sis.math.MathFunctions;
import org.apache.sis.measure.AngleFormat;
import org.apache.sis.measure.Latitude;
import org.apache.sis.measure.Longitude;
import org.apache.sis.geometry.Envelopes;
import org.apache.sis.geometry.GeneralEnvelope;
import org.apache.sis.geometry.ImmutableEnvelope;
import org.apache.sis.metadata.iso.extent.DefaultGeographicBoundingBox;
import org.apache.sis.referencing.IdentifiedObjects;
import org.apache.sis.referencing.crs.AbstractCRS;
import org.apache.sis.referencing.operation.matrix.Matrices;
import org.apache.sis.referencing.operation.matrix.MatrixSIS;
import org.apache.sis.referencing.operation.transform.MathTransforms;
import org.apache.sis.referencing.operation.transform.LinearTransform;
import org.apache.sis.referencing.operation.transform.PassThroughTransform;
import org.apache.sis.referencing.privy.ExtendedPrecisionMatrix;
import org.apache.sis.referencing.privy.DirectPositionView;
import org.apache.sis.referencing.privy.TemporalAccessor;
import org.apache.sis.referencing.privy.AxisDirections;
import org.apache.sis.metadata.privy.ReferencingServices;
import org.apache.sis.feature.internal.Resources;
import org.apache.sis.util.LenientComparable;
import org.apache.sis.util.ComparisonMode;
import org.apache.sis.util.ArgumentChecks;
import org.apache.sis.util.CharSequences;
import org.apache.sis.util.ArraysExt;
import org.apache.sis.util.Utilities;
import org.apache.sis.util.Classes;
import org.apache.sis.util.Debug;
import org.apache.sis.util.privy.DoubleDouble;
import org.apache.sis.util.privy.Numerics;
import org.apache.sis.util.collection.TreeTable;
import org.apache.sis.util.collection.TableColumn;
import org.apache.sis.util.collection.DefaultTreeTable;
import org.apache.sis.util.collection.BackingStoreException;
import org.apache.sis.util.resources.Vocabulary;
import org.apache.sis.util.resources.Messages;
import org.apache.sis.util.resources.Errors;
import org.apache.sis.util.logging.Logging;
import org.apache.sis.io.TableAppender;
import org.apache.sis.xml.NilObject;
import org.apache.sis.xml.NilReason;
import static org.apache.sis.referencing.CRS.findOperation;
// Specific to the geoapi-3.1 and geoapi-4.0 branches:
import org.opengis.coordinate.MismatchedDimensionException;
/**
* Valid extent of grid coordinates together with the transform from those grid coordinates
* to real world coordinates. {@code GridGeometry} contains:
*
* <ul class="verbose">
* <li>A {@linkplain #getExtent() grid extent} (a.k.a. <i>grid envelope</i>),
* often inferred from the {@link RenderedImage} size.</li>
* <li>A {@linkplain #getGridToCRS grid to CRS} transform,
* which may be inferred from the grid extent and the georeferenced envelope.</li>
* <li>A {@linkplain #getEnvelope() georeferenced envelope}, which can be inferred
* from the grid extent and the <i>grid to CRS</i> transform.</li>
* <li>An optional {@linkplain #getCoordinateReferenceSystem() Coordinate Reference System} (CRS)
* specified as part of the georeferenced envelope.
* This CRS is the target of the <i>grid to CRS</i> transform.</li>
* <li>An <em>estimation</em> of {@link #getResolution(boolean) grid resolution} along each CRS axes,
* computed from the <i>grid to CRS</i> transform and eventually from the grid extent.</li>
* <li>An {@linkplain #isConversionLinear indication of whether conversion for some axes is linear or not}.</li>
* </ul>
*
* The first three properties should be mandatory,
* but are allowed to be temporarily absent during grid coverage construction.
* Temporarily absent properties are allowed because they may be inferred from a wider context.
* For example, a {@code GridGeometry} knows nothing about {@link RenderedImage},
* but {@code GridCoverage2D} has this information and may use it for providing a missing grid extent.
* By default, any request for an undefined property will throw an {@link IncompleteGridGeometryException}.
* In order to check if a property is defined, use {@link #isDefined(int)}.
*
* <h2>Non-linear referencing</h2>
* A key property is the {@linkplain #getGridToCRS(PixelInCell) "grid to CRS"} conversion,
* which defines how to map pixel coordinates to "real world" coordinates such as latitudes and longitudes.
* This relationship is often linear (an affine transform), but {@linkplain #isConversionLinear not necessarily};
* {@code GridGeometry} accepts non-linear conversions as well. Non-linear conversions may occur with images
* using {@linkplain org.apache.sis.referencing.operation.builder.LocalizationGridBuilder localization grids},
* but non-linear conversions should not be used for expressing map projections (projections should be specified
* in the {@linkplain #getCoordinateReferenceSystem() Coordinate Reference System} (CRS) instead).
*
* <p>Some applications cannot handle non-linear "grid to CRS" conversions.
* For example, encoding an image in a GeoTIFF file is much simpler if the "grid to CRS" conversion is linear.
* The {@link DomainLinearizer} class can be used for replacing non-linear conversions by linear approximations.</p>
*
* <h2>Multi-threading</h2>
* {@code GridGeometry} instances are immutable and thread-safe.
* The same instance can be shared by different {@link GridCoverage} instances.
*
* @author Martin Desruisseaux (IRD, Geomatys)
* @author Johann Sorel (Geomatys)
* @version 1.5
* @since 1.0
*/
public class GridGeometry implements LenientComparable, Serializable {
/**
* Serial number for inter-operability with different versions.
*/
private static final long serialVersionUID = -954786616001606624L;
/**
* A bitmask to specify the validity of the Coordinate Reference System property.
*
* @see #isDefined(int)
* @see #getCoordinateReferenceSystem()
*/
public static final int CRS = 1;
/**
* A bitmask to specify the validity of the geodetic envelope property.
*
* @see #isDefined(int)
* @see #getEnvelope()
*/
public static final int ENVELOPE = 2;
/**
* A bitmask to specify the validity of the grid extent property.
*
* @see #isDefined(int)
* @see #getExtent()
*/
public static final int EXTENT = 4;
/**
* A bitmask to specify the validity of the <q>grid to CRS</q> transform.
*
* @see #isDefined(int)
* @see #getGridToCRS(PixelInCell)
*/
public static final int GRID_TO_CRS = 8;
/**
* A bitmask to specify the validity of the grid resolution.
*
* @see #isDefined(int)
* @see #getResolution(boolean)
*/
public static final int RESOLUTION = 16;
/**
* A bitmask to specify the validity of the geographic bounding box.
* This information can sometimes be derived from the envelope and the CRS.
* It is an optional element even with a complete grid geometry since the
* coordinate reference system is not required to have an horizontal component.
*
* @see #getGeographicExtent()
*/
public static final int GEOGRAPHIC_EXTENT = 32;
/**
* A bitmask to specify the validity of the temporal period.
* This information can sometimes be derived from the envelope and the CRS.
* It is an optional element even with a complete grid geometry since the
* coordinate reference system is not required to have a temporal component.
*
* @see #getTemporalExtent()
*/
public static final int TEMPORAL_EXTENT = 64;
/**
* The valid domain of a grid coverage, or {@code null} if unknown. The lowest valid grid coordinate is zero
* for {@link java.awt.image.BufferedImage}, but may be non-zero for arbitrary {@link RenderedImage}.
* A grid with 512 cells can have a minimum coordinate of 0 and maximum of 511.
*
* @see #EXTENT
* @see #getExtent()
*/
protected final GridExtent extent;
/**
* The geodetic envelope, or {@code null} if unknown. If non-null, this envelope is usually the grid {@link #extent}
* {@linkplain #gridToCRS transformed} to real world coordinates. The Coordinate Reference System} (CRS) of this
* envelope defines the "real world" CRS of this grid geometry.
*
* @see #ENVELOPE
* @see #getEnvelope()
*/
protected final ImmutableEnvelope envelope;
/**
* The conversion from grid indices to "real world" coordinates, or {@code null} if unknown.
* If non-null, the conversion shall map {@linkplain PixelInCell#CELL_CENTER cell center}.
* This conversion is usually, but not necessarily, affine.
*
* @see #CRS
* @see #getGridToCRS(PixelInCell)
* @see PixelInCell#CELL_CENTER
*/
@SuppressWarnings("serial") // Most SIS implementations are serializable.
protected final MathTransform gridToCRS;
/**
* Same conversion as {@link #gridToCRS} but from {@linkplain PixelInCell#CELL_CORNER cell corner}
* instead of center. This transform is preferable to {@code gridToCRS} for transforming envelopes.
*
* @serial This field is serialized because it may be a value specified explicitly at construction time,
* in which case it can be more accurate than a computed value.
*/
@SuppressWarnings("serial") // Most SIS implementations are serializable.
final MathTransform cornerToCRS;
/**
* An <em>estimation</em> of the grid resolution, in units of the CRS axes.
* Computed from {@link #gridToCRS}, eventually together with {@link #extent}.
* May be {@code null} if unknown. If non-null, the array length is equal to
* the number of CRS dimensions.
*
* @see #RESOLUTION
* @see #getResolution(boolean)
*/
protected final double[] resolution;
/**
* Whether the conversions from grid coordinates to the CRS are linear, for each target axis.
* The bit located at {@code 1L << dimension} is set to 1 when the conversion at that dimension is non-linear.
* The dimension indices are those of the CRS, not the grid. The use of {@code long} type limits the capacity
* to 64 dimensions. But actually {@code GridGeometry} can contain more dimensions provided that index of the
* last non-linear dimension is not greater than 64.
*
* @see #isConversionLinear(int...)
*/
final long nonLinears;
/**
* The geographic bounding box as an unmodifiable metadata instance, or {@code null} if not yet computed.
* If no geographic extent can be computed for the {@linkplain #envelope}, then this is set to {@link NilObject}.
*
* @see #GEOGRAPHIC_EXTENT
* @see #getGeographicExtent()
*/
private transient volatile GeographicBoundingBox geographicBBox;
/**
* The start time and end time, or {@code null} if not yet computed. If there is no time range, then the array is empty.
* If there is only a start time or an end time, then the array length is 1. Otherwise the array length is 2.
*
* @see #TEMPORAL_EXTENT
* @see #getTemporalExtent()
*/
@SuppressWarnings("VolatileArrayField") // Safe because array will not be modified after construction.
private transient volatile Instant[] timeRange;
/**
* An "empty" grid geometry with no value defined. All getter methods invoked on this instance will cause
* {@link IncompleteGridGeometryException} to be thrown. This instance can be used as a place-holder when
* the grid geometry cannot be obtained.
*/
public static final GridGeometry UNDEFINED = new GridGeometry();
/**
* Constructor for {@link #UNDEFINED} singleton only.
*/
private GridGeometry() {
extent = null;
gridToCRS = null;
cornerToCRS = null;
envelope = null;
resolution = null;
nonLinears = 0;
}
/**
* Creates a new grid geometry with the same values as the given grid geometry.
* This is a copy constructor for subclasses.
*
* @param other the other grid geometry to copy.
*/
protected GridGeometry(final GridGeometry other) {
extent = other.extent;
gridToCRS = other.gridToCRS;
cornerToCRS = other.cornerToCRS;
envelope = other.envelope;
resolution = other.resolution;
nonLinears = other.nonLinears;
}
/**
* Creates a new grid geometry derived from the given grid geometry with a new extent and a modified transform.
* This constructor is used for creating a grid geometry over a subregion (for example with the grid extent
* computed by {@link GridDerivation#subgrid(Envelope, double...)}) or grid geometry for a subsampled raster.
*
* <h4>Conversion between old and new grid geometry</h4>
* If {@code toOther} is non-null, it defines the conversion from grid coordinates of given {@code extent} to grid
* coordinates of {@code other}. That transform should be a {@linkplain MathTransforms#scale(double...) scale} and
* a {@linkplain MathTransforms#translation(double...) translation} only, but more complex transforms are accepted.
* The {@link #cornerToCRS} transform of the new grid geometry will be set to the following concatenation:
*
* <blockquote>{@code this.cornerToCRS} = {@code toOther} → {@code other.cornerToCRS}</blockquote>
*
* The new {@linkplain #getEnvelope() grid geometry envelope} will be computed from the new extent and transform,
* then {@linkplain GeneralEnvelope#intersect(Envelope) clipped} to the envelope of the other grid geometry.
* This clip is for preventing the envelope to become larger under the effect of subsampling because
* {@linkplain GridExtent#subsample(int[]) each cell become larger}. The clip is not applied when {@code toOther}
* is {@code null} because in such case, we presume that the grid extent has been changed for another reason than
* subsampling (e.g. application of a margin, in which case we want the envelope to be expanded).
*
* @param other the other grid geometry to copy.
* @param extent the new extent for the grid geometry to construct, or {@code null} if none.
* @param toOther transform from this grid coordinates to {@code other} grid coordinates, or {@code null} if none.
* @throws NullPointerException if {@code extent} is {@code null} and the other grid geometry contains no other information.
* @throws TransformException if the math transform cannot compute the geospatial envelope from the grid extent.
*
* @see GridDerivation#subgrid(Envelope, double...)
*
* @since 1.2
*/
public GridGeometry(final GridGeometry other, final GridExtent extent, final MathTransform toOther) throws TransformException {
final int dimension = other.getDimension();
this.extent = extent;
ensureDimensionMatches(dimension, extent);
if (toOther == null || toOther.isIdentity()) {
gridToCRS = other.gridToCRS;
cornerToCRS = other.cornerToCRS;
resolution = other.resolution;
nonLinears = other.nonLinears;
} else {
/*
* The `toOther` transform applies on `cornerToCRS` because the corner of upper-left pixel before scaling
* is still the corner of upper-left pixel after scaling, while "pixel center" is no longer the center of
* the same pixel. We adjust `toOther` instead of invoking `PixelTranslation.translate(cornerToCRS, …)`
* because we do not know which of `cornerToCRS` or `gridToCRS` has less NaN values.
*/
if (other.gridToCRS != null) {
final MathTransform centerShift = MathTransforms.concatenate(
MathTransforms.uniformTranslation(dimension, +0.5), toOther,
MathTransforms.uniformTranslation(dimension, -0.5));
cornerToCRS = MathTransforms.concatenate(toOther, other.cornerToCRS);
gridToCRS = MathTransforms.concatenate(centerShift, other.gridToCRS);
resolution = resolution(gridToCRS, extent, PixelInCell.CELL_CENTER);
nonLinears = findNonLinearTargets(gridToCRS);
} else {
cornerToCRS = null;
gridToCRS = null;
resolution = resolution(toOther, extent, PixelInCell.CELL_CENTER); // Save resolution even if `gridToCRS` is null.
nonLinears = findNonLinearTargets(toOther);
}
}
/*
* Recompute the envelope and clip only if a sub-sampling may have been applied (toOther != null).
* The reason for clipping is because subsampling may cause cells to appear larger.
*/
@SuppressWarnings("LocalVariableHidesMemberVariable")
ImmutableEnvelope envelope = other.envelope; // We will share the same instance if possible.
ImmutableEnvelope computed = computeEnvelope(gridToCRS, getCoordinateReferenceSystem(envelope),
toOther == null ? null : envelope); // Clip.
if (computed == null || !computed.equals(envelope)) {
envelope = computed;
}
this.envelope = envelope;
if (envelope == null && gridToCRS == null) {
ArgumentChecks.ensureNonNull("extent", extent);
}
}
/*
* Do not provide convenience constructor without PixelInCell or PixelOrientation argument.
* Experience shows that 0.5 pixel offsets in image georeferencing is a recurrent problem.
* We really want to force developers to think about whether their `gridToCRS` transform
* locates pixel corner or pixel center.
*/
/**
* Creates a new grid geometry from a grid extent and a mapping from cell coordinates to "real world" coordinates.
* At least one of {@code extent}, {@code gridToCRS} or {@code crs} arguments shall be non-null.
* If {@code gridToCRS} is non-null, then {@code anchor} shall be non-null too with one of the following values:
*
* <ul>
* <li>{@link PixelInCell#CELL_CENTER} if conversions of cell indices by {@code gridToCRS} give "real world"
* coordinates close to the center of each cell.</li>
* <li>{@link PixelInCell#CELL_CORNER} if conversions of cell indices by {@code gridToCRS} give "real world"
* coordinates at the corner of each cell. The cell corner is the one for which all grid indices have the
* smallest values (closest to negative infinity).</li>
* </ul>
*
* <h4>API note</h4>
* There is no default value for {@code anchor} because experience shows that images shifted by ½ pixel
* (with pixels that may be tens of kilometres large) is a recurrent problem. We want to encourage developers
* to always think about wether their <i>grid to CRS</i> transform is mapping pixel corner or center.
*
* <div class="warning"><b>Upcoming API generalization:</b>
* the {@code extent} type of this method may be changed to {@code GridEnvelope} interface in a future Apache SIS version.
* This is pending <a href="https://github.com/opengeospatial/geoapi/issues/36">GeoAPI update</a>.
* In addition, the {@code PixelInCell} code list currently defined in the {@code org.opengis.referencing.datum} package
* may move in another package in a future GeoAPI version because this type is no longer defined by the ISO 19111 standard
* after the 2018 revision.</div>
*
* @param extent the valid extent of grid coordinates, or {@code null} if unknown.
* @param anchor {@linkplain PixelInCell#CELL_CENTER Cell center} for OGC conventions or
* {@linkplain PixelInCell#CELL_CORNER cell corner} for Java2D/JAI conventions.
* @param gridToCRS the mapping from grid coordinates to "real world" coordinates, or {@code null} if unknown.
* @param crs the coordinate reference system of the "real world" coordinates, or {@code null} if unknown.
* @throws NullPointerException if {@code extent}, {@code gridToCRS} and {@code crs} arguments are all null.
* @throws MismatchedDimensionException if the math transform and the CRS do not have consistent dimensions.
* @throws IllegalGridGeometryException if the math transform cannot compute the geospatial envelope or resolution from the grid extent.
*/
public GridGeometry(final GridExtent extent, final PixelInCell anchor, final MathTransform gridToCRS, final CoordinateReferenceSystem crs) {
if (gridToCRS != null) {
ensureDimensionMatches(gridToCRS.getSourceDimensions(), extent);
ArgumentChecks.ensureDimensionMatches("crs", gridToCRS.getTargetDimensions(), crs);
} else if (crs == null) {
ArgumentChecks.ensureNonNull("extent", extent);
}
try {
this.extent = extent;
this.gridToCRS = PixelTranslation.translate(gridToCRS, anchor, PixelInCell.CELL_CENTER);
this.cornerToCRS = PixelTranslation.translate(gridToCRS, anchor, PixelInCell.CELL_CORNER);
this.envelope = computeEnvelope(gridToCRS, crs, null); // `gridToCRS` specified by the user, not `this.gridToCRS`.
this.resolution = resolution(gridToCRS, extent, anchor); // `gridToCRS` or `cornerToCRS` does not matter here.
this.nonLinears = findNonLinearTargets(gridToCRS);
} catch (TransformException e) {
throw new IllegalGridGeometryException(e, "gridToCRS");
}
}
/**
* Computes the envelope with the given coordinate reference system. This method is invoked from constructors.
* The {@link #extent}, {@link #gridToCRS} and {@link #cornerToCRS} fields must be set before this method is invoked.
*
* @param specified the transform specified by the user. This is not necessarily {@link #gridToCRS}.
* @param crs the coordinate reference system to declare in the envelope. May be {@code null}.
* @param limits if non-null, intersect with that envelope. The CRS must be the same as {@code crs}.
*/
private ImmutableEnvelope computeEnvelope(final MathTransform specified, final CoordinateReferenceSystem crs,
final Envelope limits) throws TransformException
{
final GeneralEnvelope env;
if (extent != null && cornerToCRS != null) {
env = extent.toEnvelope(cornerToCRS, specified, limits);
env.setCoordinateReferenceSystem(crs);
if (limits != null) {
env.intersect(limits);
}
} else if (crs != null) {
env = new GeneralEnvelope(crs);
env.setToNaN();
} else {
return null;
}
return new ImmutableEnvelope(env);
}
/**
* Creates a new grid geometry from a geospatial envelope and a mapping from cell coordinates to "real world" coordinates.
* At least one of {@code gridToCRS} or {@code envelope} arguments shall be non-null.
* If {@code gridToCRS} is non-null, then {@code anchor} shall be non-null too with one of the values documented in the
* {@link #GridGeometry(GridExtent, PixelInCell, MathTransform, CoordinateReferenceSystem) constructor expecting a grid
* extent}.
*
* <p>The given envelope shall encompass all cell surfaces, from the left border of leftmost cell to the right border
* of the rightmost cell and similarly along other axes. This constructor tries to store a geospatial envelope close
* to the specified envelope, but there is no guarantee that the envelope returned by {@link #getEnvelope()} will be
* equal to the given envelope. The envelope stored in the new {@code GridGeometry} may be slightly smaller, larger or
* shifted because the floating point values used in geospatial envelope cannot always be mapped to the integer
* coordinates used in {@link GridExtent}.
* The rules for deciding whether coordinates should be rounded toward nearest integers,
* to {@linkplain Math#floor(double) floor} or to {@linkplain Math#ceil(double) ceil} values
* are specified by the {@link GridRoundingMode} argument.</p>
*
* <p>Because of the uncertainties explained in above paragraph, this constructor should be used only in last resort,
* when the grid extent is unknown. For determinist results, developers should prefer the
* {@linkplain #GridGeometry(GridExtent, PixelInCell, MathTransform, CoordinateReferenceSystem) constructor using grid extent}
* as much as possible. In particular, this constructor is not suitable for computing grid geometry of tiles in a tiled image,
* because the above-cited uncertainties may result in apparently random black lines between tiles.</p>
*
* <div class="warning"><b>Upcoming API change:</b>
* The {@code PixelInCell} code list currently defined in the {@code org.opengis.referencing.datum} package
* may move in another package in a future GeoAPI version because this type is no longer defined by the
* ISO 19111 standard after the 2018 revision. This code list may be taken by ISO 19123 in a future revision.</div>
*
* @param anchor {@linkplain PixelInCell#CELL_CENTER Cell center} for OGC conventions or
* {@linkplain PixelInCell#CELL_CORNER cell corner} for Java2D/JAI conventions.
* @param gridToCRS the mapping from grid coordinates to "real world" coordinates, or {@code null} if unknown.
* @param envelope the geospatial envelope, including its coordinate reference system if available.
* There is no guarantee that the envelope actually stored in the {@code GridGeometry}
* will be equal to this specified envelope.
* @param rounding controls behavior of rounding from floating point values to integers.
* @throws IllegalGridGeometryException if the math transform cannot compute the grid extent or the resolution.
*/
public GridGeometry(final PixelInCell anchor, final MathTransform gridToCRS, final Envelope envelope, final GridRoundingMode rounding) {
if (gridToCRS == null) {
ArgumentChecks.ensureNonNull("envelope", envelope);
} else {
ArgumentChecks.ensureDimensionMatches("envelope", gridToCRS.getTargetDimensions(), envelope);
}
ArgumentChecks.ensureNonNull("rounding", rounding);
this.gridToCRS = PixelTranslation.translate(gridToCRS, anchor, PixelInCell.CELL_CENTER);
this.cornerToCRS = PixelTranslation.translate(gridToCRS, anchor, PixelInCell.CELL_CORNER);
Matrix scales = MathTransforms.getMatrix(gridToCRS);
int numToIgnore = 1;
if (envelope != null && cornerToCRS != null) {
GeneralEnvelope env;
try {
env = Envelopes.transform(cornerToCRS.inverse(), envelope);
extent = new GridExtent(env, rounding, GridClippingMode.STRICT, null, null, null, null);
env = extent.toEnvelope(cornerToCRS, gridToCRS, envelope);
// Use `gridToCRS` specified by the user, not `this.gridToCRS`.
} catch (TransformException e) {
throw new IllegalGridGeometryException(e, "gridToCRS");
}
env.setCoordinateReferenceSystem(envelope.getCoordinateReferenceSystem());
this.envelope = new ImmutableEnvelope(env);
if (scales == null) try {
// `gridToCRS` cannot be null if `cornerToCRS` is non-null.
scales = gridToCRS.derivative(new DirectPositionView.Double(extent.getPointOfInterest(anchor)));
numToIgnore = 0;
} catch (TransformException e) {
recoverableException("<init>", e);
}
} else {
this.extent = null;
this.envelope = ImmutableEnvelope.castOrCopy(envelope);
}
resolution = (scales != null) ? resolution(scales, numToIgnore) : null;
nonLinears = findNonLinearTargets(gridToCRS);
}
/**
* Ensures that the given dimension is equal to the expected value. If not, throws an exception.
* This method assumes that the argument name is {@code "extent"}.
*
* @param extent the extent to validate, or {@code null} if none.
* @param expected the expected number of dimension.
*/
private static void ensureDimensionMatches(final int expected, final GridExtent extent) throws MismatchedDimensionException {
if (extent != null) {
final int dimension = extent.getDimension();
if (dimension != expected) {
throw new MismatchedDimensionException(Errors.format(
Errors.Keys.MismatchedDimension_3, "extent", expected, dimension));
}
}
}
/**
* Invoked when a recoverable exception occurred. Those exceptions must be minor enough
* that they can be silently ignored in most cases.
*
* @param caller the method where exception occurred.
* @param exception the exception that occurred.
*/
static void recoverableException(final String caller, final TransformException exception) {
Logging.recoverableException(GridExtent.LOGGER, GridGeometry.class, caller, exception);
}
/**
* Creates an axis-aligned grid geometry with an extent and an envelope.
* This constructor can be used when the <i>grid to CRS</i> transform is unknown.
* If only the coordinate reference system is known, then the envelope coordinates can be
* {@linkplain GeneralEnvelope#isAllNaN() all NaN}.
*
* <p>The main purpose of this constructor is to create "desired" grid geometries, for example for use in
* {@linkplain org.apache.sis.storage.GridCoverageResource#read(GridGeometry, int...) read} or
* {@linkplain org.apache.sis.coverage.grid.GridCoverageProcessor#resample resample} operations.
* For grid geometries describing preexisting data, it is safer and more flexible to use one of
* the constructors expecting a {@link MathTransform} argument.</p>
*
* <h4>Dimension order</h4>
* The given envelope shall always declare dimensions in same order as the given extent.
* This constructor may reorder axes if {@code orientation} is (for example) {@link GridOrientation#DISPLAY},
* but in such case this constructor will derive itself an envelope and a CRS with reordered axes.
*
* @param extent the valid extent of grid coordinates, or {@code null} if unknown.
* @param envelope the envelope together with CRS of the "real world" coordinates, or {@code null} if unknown.
* @param orientation high-level description of desired characteristics of the {@code gridToCRS} transform.
* Ignored (can be null) if {@code envelope} is null.
* @throws NullPointerException if {@code extent} and {@code envelope} arguments are both null,
* or if {@code envelope} is non-null but {@code orientation} is null.
*
* @see <a href="https://en.wikipedia.org/wiki/Axis-aligned_object">Axis-aligned object on Wikipedia</a>
*
* @since 1.1
*/
public GridGeometry(GridExtent extent, final Envelope envelope, final GridOrientation orientation) {
nonLinears = 0;
/*
* Potentially change axis order and orientation according the given `GridOrientation` (which may be null).
* Current code assumes that units of measurement are unchanged, which should be true for CRSs built using
* AxisConvention.DISPLAY_ORIENTED.
*/
long flip = 0; // Bitmask specifying whether to reverse axis in each dimension.
ImmutableEnvelope target = null; // May have different axis order than the specified `envelope` CRS.
int[] sourceDimensions = null; // Indices in source envelope of axes colinear with the target envelope.
if (envelope != null) {
ArgumentChecks.ensureNonNull("orientation", orientation);
if (orientation.crsVariant != null) {
final AbstractCRS sourceCRS = AbstractCRS.castOrCopy(envelope.getCoordinateReferenceSystem());
if (sourceCRS != null) {
final AbstractCRS targetCRS = sourceCRS.forConvention(orientation.crsVariant);
if (targetCRS != sourceCRS) {
final CoordinateSystem sourceCS = sourceCRS.getCoordinateSystem();
final CoordinateSystem targetCS = targetCRS.getCoordinateSystem();
sourceDimensions = AxisDirections.indicesOfLenientMapping(sourceCS, targetCS);
if (sourceDimensions != null) {
final double[] lowerCorner = new double[sourceDimensions.length];
final double[] upperCorner = new double[sourceDimensions.length];
for (int i=0; i < sourceDimensions.length; i++) {
final int s = sourceDimensions[i];
lowerCorner[i] = envelope.getMinimum(s);
upperCorner[i] = envelope.getMaximum(s);
if (sourceCS.getAxis(s).getDirection() != targetCS.getAxis(i).getDirection()) {
flip |= Numerics.bitmask(i);
}
}
target = new ImmutableEnvelope(lowerCorner, upperCorner, targetCRS);
}
}
}
}
}
if (target == null) {
target = ImmutableEnvelope.castOrCopy(envelope);
}
/*
* If the envelope contains no useful information, then the grid extent is mandatory.
* We do that for forcing grid geometries to contain at least one information.
*/
boolean nilEnvelope = true;
if (target == null || ((nilEnvelope = target.isAllNaN()) && target.getCoordinateReferenceSystem() == null)) {
if (target != null && nilEnvelope) {
throw new NullPointerException(Errors.format(Errors.Keys.UnspecifiedCRS));
}
ArgumentChecks.ensureNonNull("extent", extent);
this.envelope = null;
} else {
this.envelope = target;
if (extent != null) {
// A non-null `sourceDimensions` implies non-null `orientation`.
if (sourceDimensions != null && orientation.canReorderGridAxis) {
if (!ArraysExt.isRange(0, sourceDimensions)) {
extent = extent.reorder(sourceDimensions);
}
sourceDimensions = null;
}
if (!nilEnvelope) {
/*
* If we have both the extent and an envelope with at least one non-NaN coordinates,
* create the `cornerToCRS` transform. The `gridToCRS` calculation uses the knowledge
* that all scale factors are on diagonal, which allows simpler calculation than full
* matrix multiplication. Use double-double arithmetic everywhere.
*/
@SuppressWarnings("null") // `!nilEnvelope` implies non-null `orientation`.
final MatrixSIS affine = extent.cornerToCRS(target, orientation.flippedAxes ^ flip, sourceDimensions);
cornerToCRS = MathTransforms.linear(affine);
final int srcDim = cornerToCRS.getSourceDimensions(); // Translation column in matrix.
final int tgtDim = cornerToCRS.getTargetDimensions(); // Number of matrix rows before last row.
resolution = new double[tgtDim];
for (int j=0; j<tgtDim; j++) {
final int i = (sourceDimensions != null) ? sourceDimensions[j] : j;
DoubleDouble scale = DoubleDouble.of(affine.getNumber(j, i), true);
DoubleDouble offset = DoubleDouble.of(affine.getNumber(j, srcDim), true);
resolution[j] = Math.abs(scale.doubleValue());
offset = offset.add(scale.scalb(-1));
affine.setNumber(j, srcDim, offset);
}
gridToCRS = MathTransforms.linear(affine);
this.extent = extent;
return;
}
}
}
this.extent = extent;
gridToCRS = null;
cornerToCRS = null;
resolution = null;
}
/**
* Creates a grid geometry with only an envelope.
*
* @param envelope the envelope together with CRS of the "real world" coordinates.
* @throws IllegalArgumentException if the envelope has no CRS and only NaN coordinate values.
*
* @since 1.4
*/
public GridGeometry(final Envelope envelope) {
// Implicit null value check below.
this.envelope = ImmutableEnvelope.castOrCopy(envelope);
if (this.envelope.isAllNaN() && this.envelope.getCoordinateReferenceSystem() == null) {
throw new IllegalArgumentException(Errors.format(Errors.Keys.EmptyArgument_1, "envelope"));
}
extent = null;
gridToCRS = null;
cornerToCRS = null;
resolution = null;
nonLinears = 0;
}
/**
* Creates a new grid geometry from the given components.
* This constructor performs no verification (unless assertions are enabled).
*/
GridGeometry(final GridExtent extent, final MathTransform gridToCRS, final MathTransform cornerToCRS,
final ImmutableEnvelope envelope, final double[] resolution, final long nonLinears)
{
this.extent = extent;
this.gridToCRS = gridToCRS;
this.cornerToCRS = cornerToCRS;
this.envelope = envelope;
this.resolution = resolution;
this.nonLinears = nonLinears;
if (gridToCRS != null) {
assert (extent == null) || gridToCRS.getSourceDimensions() == extent.getDimension();
assert (envelope == null) || gridToCRS.getTargetDimensions() == envelope.getDimension();
assert (resolution == null) || gridToCRS.getTargetDimensions() == resolution.length;
}
}
/**
* Creates a new grid geometry as the concatenation of the two specified grid geometries.
* The number of dimensions of the new grid geometry is the sum of the number of dimensions
* of the two specified grid geometries. The dimensions of the lower grid geometry are first,
* followed by the dimensions of the upper grid geometry.
*
* @param lower the grid geometry providing the lowest dimensions.
* @param upper the grid geometry providing the highest dimensions.
* @throws IncompleteGridGeometryException if a property presents in one grid geometry is absent in the other.
* @throws IllegalArgumentException if the concatenation results in duplicated
* {@linkplain GridExtent#getAxisType(int) grid axis types}.
* @throws FactoryException if the geodetic factory failed to create the compound CRS.
*
* @see #selectDimensions(int...)
*
* @since 1.5
*/
public GridGeometry(final GridGeometry lower, final GridGeometry upper) throws FactoryException {
/*
* Concatenate the extents first because they perform a check for axis duplication.
* We set the extent to null only if both `lower` and `upper` have a null extent.
* Otherwise, invoke `getExtent()` for causing an `IncompleteGridGeometryException`
* to be thrown if only one grid geometry has an extent.
*/
extent = (lower.extent == null && upper.extent == null) ? null
: new GridExtent(lower.getExtent(), upper.getExtent());
/*
* Concatenate the "grid to CRS" transforms and derived information
* (resolutions and non-linearity).
*/
if (lower.gridToCRS == null && upper.gridToCRS == null) {
gridToCRS = null;
cornerToCRS = null;
} else {
gridToCRS = compound(lower, upper, PixelInCell.CELL_CENTER);
cornerToCRS = compound(lower, upper, PixelInCell.CELL_CORNER);
}
nonLinears = lower.nonLinears | (upper.nonLinears << lower.getDimension());
resolution = (lower.resolution == null && upper.resolution == null) ? null
: ArraysExt.concatenate(lower.getResolution(true), upper.getResolution(true));
/*
* The check for presence/absence of envelope is more complex because an envelope
* may be considered missing, but still be non-null in order to store the CRS.
*/
if (lower.envelope == null && upper.envelope == null) {
envelope = null;
} else {
Envelope e1 = lower.envelope; if (e1 == null) e1 = lower.getEnvelope();
Envelope e2 = upper.envelope; if (e2 == null) e2 = upper.getEnvelope();
envelope = ImmutableEnvelope.castOrCopy(Envelopes.compound(e1, e2));
}
}
/**
* Aggregates the dimensions of the "grid to CRS" transforms of the given grid geometries.
*/
private static MathTransform compound(final GridGeometry lower, final GridGeometry upper, final PixelInCell anchor) {
return MathTransforms.compound(lower.getGridToCRS(anchor), upper.getGridToCRS(anchor));
}
/**
* Converts a "grid to CRS" transform from "cell corner" convention to "cell center" convention.
* This is a helper method for use of {@link #GridGeometry(GridExtent, MathTransform, MathTransform,
* ImmutableEnvelope, double[], long)} constructor.
*
* @param gridToCRS the transform to convert, or {@code null} if none.
* @return the converted transform, or {@code null} if the given transform was null.
*/
private static MathTransform cornerToCenter(final MathTransform gridToCRS) {
return PixelTranslation.translate(gridToCRS, PixelInCell.CELL_CORNER, PixelInCell.CELL_CENTER);
}
/**
* Returns the number of dimensions of the <em>grid</em>. This is typically the same
* than the number of {@linkplain #getEnvelope() envelope} dimensions or the number of
* {@linkplain #getCoordinateReferenceSystem() coordinate reference system} dimensions,
* but not necessarily.
*
* @return the number of grid dimensions.
*
* @see #selectDimensions(int[])
* @see GridExtent#getDimension()
*/
public final int getDimension() {
if (extent != null) {
return extent.getDimension(); // Most reliable source since that method is final.
} else if (gridToCRS != null) {
return gridToCRS.getSourceDimensions();
} else if (envelope != null) {
/*
* Last resort only since we have no guarantee that the envelope dimension is the same
* than the grid dimension (see above javadoc). The envelope should never be null at
* this point since the constructor verified that at least one argument was non-null,
* except if this method is invoked on UNDEFINED.
*/
return envelope.getDimension();
} else {
throw incomplete(EXTENT, Resources.Keys.UnspecifiedGridExtent);
}
}
/**
* Returns the number of dimensions of the <em>CRS</em>. This is typically the same as the
* number of {@linkplain #getDimension() grid dimensions}, but not necessarily.
*/
final int getTargetDimension() {
if (envelope != null) {
return envelope.getDimension(); // Most reliable source since that class is final.
} else if (gridToCRS != null) {
return gridToCRS.getTargetDimensions();
} else {
/*
* Last resort only since we have no guarantee that the grid dimension is the same
* then the CRS dimension (converse of the rational in getDimension() method).
*/
return extent.getDimension();
}
}
/**
* Returns the valid coordinate range of a grid coverage. The lowest valid grid coordinate is zero
* for {@link java.awt.image.BufferedImage}, but may be non-zero for arbitrary {@link RenderedImage}.
* A grid with 512 cells can have a minimum coordinate of 0 and maximum of 511.
*
* <div class="warning"><b>Upcoming API generalization:</b>
* the return type of this method may be changed to {@code GridEnvelope} interface in a future Apache SIS version.
* This is pending <a href="https://github.com/opengeospatial/geoapi/issues/36">GeoAPI update</a>.</div>
*
* @return the valid extent of grid coordinates (never {@code null}).
* @throws IncompleteGridGeometryException if this grid geometry has no extent —
* i.e. <code>{@linkplain #isDefined(int) isDefined}({@linkplain #EXTENT})</code> returned {@code false}.
*/
public GridExtent getExtent() {
if (extent != null) {
return extent;
}
throw incomplete(EXTENT, Resources.Keys.UnspecifiedGridExtent);
}
/**
* Returns the conversion from grid coordinates to "real world" coordinates.
* The conversion is often an affine transform, but not necessarily.
* Conversions from cell indices to geospatial coordinates can be performed for example as below:
*
* {@snippet lang="java" :
* MathTransform gridToCRS = gridGeometry.getGridToCRS(PixelInCell.CELL_CENTER);
* DirectPosition indicesOfCell = new GeneralDirectPosition(2, 3, 4):
* DirectPosition aPixelCenter = gridToCRS.transform(indicesOfCell, null);
* }
*
* Callers must specify whether they want the "real world" coordinates of cell center or cell corner.
* The cell corner is the one for which all grid indices have the smallest values (closest to negative infinity).
* As a rule of thumb:
*
* <ul>
* <li>Use {@link PixelInCell#CELL_CENTER} for transforming <em>points</em>.</li>
* <li>Use {@link PixelInCell#CELL_CORNER} for transforming <em>envelopes</em>
* with inclusive lower coordinates and <strong>exclusive</strong> upper coordinates.</li>
* </ul>
*
* <h4>API note</h4>
* There is no default value for {@code anchor} because experience shows that images shifted by ½ pixel
* (with pixels that may be tens of kilometres large) is a recurrent problem. We want to encourage developers
* to always think about wether the desired <i>grid to CRS</i> transform shall map pixel corner or center.
*
* @param anchor the cell part to map (center or corner).
* @return the conversion from grid coordinates to "real world" coordinates (never {@code null}).
* @throws IllegalArgumentException if the given {@code anchor} is not a known code list value.
* @throws IncompleteGridGeometryException if this grid geometry has no transform —
* i.e. <code>{@linkplain #isDefined(int) isDefined}({@linkplain #GRID_TO_CRS})</code> returned {@code false}.
*/
public MathTransform getGridToCRS(final PixelInCell anchor) {
final MathTransform mt;
switch (anchor) {
case CELL_CENTER: mt = gridToCRS; break;
case CELL_CORNER: mt = cornerToCRS; break;
default: mt = PixelTranslation.translate(gridToCRS, PixelInCell.CELL_CENTER, anchor); break;
}
if (mt != null) {
return mt;
}
throw incomplete(GRID_TO_CRS, Resources.Keys.UnspecifiedTransform);
}
/**
* Returns a linear approximation of the conversion from grid coordinates to "real world" coordinates.
* If the value returned by {@link #getGridToCRS(PixelInCell)} is already an instance of {@link LinearTransform},
* then it is returned as is. Otherwise this method computes the tangent of the transform at the grid extent
* {@linkplain GridExtent#getPointOfInterest(PixelInCell) point of interest} (usually the center of the grid).
*
* @param anchor the cell part to map (center or corner).
* @return linear approximation of the conversion from grid coordinates to "real world" coordinates.
* @throws IllegalArgumentException if the given {@code anchor} is not a known code list value.
* @throws IncompleteGridGeometryException if this grid geometry has no transform,
* or if the transform is non-linear but this grid geometry has no extent.
* @throws TransformException if an error occurred while computing the tangent.
*
* @since 1.3
*/
public LinearTransform getLinearGridToCRS(final PixelInCell anchor) throws TransformException {
final MathTransform tr = getGridToCRS(anchor);
if (tr instanceof LinearTransform) {
return (LinearTransform) tr;
}
return MathTransforms.linear(MathTransforms.getMatrix(tr,
new DirectPositionView.Double(getExtent().getPointOfInterest(anchor))));
}
/*
* Do not provide a convenience `getGridToCRS()` method without PixelInCell or PixelOrientation argument.
* Experience shows that 0.5 pixel offset in image localization is a recurrent problem. We really want to
* force developers to think about whether they want a gridToCRS transform locating pixel corner or center.
*/
/**
* Returns the coordinate reference system of the given envelope if defined, or {@code null} if none.
* Contrarily to {@link #getCoordinateReferenceSystem()}, this method does not throw exception.
*/
private static CoordinateReferenceSystem getCoordinateReferenceSystem(final Envelope envelope) {
return (envelope != null) ? envelope.getCoordinateReferenceSystem() : null;
}
/**
* Returns the "real world" coordinate reference system.
*
* @return the coordinate reference system (never {@code null}).
* @throws IncompleteGridGeometryException if this grid geometry has no CRS —
* i.e. <code>{@linkplain #isDefined isDefined}({@linkplain #CRS})</code> returned {@code false}.
*/
public CoordinateReferenceSystem getCoordinateReferenceSystem() {
final CoordinateReferenceSystem crs = getCoordinateReferenceSystem(envelope);
if (crs != null) return crs;
throw incomplete(CRS, Resources.Keys.UnspecifiedCRS);
}
/**
* Returns the bounding box of "real world" coordinates for this grid geometry.
* This envelope is computed from the {@linkplain #getExtent() grid extent}, which is
* {@linkplain #getGridToCRS(PixelInCell) transformed} to the "real world" coordinate system.
* The initial envelope encompasses all cell surfaces, from the left border of leftmost cell
* to the right border of the rightmost cell and similarly along other axes.
* If this grid geometry is a {@linkplain GridDerivation#subgrid(Envelope, double...) subgrid}, then the envelope is also
* {@linkplain GeneralEnvelope#intersect(Envelope) clipped} to the envelope of the original (non subsampled) grid geometry.
*
* @return the bounding box in "real world" coordinates (never {@code null}).
* @throws IncompleteGridGeometryException if this grid geometry has no envelope —
* i.e. <code>{@linkplain #isDefined(int) isDefined}({@linkplain #ENVELOPE})</code> returned {@code false}.
*/
public Envelope getEnvelope() {
if (envelope != null && !envelope.isAllNaN()) {
return envelope;
}
throw incomplete(ENVELOPE, (extent == null) ? Resources.Keys.UnspecifiedGridExtent : Resources.Keys.UnspecifiedTransform);
}
/**
* Returns the "real world" bounding box of this grid geometry transformed to the given CRS.
* This envelope is computed from the {@linkplain #getExtent() grid extent} if available,
* or from the {@linkplain #getEnvelope() envelope} otherwise.
*
* @param crs the desired coordinate reference system for the returned envelope.
* @return the bounding box in "real world" coordinates (never {@code null}).
* @throws IncompleteGridGeometryException if this grid geometry has no extent and no envelope.
* @throws TransformException if the envelope cannot be transformed to the specified CRS.
*
* @since 1.2
*/
public Envelope getEnvelope(final CoordinateReferenceSystem crs) throws TransformException {
ArgumentChecks.ensureNonNull("crs", crs);
final int bitmask; // CRS, EXTENT or GRID_TO_CRS
final short errorKey; // Resource key for error message.
final CoordinateReferenceSystem sourceCRS = getCoordinateReferenceSystem(envelope);
if (Utilities.equalsIgnoreMetadata(sourceCRS, crs)) {
return envelope;
} else if (sourceCRS == null) {
bitmask = CRS;
errorKey = Resources.Keys.UnspecifiedCRS;
} else if (extent == null && envelope == null) {
bitmask = EXTENT;
errorKey = Resources.Keys.UnspecifiedGridExtent;
} else if (cornerToCRS == null && envelope == null) {
bitmask = GRID_TO_CRS;
errorKey = Resources.Keys.UnspecifiedTransform;
} else try {
/*
* At this point the envelope should never be null because of invariants enforced by constructors.
* But we nevertheless perform some paranoiac checks. If we fail to transform the envelope, its okay.
* The main transform is the one operating on grid extent. The envelope transformation is for taking
* in account singularity points (mostly poles) and in case this grid geometry is a sub-grid geometry,
* in which case the envelope may have been clipped and we want to keep that clip.
*/
final boolean onlyEnvelope = (extent == null || cornerToCRS == null);
final CoordinateOperation op = findOperation(sourceCRS, crs, geographicBBox());
Envelope clip;
try {
clip = Envelopes.transform(op, envelope);
if (onlyEnvelope) return clip;
} catch (TransformException e) {
if (onlyEnvelope) throw e;
recoverableException("getEnvelope", e);
clip = null;
}
MathTransform tr = MathTransforms.concatenate(cornerToCRS, op.getMathTransform());
final GeneralEnvelope env = extent.toEnvelope(tr, tr, clip);
env.setCoordinateReferenceSystem(op.getTargetCRS());
env.normalize();
if (clip != null) {
env.intersect(clip);
}
return env;
} catch (FactoryException e) {
throw new TransformException(e);
}
throw incomplete(bitmask, errorKey);
}
/**
* Returns the approximate latitude and longitude coordinates of the grid.
* The prime meridian is Greenwich, but the geodetic reference frame is not necessarily WGS 84.
* This is computed from the {@linkplain #getEnvelope() envelope} if the coordinate reference system
* contains an horizontal component such as a geographic or projected CRS.
*
* <h4>API note</h4>
* This method does not throw {@link IncompleteGridGeometryException} because the geographic extent
* may be absent even with a complete grid geometry. Grid geometries are not required to have a
* spatial component on Earth surface; a raster could be a vertical profile for example.
*
* @return the geographic bounding box in "real world" coordinates.
*/
public Optional<GeographicBoundingBox> getGeographicExtent() {
return Optional.ofNullable(geographicBBox());
}
/**
* Returns the {@link #geographicBBox} value or {@code null} if none.
* This method computes the box when first needed.
*/
final GeographicBoundingBox geographicBBox() {
GeographicBoundingBox bbox = geographicBBox;
if (bbox == null) {
if (getCoordinateReferenceSystem(envelope) != null && !envelope.isAllNaN()) {
DefaultGeographicBoundingBox db;
NilReason reason = null;
try {
db = ReferencingServices.getInstance().setBounds(envelope, null, null);
} catch (TransformException e) {
db = null;
reason = NilReason.INAPPLICABLE;
}
if (db != null) {
db.transitionTo(DefaultGeographicBoundingBox.State.FINAL);
bbox = db;
} else {
if (reason == null) reason = NilReason.MISSING;
bbox = reason.createNilObject(GeographicBoundingBox.class);
}
geographicBBox = bbox;
}
}
return (bbox instanceof NilObject) ? null : bbox;
}
/**
* Returns the start time and end time of coordinates of the grid.
* If the grid has no temporal dimension, then this method returns an empty array.
* If only the start time or end time is defined, then returns an array of length 1.
* Otherwise this method returns an array of length 2 with the start time in the first element
* and the end time in the last element.
*
* @return time range as an array of length 0 (if none), 1 or 2.
*/
public Instant[] getTemporalExtent() {
Instant[] times = timeRange();
if (times.length != 0) {
times = times.clone();
}
return times;
}
/**
* Returns the {@link #timeRange} value without cloning.
* This method computes the time range when first needed.
*/
private Instant[] timeRange() {
Instant[] times = timeRange;
if (times == null) {
final TemporalAccessor t = TemporalAccessor.of(getCoordinateReferenceSystem(envelope), 0);
times = (t != null) ? t.getTimeBounds(envelope) : TemporalAccessor.EMPTY;
timeRange = times;
}
return times;
}
/**
* Returns an <em>estimation</em> of the grid resolution, in units of the coordinate reference system axes.
* The length of the returned array is the number of CRS dimensions, with {@code resolution[0]}
* being the resolution along the first CRS axis, {@code resolution[1]} the resolution along the second CRS
* axis, <i>etc</i>. Note that this axis order is not necessarily the same as grid axis order.
*
* <p>If the resolution at CRS dimension <var>i</var> is not a constant factor
* (i.e. the <code>{@linkplain #isConversionLinear(int...) isConversionLinear}(i)</code> returns {@code false}),
* then {@code resolution[i]} is set to one of the following values:</p>
*
* <ul>
* <li>{@link Double#NaN} if {@code allowEstimates} is {@code false}.</li>
* <li>An arbitrary representative resolution otherwise.
* Current implementation computes the resolution at {@link GridExtent#getPointOfInterest(PixelInCell) grid center},
* but different implementations may use alternative algorithms.</li>
* </ul>
*
* @param allowEstimates whether to provide some values even for resolutions that are not constant factors.
* @return an <em>estimation</em> of the grid resolution (never {@code null}).
* @throws IncompleteGridGeometryException if this grid geometry has no resolution —
* i.e. <code>{@linkplain #isDefined(int) isDefined}({@linkplain #RESOLUTION})</code> returned {@code false}.
*/
public double[] getResolution(final boolean allowEstimates) {
if (resolution != null) {
final double[] res = resolution.clone();
if (!allowEstimates) {
long nonLinearDimensions = nonLinears;
while (nonLinearDimensions != 0) {
final int i = Long.numberOfTrailingZeros(nonLinearDimensions);
nonLinearDimensions &= ~(1L << i);
res[i] = Double.NaN;
}
}
return res;
}
throw incomplete(RESOLUTION, (gridToCRS == null) ? Resources.Keys.UnspecifiedTransform : Resources.Keys.UnspecifiedGridExtent);
}
/**
* Computes the resolution for the given grid extent and transform, or returns {@code null} if unknown.
* Resolutions are given in order of target axes and give a scale factor from source to target coordinates.
* If the {@code gridToCRS} transform is linear, we do not even need to check the grid extent; it can be null.
* Otherwise (if the transform is non-linear) the extent is necessary. The easiest way to estimate a resolution
* is then to ask for the derivative at some arbitrary point (the point of interest).
*
* <p>Note that for this computation, it does not matter if {@code gridToCRS} is the user-specified
* transform or the {@code this.gridToCRS} field value; both should produce equivalent results.</p>
*
* @param gridToCRS a transform for which to compute the resolution, or {@code null} if none.
* @param domain the domain for which to get a resolution, or {@code null} if none.
* If non-null, must be the source of {@code gridToCRS}.
* @param anchor the pixel corner versus pixel center convention to use.
* @return the resolutions as positive numbers. May contain NaN values.
*/
static double[] resolution(final MathTransform gridToCRS, final GridExtent domain, final PixelInCell anchor) {
final Matrix matrix = MathTransforms.getMatrix(gridToCRS);
if (matrix != null) {
return resolution(matrix, 1);
} else if (domain != null && gridToCRS != null) try {
return resolution(gridToCRS.derivative(new DirectPositionView.Double(domain.getPointOfInterest(anchor))), 0);
} catch (TransformException e) {
recoverableException("resolution", e);
}
return null;
}
/**
* Computes the resolutions from the given matrix. This is the magnitude of each row vector.
* Resolutions are given in order of target axes.
*
* @param gridToCRS Jacobian matrix or affine transform for which to compute the resolution.
* @param numToIgnore number of rows and columns to ignore at the end of the matrix.
* This is 0 if the matrix is a derivative (i.e. we ignore nothing), or 1 if the matrix
* is an affine transform (i.e. we ignore the translation column and the [0 0 … 1] row).
* @return the resolutions as positive numbers. May contain NaN values.
*/
private static double[] resolution(final Matrix gridToCRS, final int numToIgnore) {
final double[] resolution = new double[gridToCRS.getNumRow() - numToIgnore];
final double[] buffer = new double[gridToCRS.getNumCol() - numToIgnore];
for (int j=0; j<resolution.length; j++) {
for (int i=0; i<buffer.length; i++) {
buffer[i] = gridToCRS.getElement(j,i);
}
resolution[j] = MathFunctions.magnitude(buffer);
}
return resolution;
}
/**
* Indicates whether the <i>grid to CRS</i> conversion is linear for all the specified CRS axes.
* The conversion from grid coordinates to real world coordinates is often linear for some dimensions,
* typically the horizontal ones at indices 0 and 1. But the vertical dimension (usually at index 2)
* is often non-linear, for example with data at 0, 5, 10, 100 and 1000 metres.
*
* @param targets indices of CRS axes. This is not necessarily the same as indices of grid axes.
* @return {@code true} if the conversion from grid coordinates to "real world" coordinates is linear
* for all the given CRS dimension.
*/
public boolean isConversionLinear(final int... targets) {
final int dimension = getTargetDimension();
long mask = 0;
for (final int d : targets) {
mask |= Numerics.bitmask(Objects.checkIndex(d, dimension));
}
return (nonLinears & mask) == 0;
}
/**
* Guesses which target dimensions may be non-linear. We currently don't have an API for finding the non-linear dimensions.
* Current implementation assumes that everything else than {@code LinearTransform} and pass-through dimensions are non-linear.
* This is not always true (e.g. in a Mercator projection, the "longitude → easting" part is linear too), but should be okay
* for {@code GridGeometry} purposes.
*
* <p>We keep trace of non-linear dimensions in a bitmask, with bits of non-linear dimensions set to 1.
* This limit us to 64 dimensions, which is assumed more than enough. Note that {@code GridGeometry} can
* contain more dimensions provided that index of the last non-linear dimension is not greater than 64.</p>
*
* @param gridToCRS the transform to "real world" coordinates, or {@code null} if unknown.
* @return a bitmask of dimensions, or 0 (i.e. conversion assumed fully linear) if the given transform was null.
*/
private static long findNonLinearTargets(final MathTransform gridToCRS) {
long nonLinearDimensions = 0;
for (final MathTransform step : MathTransforms.getSteps(gridToCRS)) {
final Matrix mat = MathTransforms.getMatrix(step);
if (mat != null) {
/*
* For linear transforms there are no bits to set. However if some bits were set by a previous
* iteration, we may need to move them (for example the transform may swap axes). We take the
* current bitmasks as source dimensions and find what are the target dimensions for them.
*/
long mask = nonLinearDimensions;
nonLinearDimensions = 0;
while (mask != 0) {
final int i = Long.numberOfTrailingZeros(mask); // Source dimension of non-linear part
for (int j = mat.getNumRow() - 1; --j >= 0;) { // Possible target dimensions
if (mat.getElement(j, i) != 0) {
if (j >= Long.SIZE) {
throw excessiveDimension(gridToCRS);
}
nonLinearDimensions |= (1L << j);
}
}
mask &= ~(1L << i);
}
} else if (step instanceof PassThroughTransform) {
/*
* Assume that all modified coordinates use non-linear transform. We do not inspect the
* sub-transform recursively because if it had a non-linear step, PassThroughTransform
* should have moved that step outside the sub-transform for easier concatenation with
* the LinearTransforms before of after that PassThroughTransform.
*/
long mask = 0;
final int dimIncrease = step.getTargetDimensions() - step.getSourceDimensions();
final int maxBits = Long.SIZE - Math.max(dimIncrease, 0);
for (final int i : ((PassThroughTransform) step).getModifiedCoordinates()) {
if (i >= maxBits) {
throw excessiveDimension(gridToCRS);
}
mask |= (1L << i);
}
/*
* The mask we just computed identifies non-linear source dimensions, but we need target
* dimensions. They are usually the same (the pass-through coordinate values do not have
* their order changed). However, we have a difficulty if the number of dimensions changes.
* We know that the change happen in the sub-transform, but we do not know where exactly.
* For example if the mask is 001010 and the number of dimensions increases by 1, we know
* that we still have "00" at the beginning and "0" at the end of the mask, but we don't
* know what happen between the two. Does "101" become "1101" or "1011"? We conservatively
* take "1111", i.e. we unconditionally set all bits in the middle to 1.
*
* Mathematics:
* (Long.highestOneBit(mask) << 1) - 1
* is a mask identifying all source dimensions before trailing pass-through dimensions.
*
* maskHigh = (Long.highestOneBit(mask) << (dimIncrease + 1)) - 1
* is a mask identifying all target dimensions before trailing pass-through dimensions.
*
* maskLow = Long.lowestOneBit(mask) - 1
* is a mask identifying all leading pass-through dimensions (both source and target).
*
* maskHigh & ~maskLow
* is a mask identifying only target dimensions after leading pass-through and before
* trailing pass-through dimensions. In our case, all 1 bits in maskLow are also 1 bits
* in maskHigh. So we can rewrite as
*
* maskHigh - maskLow
* and the -1 terms cancel each other.
*/
if (dimIncrease != 0) {
mask = (Long.highestOneBit(mask) << (dimIncrease + 1)) - Long.lowestOneBit(mask);
}
nonLinearDimensions |= mask;
} else {
/*
* Not a known transform. Assume all dimensions may become non-linear.
*/
final int dimension = gridToCRS.getTargetDimensions();
if (dimension > Long.SIZE) {
throw excessiveDimension(gridToCRS);
}
return Numerics.bitmask(dimension) - 1;
}
}
return nonLinearDimensions;
}
/**
* Invoked when the number of non-linear dimensions exceeds the {@code GridGeometry} capacity.
*/
static ArithmeticException excessiveDimension(final MathTransform gridToCRS) {
return new ArithmeticException(Errors.format(Errors.Keys.ExcessiveNumberOfDimensions_1, gridToCRS.getTargetDimensions()));
}
/**
* Invoked when a property has been requested for which we have for information.
*
* @param bitmask the requested property, for assertion purpose.
* @param errorKey the resource key to use in error message.
* @return the exception to be thrown by the caller.
*/
private IncompleteGridGeometryException incomplete(final int bitmask, final short errorKey) {
assert getClass() != GridGeometry.class || !isDefined(bitmask);
return new IncompleteGridGeometryException(Resources.format(errorKey));
}
/**
* Verifies that this grid geometry defines an {@linkplain #extent} and a {@link #cornerToCRS} transform.
* They are the information required for mapping the grid to a spatiotemporal envelope or position.
* Note that this implies that {@link #envelope} is non-null (but not necessarily that its CRS is non-null).
*
* @param center {@code true} for "center to CRS" transform, {@code false} for "corner to CRS" transform.
* @return {@link #gridToCRS} or {@link #cornerToCRS}.
*/
final MathTransform requireGridToCRS(final boolean center) throws IncompleteGridGeometryException {
if (extent == null) {
throw incomplete(EXTENT, Resources.Keys.UnspecifiedGridExtent);
}
final MathTransform mt = center ? gridToCRS : cornerToCRS;
if (mt == null) {
throw incomplete(GRID_TO_CRS, Resources.Keys.UnspecifiedTransform);
}
return mt;
}
/**
* Returns {@code true} if all the properties specified by the argument are set.
* If this method returns {@code true}, then invoking the corresponding getter
* methods will not throw {@link IncompleteGridGeometryException}.
*
* @param bitmask any combination of {@link #CRS}, {@link #ENVELOPE}, {@link #EXTENT},
* {@link #GRID_TO_CRS} and {@link #RESOLUTION}.
* @return {@code true} if all specified properties are defined (i.e. invoking the
* corresponding getter methods will not throw {@link IncompleteGridGeometryException}).
* @throws IllegalArgumentException if the specified bitmask is not a combination of known masks.
*
* @see #getCoordinateReferenceSystem()
* @see #getEnvelope()
* @see #getExtent()
* @see #getResolution(boolean)
* @see #getGridToCRS(PixelInCell)
*/
public boolean isDefined(final int bitmask) {
if ((bitmask & ~(CRS | ENVELOPE | EXTENT | GRID_TO_CRS | RESOLUTION | GEOGRAPHIC_EXTENT | TEMPORAL_EXTENT)) != 0) {
throw new IllegalArgumentException(Errors.format(Errors.Keys.IllegalArgumentValue_2, "bitmask", bitmask));
}
return ((bitmask & CRS) == 0 || (null != getCoordinateReferenceSystem(envelope)))
&& ((bitmask & ENVELOPE) == 0 || (null != envelope && !envelope.isAllNaN()))
&& ((bitmask & EXTENT) == 0 || (null != extent))
&& ((bitmask & GRID_TO_CRS) == 0 || (null != gridToCRS))
&& ((bitmask & RESOLUTION) == 0 || (null != resolution))
&& ((bitmask & GEOGRAPHIC_EXTENT) == 0 || (null != geographicBBox()))
&& ((bitmask & TEMPORAL_EXTENT) == 0 || (timeRange().length != 0));
}
/**
* Returns {@code true} if this grid geometry contains only a grid extent and no other information.
* Note: if {@link #gridToCRS} is {@code null}, then {@link #cornerToCRS} and {@link #resolution}
* should be null as well.
*/
final boolean isExtentOnly() {
return gridToCRS == null && envelope == null && extent != null;
}
/**
* Returns {@code true} if this grid geometry contains only an envelope and no other information.
* Note: if {@link #gridToCRS} is {@code null}, then {@link #cornerToCRS} and {@link #resolution}
* should be null as well.
*/
final boolean isEnvelopeOnly() {
return gridToCRS == null && extent == null && envelope != null;
}
/**
* Returns an object that can be used for creating a new grid geometry derived from this grid geometry.
* {@code GridDerivation} does not change the state of this {@code GridGeometry} but instead creates
* new instances as needed. Examples of modifications include clipping to a sub-area or applying a sub-sampling.
*
* <p>Each {@code GridDerivation} instance can be used only once and should be used in a single thread.
* {@code GridDerivation} preserves the number of dimensions. For example, {@linkplain GridDerivation#slice slicing}
* sets the {@linkplain GridExtent#getSize(int) grid size} to 1 in all dimensions specified by a <i>slice point</i>,
* but does not remove those dimensions from the grid geometry. For dimensionality reduction, see {@link #selectDimensions(int[])}.</p>
*
* <h4>Example</h4>
* For clipping this grid geometry to a sub-area, one can use:
*
* {@snippet lang="java" :
* GridGeometry gg = ...;
* Envelope areaOfInterest = ...;
* gg = gg.derive().rounding(GridRoundingMode.ENCLOSING)
* .subgrid(areaOfInterest).build();
* }
*
* @return an object for deriving a grid geometry from {@code this}.
*/
public GridDerivation derive() {
return new GridDerivation(this);
}
/**
* Creates a new grid geometry upsampled by the given number of cells along each grid dimensions.
* This method multiplies {@linkplain GridExtent#getLow(int) low} and {@linkplain GridExtent#getHigh(int) high}
* coordinates by the given periods, then scales the {@link #getGridToCRS(PixelInCell) grid to CRS} transform
* for compensating the grid change.
* The grid geometry {@link #getEnvelope() envelope} is preserved after upsampling.
*
* <h4>Number of arguments</h4>
* The {@code periods} array length should be equal to the {@linkplain #getDimension() number of grid dimensions}.
* If the array is shorter, missing values default to 1 (i.e. samplings in unspecified dimensions are unchanged).
* If the array is longer, extraneous values are ignored.
*
* @param periods the upsampling. Length shall be equal to the number of dimension and all values shall be greater than zero.
* @return the upsampled grid geometry, or {@code this} is upsampling results in the same extent.
* @throws IllegalArgumentException if a period is not greater than zero.
*
* @see GridExtent#upsample(int...)
* @since 1.3
*/
public GridGeometry upsample(final int... periods) {
GridExtent newExtent = extent;
if (newExtent != null) {
newExtent = newExtent.upsample(periods);
if (newExtent == extent) return this; // Unchanged.
}
boolean changed = false; // Additional check in case `extent` was null.
MathTransform newGridToCRS = null;
double[] newResolution = null;
if (cornerToCRS != null) {
final int tgtDim = cornerToCRS.getTargetDimensions();
final int srcDim = cornerToCRS.getSourceDimensions();
MatrixSIS matrix = Matrices.copy(MathTransforms.getMatrix(cornerToCRS));
final boolean isNonLinear = (matrix == null);
if (isNonLinear) {
matrix = Matrices.create(tgtDim+1, srcDim+1, ExtendedPrecisionMatrix.CREATE_IDENTITY);
}
/*
* By dividing the matrix elements directly, we avoid some numerical errors.
* This is because 100*0.1 is not exactly equal to 100/10.
* Each period value must be multiplied by a full column.
*/
newResolution = resolution.clone();
for (int i = Math.min(srcDim, periods.length); --i >= 0;) {
final double p = periods[i];
if (p != 1) {
newResolution[i] /= p;
final DoubleDouble pd = DoubleDouble.of(p, true);
for (int j=0; j<tgtDim; j++) {
DoubleDouble e = DoubleDouble.of(matrix.getNumber(j, i), true);
matrix.setNumber(j, i, e.divide(pd));
changed = true;
}
}
}
newGridToCRS = MathTransforms.linear(matrix);
if (isNonLinear) {
newGridToCRS = MathTransforms.concatenate(newGridToCRS, cornerToCRS);
}
}
if (!changed) return this;
return new GridGeometry(newExtent, cornerToCenter(newGridToCRS), newGridToCRS, envelope, newResolution, nonLinears);
}
/**
* Translates grid coordinates by the given number of cells without changing "real world" coordinates.
* The returned grid has the same {@linkplain GridExtent#getSize(int) size} than this grid,
* i.e. both low and high grid coordinates are displaced by the same number of cells.
* The "grid to CRS" transforms are adjusted accordingly in order to map to the same
* "real world" coordinates.
*
* <h4>Number of arguments</h4>
* The {@code translation} array length should be equal to the {@linkplain #getDimension() number of dimensions}.
* If the array is shorter, missing values default to 0 (i.e. no translation in unspecified dimensions).
* If the array is longer, extraneous values are ignored.
*
* @param translation translation to apply on each grid axis in order.
* @return a grid geometry whose coordinates (both low and high ones) and
* the "grid to CRS" transforms have been translated by given amounts.
* If the given translation is a no-op (no value or only 0 ones), then this grid is returned as is.
* @throws ArithmeticException if the translation results in coordinates that overflow 64-bits integer.
*
* @see GridExtent#translate(long...)
*
* @since 1.3
*/
public GridGeometry shiftGrid(final long... translation) {
ArgumentChecks.ensureNonNull("translation", translation);
GridExtent newExtent = extent;
if (newExtent != null) {
newExtent = newExtent.translate(translation);
if (newExtent == extent) return this;
}
MathTransform t1 = gridToCRS;
MathTransform t2 = cornerToCRS;
if (t1 != null || t2 != null) {
boolean isZero = true;
final double[] vector = new double[getDimension()];
for (int i=Math.min(vector.length, translation.length); --i >= 0;) {
isZero &= (translation[i] == 0);
vector[i] = Math.negateExact(translation[i]);
}
if (isZero) return this;
final MathTransform t = MathTransforms.translation(vector);
t1 = MathTransforms.concatenate(t, t1);
t2 = MathTransforms.concatenate(t, t2);
}
return new GridGeometry(newExtent, t1, t2, envelope, resolution, nonLinears);
}
/**
* Returns a grid geometry with the given grid extent, which implies a new "real world" computation.
* The "grid to CRS" transforms and the resolution stay the same as this {@code GridGeometry}.
* The "real world" envelope is recomputed for the new grid extent using the "grid to CRS" transforms.
*
* <p>The given extent is taken verbatim; this method does no clipping.
* The given extent does not need to intersect the extent of this grid geometry.</p>
*
* @param newExtent extent of the grid geometry to return.
* @return grid geometry with the given extent. May be {@code this} if there is no change.
* @throws TransformException if the geospatial envelope cannot be recomputed with the new grid extent.
*
* @since 1.3
*/
public GridGeometry relocate(final GridExtent newExtent) throws TransformException {
if (newExtent.equals(extent)) {
return this;
}
ensureDimensionMatches(getDimension(), newExtent);
final ImmutableEnvelope relocated;
if (cornerToCRS != null) {
final GeneralEnvelope env = newExtent.toEnvelope(cornerToCRS, gridToCRS, null);
env.setCoordinateReferenceSystem(getCoordinateReferenceSystem(envelope));
relocated = new ImmutableEnvelope(env);
} else {
relocated = envelope; // Either null or contains only the CRS.
}
return new GridGeometry(newExtent, gridToCRS, cornerToCRS, relocated, resolution, nonLinears);
}
/**
* Returns a grid geometry that encompass only some dimensions of this grid geometry.
* The specified dimensions will be copied into a new grid geometry if necessary.
* The selection is applied on {@linkplain #getExtent() grid extent} dimensions;
* they are not necessarily the same as the {@linkplain #getEnvelope() envelope} dimensions.
* The given dimensions must be in strictly ascending order without duplicated values.
* The number of dimensions of the sub grid geometry will be {@code dimensions.length}.
*
* <p>This method performs a <i>dimensionality reduction</i>.
* This method cannot be used for changing dimension order.
* The converse operation is the {@linkplain #GridGeometry(GridGeometry, GridGeometry) concatenation}.</p>
*
* @param indices the grid (not CRS) dimensions to select, in strictly increasing order.
* @return the sub-grid geometry, or {@code this} if the given array contains all dimensions of this grid geometry.
* @throws IndexOutOfBoundsException if an index is out of bounds.
*
* @see GridExtent#getSubspaceDimensions(int)
* @see GridExtent#selectDimensions(int[])
* @see org.apache.sis.referencing.CRS#selectDimensions(CoordinateReferenceSystem, int[])
*
* @since 1.3
*/
public GridGeometry selectDimensions(int... indices) {
indices = GridExtent.verifyDimensions(indices, getDimension());
if (indices != null) try {
return new SliceGeometry(this, null, indices, null).reduce(null, -1);
} catch (FactoryException e) {
throw new BackingStoreException(e);
}
return this;
}
/**
* Creates a one-, two- or three-dimensional coordinate reference system for cell indices in the grid.
* This method returns a CRS which is derived from the "real world" CRS or a subset of it.
* If the "real world" CRS is an instance of {@link org.opengis.referencing.crs.SingleCRS},
* then the derived CRS has the following properties:
*
* <ul>
* <li>{@link DerivedCRS#getBaseCRS()} is {@link #getCoordinateReferenceSystem()}.</li>
* <li>{@link DerivedCRS#getConversionFromBase()} is the inverse of {@link #getGridToCRS(PixelInCell)}.</li>
* </ul>
*
* Otherwise if the "real world" CRS is an instance of {@link org.opengis.referencing.crs.CompoundCRS},
* then only the first {@link org.opengis.referencing.crs.SingleCRS} (the head) is used.
* This is usually (but not necessarily) the horizontal component of the spatial CRS.
* The result is usually two-dimensional, but 1 and 3 dimensions are also possible.
*
* <p>Because of above relationship, it is possible to use the derived CRS in a chain of operations
* with (for example) {@link org.apache.sis.referencing.CRS#findOperation CRS.findOperation(…)}.</p>
*
* @param name name of the CRS to create.
* @param anchor the cell part to map (center or corner).
* @return a derived CRS for coordinates (cell indices) associated to the grid extent.
* @throws IncompleteGridGeometryException if the CRS, grid extent or "grid to CRS" transform is missing.
*
* @since 1.3
*/
public DerivedCRS createImageCRS(final String name, final PixelInCell anchor) {
ArgumentChecks.ensureNonEmpty("name", name);
try {
return GridExtentCRS.forCoverage(name, this, anchor, null);
} catch (FactoryException e) {
throw new BackingStoreException(e);
} catch (NoninvertibleTransformException e) {
throw new IllegalStateException(e);
}
}
/**
* Creates a transform from cell coordinates in this grid to cell coordinates in the given grid.
* The returned transform handles change of Coordinate Reference System and wraparound axes
* (e.g. longitude axis crossing the ±180° meridian) if applicable.
*
* <p><b>Note:</b> the transform created by this method may be non-invertible.</p>
*
* @param target the grid which will be the target of returned transform.
* @param anchor {@linkplain PixelInCell#CELL_CENTER Cell center} for OGC conventions or
* {@linkplain PixelInCell#CELL_CORNER cell corner} for Java2D/JAI conventions.
* @return transform from cell coordinates in this grid to cell coordinates in the given grid.
* @throws TransformException if the math transform cannot be created.
*
* @since 1.1
*/
public MathTransform createTransformTo(final GridGeometry target, final PixelInCell anchor) throws TransformException {
ArgumentChecks.ensureNonNull("target", target);
ArgumentChecks.ensureNonNull("anchor", anchor);
/*
* Inverse `source` and `target` because `CoordinateOperationFinder.inverse(…)` does an
* effort for using `WraparoundTransform` only if needed (contrarily to `gridToCRS(…)`).
*/
final CoordinateOperationFinder finder = new CoordinateOperationFinder(target, this);
finder.verifyPresenceOfCRS(false);
finder.setAnchor(anchor);
final MathTransform tr;
try {
tr = finder.inverse();
} catch (FactoryException e) {
throw new TransformException(e);
}
return MathTransforms.concatenate(getGridToCRS(anchor), tr);
}
/**
* Returns a hash value for this grid geometry. This value needs not to remain
* consistent between different implementations of the same class.
*/
@Override
public int hashCode() {
int code = (int) serialVersionUID;
if (gridToCRS != null) {
code += gridToCRS.hashCode();
}
if (extent != null) {
code += extent.hashCode();
}
// We do not check the envelope since it has a determinist relationship with other attributes.
return code;
}
/**
* Compares the specified object with this grid geometry for equality.
* This method delegates to {@code equals(object, ComparisonMode.STRICT)}.
*
* @param object the object to compare with.
* @return {@code true} if the given object is equal to this grid geometry.
*/
@Override
public boolean equals(final Object object) {
return equals(object, ComparisonMode.STRICT);
}
/**
* Compares the specified object with this grid geometry for equality.
* If the mode is {@link ComparisonMode#IGNORE_METADATA} or more flexible,
* then the {@linkplain GridExtent#getAxisType(int) axis types} are ignored.
*
* @param object the object to compare with this grid geometry for equality.
* @param mode the strictness level of the comparison.
* @return {@code true} if the given object is equal to this grid geometry.
*
* @since 1.1
*/
@Override
public boolean equals(final Object object, final ComparisonMode mode) {
if (object == this) {
return true;
}
if (object instanceof GridGeometry) {
final GridGeometry that = (GridGeometry) object;
if ((mode != ComparisonMode.STRICT || getClass().equals(object.getClass()))
&& Utilities.deepEquals(extent, that.extent, mode)
&& Utilities.deepEquals(gridToCRS, that.gridToCRS, mode))
{
final ImmutableEnvelope othenv = that.envelope;
if (!mode.isApproximate()) {
return Utilities.deepEquals(envelope, othenv, mode);
}
if ((envelope == null) == (othenv == null) &&
Utilities.deepEquals(getCoordinateReferenceSystem(envelope),
getCoordinateReferenceSystem(othenv), mode))
{
return equalsApproximately(othenv);
}
}
}
return false;
}
/**
* Returns whether the given envelope is equal to this grid geometry envelope with a tolerance threshold
* computed from grid resolution. If this grid geometry has no envelope, then this method arbitrarily
* returns {@code true} (this unusual behavior is required by {@link #equals(Object, ComparisonMode)}).
*/
final boolean equalsApproximately(final ImmutableEnvelope othenv) {
if (envelope != null) {
for (int i=envelope.getDimension(); --i >= 0;) {
// Arbitrary threshold of ½ pixel.
final double ε = (resolution != null) ? resolution[i] * 0.5 : 0;
if (!MathFunctions.epsilonEqual(envelope.getLower(i), othenv.getLower(i), ε) ||
!MathFunctions.epsilonEqual(envelope.getUpper(i), othenv.getUpper(i), ε))
{
return false;
}
}
}
return true;
}
/**
* Returns the default set of flags to use for {@link #toString()} implementations.
* Current implementation returns all core properties, augmented with only derived
* properties that are defined.
*/
final int defaultFlags() {
int flags = EXTENT | GRID_TO_CRS | CRS;
if (null != envelope) flags |= ENVELOPE;
if (null != resolution) flags |= RESOLUTION;
if (null != geographicBBox()) flags |= GEOGRAPHIC_EXTENT;
if (timeRange().length != 0) flags |= TEMPORAL_EXTENT;
return flags;
}
/**
* Returns a string representation of this grid geometry.
* The returned string is implementation dependent and may change in any future version.
* Current implementation is equivalent to a call to {@link #toTree(Locale, int)} with
* at least {@link #EXTENT}, {@link #ENVELOPE} and {@link #CRS} flags.
* Whether more flags are present or not is unspecified.
*/
@Override
public String toString() {
return toTree(Locale.getDefault(), defaultFlags()).toString();
}
/**
* Returns a tree representation of some elements of this grid geometry.
* The tree representation is for debugging or logging purposes
* and may change in any future SIS version.
*
* @param locale the locale to use for textual labels.
* @param bitmask combination of {@link #EXTENT}, {@link #ENVELOPE}, {@link #CRS}, {@link #GRID_TO_CRS},
* {@link #RESOLUTION}, {@link #GEOGRAPHIC_EXTENT} and {@link #TEMPORAL_EXTENT}.
* @return a tree representation of the specified elements.
*/
@Debug
public TreeTable toTree(final Locale locale, final int bitmask) {
ArgumentChecks.ensureNonNull("locale", locale);
final TreeTable tree = new DefaultTreeTable(TableColumn.VALUE_AS_TEXT);
final TreeTable.Node root = tree.getRoot();
root.setValue(TableColumn.VALUE_AS_TEXT, Classes.getShortClassName(this));
formatTo(locale, Vocabulary.forLocale(locale), bitmask, root);
return tree;
}
/**
* Formats a string representation of this grid geometry in the specified tree.
*/
final void formatTo(final Locale locale, final Vocabulary vocabulary, final int bitmask, final TreeTable.Node root) {
if ((bitmask & ~(EXTENT | ENVELOPE | CRS | GRID_TO_CRS | RESOLUTION | GEOGRAPHIC_EXTENT | TEMPORAL_EXTENT)) != 0) {
throw new IllegalArgumentException(Errors.format(
Errors.Keys.IllegalArgumentValue_2, "bitmask", bitmask));
}
try {
new Formatter(locale, vocabulary, bitmask, root).format();
} catch (IOException e) {
throw new UncheckedIOException(e);
}
}
/**
* Helper class for formatting a {@link GridGeometry} instance.
*/
private final class Formatter {
/**
* Combination of {@link #EXTENT}, {@link #ENVELOPE}, {@link #CRS}, {@link #GRID_TO_CRS},
* {@link #RESOLUTION}, {@link #GEOGRAPHIC_EXTENT} and {@link #TEMPORAL_EXTENT}.
*/
private final int bitmask;
/**
* Temporary buffer for formatting node values.
*/
private final StringBuilder buffer;
/**
* Where to write the {@link GridGeometry} string representation.
*/
private final TreeTable.Node root;
/**
* The section under the {@linkplain #root} where to write elements.
* This is updated when {@link #section(int, short, boolean, boolean)} is invoked.
*/
private TreeTable.Node section;
/**
* Localized words.
*/
private final Vocabulary vocabulary;
/**
* The locale for the texts, numbers (except grid extent) and dates.
*/
private final Locale locale;
/**
* The coordinate reference system, or {@code null} if none.
*/
private final CoordinateReferenceSystem crs;
/**
* The coordinate system, or {@code null} if none.
*/
private final CoordinateSystem cs;
/**
* Creates a new formatter for the given combination of {@link #EXTENT}, {@link #ENVELOPE},
* {@link #CRS}, {@link #GRID_TO_CRS} and {@link #RESOLUTION}.
*/
Formatter(final Locale locale, final Vocabulary vocabulary, final int bitmask, final TreeTable.Node out) {
this.root = out;
this.bitmask = bitmask;
this.buffer = new StringBuilder(256);
this.locale = locale;
this.vocabulary = vocabulary;
this.crs = getCoordinateReferenceSystem(envelope);
this.cs = (crs != null) ? crs.getCoordinateSystem() : null;
}
/**
* Formats a string representation of the enclosing {@link GridGeometry} instance
* in the buffer specified at construction time.
*/
final void format() throws IOException {
/*
* Example: Grid extent
* ├─ Dimension 0: [370 … 389] (20 cells)
* └─ Dimension 1: [ 41 … 340] (300 cells)
*
* We need to use the implementation provided by GridExtent in order
* to format correctly the unsigned long numbers for very large sizes.
*/
if (section(EXTENT, Vocabulary.Keys.GridExtent, true, false)) {
extent.appendTo(buffer, vocabulary);
writeNodes();
}
/*
* Example: Geographic extent
*  ├─Lower bound: 80°00′00″S 180°00′00″W 2019-05-01T21:00:00Z
*  └─Upper bound: 90°00′00″N 180°00′00″E 2019-05-02T00:00:00Z
*
* The angle and date/time patterns are fixed for now, with a precision equivalent to about 30 metres.
* The angles are rounded toward up and down for making sure that the box encloses fully the coverage.
*/
if (section(GEOGRAPHIC_EXTENT, Vocabulary.Keys.GeographicExtent, false, false) ||
section(TEMPORAL_EXTENT, Vocabulary.Keys.TemporalExtent, false, false))
{
final TableAppender table = new TableAppender(buffer, " ");
final AngleFormat nf = new AngleFormat("DD°MM′SS″", locale);
final GeographicBoundingBox bbox = ((bitmask & GEOGRAPHIC_EXTENT) != 0) ? geographicBBox() : null;
double westBoundLongitude = Double.NaN;
double eastBoundLongitude = Double.NaN;
final Instant[] times = ((bitmask & TEMPORAL_EXTENT) != 0) ? timeRange() : TemporalAccessor.EMPTY;
vocabulary.appendLabel(Vocabulary.Keys.LowerBound, table);
table.setCellAlignment(TableAppender.ALIGN_RIGHT);
if (bbox != null) {
nf.setRoundingMode(RoundingMode.FLOOR);
westBoundLongitude = bbox.getWestBoundLongitude();
table.nextColumn(); table.append(nf.format(new Latitude(bbox.getSouthBoundLatitude())));
table.nextColumn(); table.append(nf.format(new Longitude(westBoundLongitude)));
}
if (times.length >= 1) {
table.nextColumn();
table.append(times[0].toString());
}
table.nextLine();
table.setCellAlignment(TableAppender.ALIGN_LEFT);
vocabulary.appendLabel(Vocabulary.Keys.UpperBound, table);
table.setCellAlignment(TableAppender.ALIGN_RIGHT);
if (bbox != null) {
nf.setRoundingMode(RoundingMode.CEILING);
eastBoundLongitude = bbox.getEastBoundLongitude();
table.nextColumn(); table.append(nf.format(new Latitude(bbox.getNorthBoundLatitude())));
table.nextColumn(); table.append(nf.format(new Longitude(eastBoundLongitude)));
}
if (times.length >= 2) {
table.nextColumn();
table.append(times[1].toString());
}
table.flush();
if (Longitude.isWraparound(westBoundLongitude, eastBoundLongitude)) {
vocabulary.appendLabel(Vocabulary.Keys.Note, buffer);
buffer.append(' ')
.append(Messages.forLocale(locale).getString(Messages.Keys.BoxCrossesAntiMeridian));
}
writeNodes();
}
/*
* Example: Envelope
* ├─ Geodetic latitude: -69.75 … 80.25 ∆φ = 0.5°
* └─ Geodetic longitude: 4.75 … 14.75 ∆λ = 0.5°
*
* The minimum number of fraction digits is the number required for differentiating two consecutive cells.
* The maximum number of fraction digits avoids to print more digits than the precision of `double` type.
* Those numbers vary for each line depending on the envelope values and the resolution at that line.
*/
if (section(ENVELOPE, Vocabulary.Keys.Envelope, true, false)) {
final boolean appendResolution = (bitmask & RESOLUTION) != 0 && resolution != null;
final TableAppender table = new TableAppender(buffer, "");
final int dimension = envelope.getDimension();
final NumberFormat nf = NumberFormat.getNumberInstance(locale);
for (int i=0; i<dimension; i++) {
final double lower = envelope.getLower(i);
final double upper = envelope.getUpper(i);
final double delta = (resolution != null) ? resolution[i] : Double.NaN;
nf.setMinimumFractionDigits(Numerics.fractionDigitsForDelta(delta));
nf.setMaximumFractionDigits(Numerics.suggestFractionDigits(lower, upper) - 1); // The -1 is for rounding errors.
final CoordinateSystemAxis axis = (cs != null) ? cs.getAxis(i) : null;
final String name = (axis != null) ? axis.getName().getCode() : vocabulary.getString(Vocabulary.Keys.Dimension_1, i);
table.append(name).append(": ").nextColumn();
table.setCellAlignment(TableAppender.ALIGN_RIGHT);
table.append(nf.format(lower)).nextColumn();
table.setCellAlignment(TableAppender.ALIGN_LEFT);
table.append(" … ").append(nf.format(upper));
if (appendResolution) {
final boolean isLinear = (i < Long.SIZE) && (nonLinears & (1L << i)) == 0;
table.nextColumn();
table.append(" ∆");
if (axis != null) {
table.append(axis.getAbbreviation());
}
table.nextColumn();
table.append(' ').append(isLinear ? '=' : '≈').append(' ');
appendResolution(table, nf, delta, i);
}
table.nextLine();
}
table.flush();
writeNodes();
} else if (section(RESOLUTION, Vocabulary.Keys.Resolution, true, false)) {
/*
* Example: Resolution
* └─ 0.5° × 0.5°
*
* Formatted only as a fallback if the envelope was not formatted.
* Otherwise, this information is already part of above envelope.
*/
String separator = "";
final NumberFormat nf = NumberFormat.getNumberInstance(locale);
for (int i=0; i < resolution.length; i++) {
appendResolution(buffer.append(separator), nf, resolution[i], i);
separator = " × ";
}
writeNode();
}
/*
* Example: Coordinate reference system
* └─ EPSG:4326 — WGS 84 (φ,λ)
*/
if (section(CRS, Vocabulary.Keys.CoordinateRefSys, true, false)) {
final Identifier id = IdentifiedObjects.getIdentifier(crs, null);
if (id != null) {
buffer.append(IdentifiedObjects.toString(id)).append(" — ");
}
buffer.append(crs.getName().getCode());
writeNode();
}
/*
* Example: Conversion
* └─ 2D → 2D non linear in 2
*/
final Matrix matrix = MathTransforms.getMatrix(gridToCRS);
if (section(GRID_TO_CRS, Vocabulary.Keys.Conversion, true, matrix != null)) {
if (matrix != null) {
writeNode(Matrices.toString(matrix));
} else {
buffer.append(gridToCRS.getSourceDimensions()).append("D → ")
.append(gridToCRS.getTargetDimensions()).append("D ");
long nonLinearDimensions = nonLinears;
String separator = Resources.forLocale(locale)
.getString(Resources.Keys.NonLinearInDimensions_1, Long.bitCount(nonLinearDimensions));
while (nonLinearDimensions != 0) {
final int i = Long.numberOfTrailingZeros(nonLinearDimensions);
nonLinearDimensions &= ~(1L << i);
buffer.append(separator).append(' ')
.append(cs != null ? cs.getAxis(i).getName().getCode() : String.valueOf(i));
separator = ",";
}
writeNode();
}
}
}
/**
* Starts a new section for the given property.
*
* @param property one of {@link #EXTENT}, {@link #ENVELOPE}, {@link #CRS}, {@link #GRID_TO_CRS} and {@link #RESOLUTION}.
* @param title the {@link Vocabulary} key for the title to show for this section, if formatted.
* @param mandatory whether to write "undefined" if the property is undefined.
* @param cellCenter whether to add a "origin in cell center" text in the title. This is relevant only for conversion.
* @return {@code true} if the caller shall format the value.
*/
private boolean section(final int property, final short title, final boolean mandatory, final boolean cellCenter) {
if ((bitmask & property) != 0) {
CharSequence text = vocabulary.getString(title);
if (cellCenter) {
text = buffer.append(text).append(" (")
.append(vocabulary.getString(Vocabulary.Keys.OriginInCellCenter).toLowerCase(locale))
.append(')').toString();
buffer.setLength(0);
}
section = root.newChild();
section.setValue(TableColumn.VALUE_AS_TEXT, text);
if (isDefined(property)) {
return true;
}
if (mandatory) {
writeNode(vocabulary.getString(Vocabulary.Keys.Unspecified));
}
}
return false;
}
/**
* Appends a single line as a node in the current section.
*/
private void writeNode(final CharSequence line) {
String text = line.toString().trim();
if (!text.isEmpty()) {
section.newChild().setValue(TableColumn.VALUE_AS_TEXT, text);
}
}
/**
* Appends a node with current {@link #buffer} content as a single line, then clears the buffer.
*/
private void writeNode() {
writeNode(buffer);
buffer.setLength(0);
}
/**
* Appends nodes with current {@link #buffer} content as multi-lines text, then clears the buffer.
*/
private void writeNodes() {
for (final CharSequence line : CharSequences.splitOnEOL(buffer)) {
writeNode(line);
}
buffer.setLength(0);
}
/**
* Appends a single value on the resolution line, together with its unit of measurement.
*
* @param out where to write the resolution.
* @param nf number format to use for writing the number.
* @param res the resolution to write, or {@link Double#NaN}.
* @param dim index of the coordinate system axis of the resolution.
*/
private void appendResolution(final Appendable out, final NumberFormat nf, final double res, final int dim) throws IOException {
if (Double.isNaN(res)) {
out.append('?');
} else {
nf.setMaximumFractionDigits(Numerics.suggestFractionDigits(res) / 2);
out.append(nf.format(res));
}
if (cs != null) {
final String unit = String.valueOf(cs.getAxis(dim).getUnit());
if (unit.isEmpty() || Character.isLetterOrDigit(unit.codePointAt(0))) {
out.append(' ');
}
out.append(unit);
}
}
}
}