/*
 * 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.Arrays;
import java.util.Objects;
import java.util.Locale;
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.geometry.MismatchedDimensionException;
import org.opengis.referencing.datum.PixelInCell;
import org.opengis.referencing.operation.Matrix;
import org.opengis.referencing.operation.MathTransform;
import org.opengis.referencing.operation.TransformException;
import org.opengis.referencing.crs.CoordinateReferenceSystem;
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.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.TransformSeparator;
import org.apache.sis.referencing.operation.transform.PassThroughTransform;
import org.apache.sis.internal.referencing.DirectPositionView;
import org.apache.sis.internal.referencing.TemporalAccessor;
import org.apache.sis.internal.metadata.ReferencingServices;
import org.apache.sis.internal.system.Modules;
import org.apache.sis.internal.feature.Resources;
import org.apache.sis.internal.util.DoubleDouble;
import org.apache.sis.internal.util.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.resources.Vocabulary;
import org.apache.sis.util.resources.Errors;
import org.apache.sis.util.logging.Logging;
import org.apache.sis.util.ArgumentChecks;
import org.apache.sis.util.CharSequences;
import org.apache.sis.util.ArraysExt;
import org.apache.sis.util.Classes;
import org.apache.sis.util.Debug;
import org.apache.sis.io.TableAppender;
import org.apache.sis.xml.NilObject;
import org.apache.sis.xml.NilReason;


/**
 * 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. <cite>grid envelope</cite>),
 *       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 <cite>grid to CRS</cite> 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 <cite>grid to CRS</cite> transform.</li>
 *   <li>An <em>estimation</em> of {@link #getResolution(boolean) grid resolution} along each CRS axes,
 *       computed from the <cite>grid to CRS</cite> 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 grid geometry know nothing about {@link RenderedImage},
 * but {@code GridCoverage2D} does and may use that information 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)}.
 *
 * <p>{@code GridGeometry} instances are immutable and thread-safe.
 * The same instance can be shared by different {@link GridCoverage} instances.</p>
 *
 * @author  Martin Desruisseaux (IRD, Geomatys)
 * @version 1.0
 * @since   1.0
 * @module
 */
public class GridGeometry implements 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 <cite>"grid to CRS"</cite> 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 sometime 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 sometime 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
     */
    protected final MathTransform gridToCRS;

    /**
     * Same conversion than {@link #gridToCRS} but from {@linkplain PixelInCell#CELL_CORNER cell corner}
     * instead than 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.
     */
    private 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.
     *
     * @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...)
     */
    private 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 can not 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 than 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.
     *
     * <p>If {@code toOther} is non-null, it should be a transform from the given {@code extent} coordinates to the
     * {@code other} grid coordinates. That transform should be merely a {@linkplain MathTransforms#scale(double...)
     * scale} and {@linkplain MathTransforms#translation(double...) translation} even if more complex transforms are
     * accepted. The {@link #cornerToCRS} transform of the new grid geometry will be set to the following concatenation:</p>
     *
     * <blockquote>{@code this.cornerToCRS} = {@code toOther} → {@code other.cornerToCRS}</blockquote>
     *
     * The new {@linkplain #getEnvelope() grid geometry envelope} will be {@linkplain GeneralEnvelope#intersect(Envelope)
     * clipped} to the envelope of the other grid geometry. This is for preventing the envelope to become larger under the
     * effect of subsampling (because {@link GridExtent#subsample(int[]) each cell become larger}).
     *
     * @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 can not compute the geospatial envelope from the grid extent.
     *
     * @see GridDerivation#subgrid(Envelope, double...)
     */
    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 than invoking `PixelTranslation.translate(cornerToCRS, …)`
             * because we do not know which of `cornerToCRS` or `gridToCRS` has less NaN values.
             */
            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);
            nonLinears  = findNonLinearTargets(gridToCRS);
        }
        ImmutableEnvelope envelope = other.envelope;            // We will share the same instance if possible.
        ImmutableEnvelope computed = computeEnvelope(gridToCRS, getCoordinateReferenceSystem(envelope), envelope);
        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>
     *
     * <div class="note"><b>API note:</b>
     * 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 <cite>grid to CRS</cite> transform is mapping pixel corner or center.</div>
     *
     * <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 can not 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);           // '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.
     * @param  limits     if non-null, intersect with that envelope. The CRS must be the same than {@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.toCRS(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 guarantees 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 can not 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 guarantees 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 can not compute the grid extent or the resolution.
     */
    @SuppressWarnings("null")
    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, null, null, null);
                env = extent.toCRS(cornerToCRS, gridToCRS, envelope);     // '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' can not be null if 'cornerToCRS' is non-null.
                scales = gridToCRS.derivative(new DirectPositionView.Double(extent.getPointOfInterest()));
                numToIgnore = 0;
            } catch (TransformException e) {
                recoverableException(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 equals 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  exception  the exception that occurred.
     */
    static void recoverableException(final Exception exception) {
        Logging.recoverableException(Logging.getLogger(Modules.RASTER), GridGeometry.class, "<init>", exception);
    }

    /**
     * Creates a grid geometry with an extent and an envelope.
     * This constructor can be used when the <cite>grid to CRS</cite> transform is unknown.
     * If only the coordinate reference system is known, then the envelope coordinates can be
     * {@link GeneralEnvelope#isAllNaN() NaN}.
     *
     * <p>This constructor is generally not recommended since creating a <cite>grid to CRS</cite> from an envelope
     * requires assumption on axis order and axis directions. This constructor assumes that all grid axes are in the
     * same order than CRS axes and no axis is flipped. This straightforward approach often results in the <var>y</var>
     * axis to be oriented toward up, not down as expected in rendered images. Those assumptions are often not suitable.
     * For better control, use one of the constructors expecting a {@link MathTransform} argument instead.
     * This constructor is provided mostly as a convenience for testing purposes, or when only the extent is known.</p>
     *
     * @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.
     * @throws NullPointerException if {@code extent} and {@code envelope} arguments are both null.
     */
    public GridGeometry(final GridExtent extent, final Envelope envelope) {
        this.extent = extent;
        nonLinears  = 0;
        /*
         * 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;
        final ImmutableEnvelope env = ImmutableEnvelope.castOrCopy(envelope);
        if (env == null || ((nilEnvelope = env.isAllNaN()) && env.getCoordinateReferenceSystem() == null)) {
            ArgumentChecks.ensureNonNull("extent", extent);
            this.envelope = null;
        } else {
            this.envelope = env;
            if (extent != null && !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 with no sign reversal, which allows simpler
                 * calculation than full matrix multiplication. Use double-double arithmetic everywhere.
                 */
                final MatrixSIS affine = extent.cornerToCRS(env);
                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 DoubleDouble scale  = (DoubleDouble) affine.getNumber(j, j);
                    final DoubleDouble offset = (DoubleDouble) affine.getNumber(j, srcDim);
                    resolution[j] = scale.doubleValue();
                    scale.multiply(0.5);
                    offset.add(scale);
                    affine.setNumber(j, srcDim, offset);
                }
                gridToCRS = MathTransforms.linear(affine);
                return;
            }
        }
        gridToCRS   = null;
        cornerToCRS = null;
        resolution  = null;
    }

    /**
     * Creates a new grid geometry over the specified dimensions of the given grid geometry.
     *
     * @param  other       the grid geometry to copy.
     * @param  dimensions  the dimensions to select, in strictly increasing order (this may not be verified).
     * @throws FactoryException if an error occurred while separating the "grid to CRS" transform.
     *
     * @see #reduce(int...)
     */
    private GridGeometry(final GridGeometry other, int[] dimensions) throws FactoryException {
        extent = (other.extent != null) ? other.extent.reduce(dimensions) : null;
        final int n = dimensions.length;
        if (other.gridToCRS != null) {
            final int[] sources = dimensions;
            TransformSeparator sep = new TransformSeparator(other.gridToCRS);
            sep.addSourceDimensions(sources);
            gridToCRS  = sep.separate();
            dimensions = sep.getTargetDimensions();
            assert dimensions.length == n : Arrays.toString(dimensions);
            if (!ArraysExt.isSorted(dimensions, true)) {
                throw new IllegalGridGeometryException(Resources.format(Resources.Keys.IllegalGridGeometryComponent_1, "dimensions"));
            }
            /*
             * We redo a separation for 'cornerToCRS' instead than applying a translation of the 'gridToCRS'
             * computed above because we don't know which of 'gridToCRS' and 'cornerToCRS' has less NaN values.
             */
            sep = new TransformSeparator(other.cornerToCRS);
            sep.addSourceDimensions(sources);
            sep.addTargetDimensions(dimensions);
            cornerToCRS = sep.separate();
            assert Arrays.equals(sep.getSourceDimensions(), dimensions) : Arrays.toString(dimensions);
        } else {
            gridToCRS   = null;
            cornerToCRS = null;
        }
        final ImmutableEnvelope env = other.envelope;
        if (env != null) {
            CoordinateReferenceSystem crs = env.getCoordinateReferenceSystem();
            crs = org.apache.sis.referencing.CRS.reduce(crs, dimensions);
            final double[] min = new double[n];
            final double[] max = new double[n];
            for (int i=0; i<n; i++) {
                final int j = dimensions[i];
                min[i] = env.getLower(j);
                max[i] = env.getUpper(j);
            }
            envelope = new ImmutableEnvelope(min, max, crs);
        } else {
            envelope = null;
        }
        long     nonLinears = 0;
        double[] resolution = other.resolution;
        if (resolution != null) {
            resolution = new double[n];
        }
        for (int i=0; i<n; i++) {
            final int j = dimensions[i];
            if (resolution != null) {
                resolution[i] = other.resolution[j];
            }
            nonLinears |= ((other.nonLinears >>> j) & 1L) << i;
        }
        this.resolution = resolution;
        this.nonLinears = nonLinears;
    }

    /**
     * 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 #reduce(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 {
            /*
             * 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.
             */
            return envelope.getDimension();
        }
    }

    /**
     * Returns the number of dimensions of the <em>CRS</em>. This is typically the same than the
     * number of {@linkplain #getDimension() grid dimensions}, but not necessarily.
     */
    private 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:
     *
     * {@preformat 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>
     *
     * <div class="note"><b>API note:</b>
     * 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 <cite>grid to CRS</cite> transform shall map pixel corner or center.</div>
     *
     * @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;
        if (PixelInCell.CELL_CENTER.equals(anchor)) {
            mt = gridToCRS;
        } else if (PixelInCell.CELL_CORNER.equals(anchor)) {
            mt = cornerToCRS;
        }  else {
            mt = PixelTranslation.translate(gridToCRS, PixelInCell.CELL_CENTER, anchor);
        }
        if (mt != null) {
            return mt;
        }
        throw incomplete(GRID_TO_CRS, Resources.Keys.UnspecifiedTransform);
    }

    /*
     * 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 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.
     *
     * <div class="note"><b>API note:</b>
     * this method does not throw {@link IncompleteGridGeometryException} because the geographic extent may be absent
     * even with a complete grid geometry. This is because grid geometries are not required to have an spatial component
     * on Earth surface; a raster could be a vertical profile for example.</div>
     *
     * @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.
     */
    private GeographicBoundingBox geographicBBox() {
        GeographicBoundingBox bbox = geographicBBox;
        if (bbox == null) {
            if (getCoordinateReferenceSystem(envelope) != null && !envelope.isAllNaN()) {
                try {
                    final DefaultGeographicBoundingBox db = ReferencingServices.getInstance().setBounds(envelope, null, null);
                    db.transitionTo(DefaultGeographicBoundingBox.State.EDITABLE);
                    bbox = db;
                } catch (TransformException e) {
                    bbox = NilReason.INAPPLICABLE.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.getTimeRange(envelope) : new Instant[0];
            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 than 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() 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}.
     * @return the resolutions as positive numbers. May contain NaN values.
     */
    static double[] resolution(final MathTransform gridToCRS, final GridExtent domain) {
        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())), 0);
        } catch (TransformException e) {
            recoverableException(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 <cite>grid to CRS</cite> 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 than 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) {
            ArgumentChecks.ensureValidIndex(dimension, d);
            if (d < Long.SIZE) mask |= (1L << d);
        }
        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 is 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 (dimension >= Long.SIZE) ? -1 : (1L << dimension) - 1;
            }
        }
        return nonLinearDimensions;
    }

    /**
     * Invoked when the number of non-linear dimensions exceeds the {@code GridGeometry} capacity.
     */
    private 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 parameters 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 attributes are defined (i.e. invoking the
     *         corresponding method will not thrown an {@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 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.
     *
     * <div class="note"><b>Example:</b>
     * for clipping this grid geometry to a sub-area, one can use:
     *
     * {@preformat java
     *     GridGeometry gg = ...;
     *     Envelope areaOfInterest = ...;
     *     gg = gg.derive().rounding(GridRoundingMode.ENCLOSING)
     *                     .subgrid(areaOfInterest).build();
     * }
     * </div>
     *
     * 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 <cite>slice point</cite>,
     * but does not remove those dimensions from the grid geometry. For dimensionality reduction, see {@link #reduce(int...)}.
     *
     * @return an object for deriving a grid geometry from {@code this}.
     */
    public GridDerivation derive() {
        return new GridDerivation(this);
    }

    /**
     * Returns a grid geometry that encompass only some dimensions of the grid geometry.
     * The specified dimensions will be copied into a new grid geometry.
     * The selection is applied on {@linkplain #getExtent() grid extent} dimensions;
     * they are not necessarily the same than 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 <cite>dimensionality reduction</cite>.
     * This method can not be used for changing dimension order.</p>
     *
     * @param  dimensions  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#reduce(int...)
     * @see org.apache.sis.referencing.CRS#reduce(CoordinateReferenceSystem, int...)
     */
    public GridGeometry reduce(int... dimensions) {
        dimensions = GridExtent.verifyDimensions(dimensions, getDimension());
        if (dimensions != null) try {
            return new GridGeometry(this, dimensions);
        } catch (FactoryException e) {
            throw new IllegalGridGeometryException(e, "dimensions");
        }
        return this;
    }

    /**
     * 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.
     *
     * @param  object  the object to compare with.
     * @return {@code true} if the given object is equals to this grid geometry.
     */
    @Override
    public boolean equals(final Object object) {
        if (object != null && object.getClass() == getClass()) {
            final GridGeometry that = (GridGeometry) object;
            return Objects.equals(this.extent,    that.extent)    &&
                   Objects.equals(this.gridToCRS, that.gridToCRS) &&
                   Objects.equals(this.envelope,  that.envelope);
        }
        return false;
    }

    /**
     * Returns the default set of flags to use for {@link #toString()} implementations.
     * Current implementation returns all flags, but future implementation may omit some
     * flags if experience suggests that they are too verbose in practice.
     */
    static int defaultFlags() {
        return EXTENT | ENVELOPE | CRS | GRID_TO_CRS | RESOLUTION | GEOGRAPHIC_EXTENT | TEMPORAL_EXTENT;
    }

    /**
     * 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 purpose only 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.getResources(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
             *  ├─Upper bound:  90°00′00″N  180°00′00″E  2019-05-02T00:00:00Z
             *  └─Lower bound:  80°00′00″S  180°00′00″W  2019-05-01T21: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 = geographicBBox();
                final Instant[] times = timeRange();
                vocabulary.appendLabel(Vocabulary.Keys.UpperBound, table);
                table.setCellAlignment(TableAppender.ALIGN_RIGHT);
                if (bbox != null) {
                    nf.setRoundingMode(RoundingMode.CEILING);
                    table.nextColumn(); table.append(nf.format(new  Latitude(bbox.getNorthBoundLatitude())));
                    table.nextColumn(); table.append(nf.format(new Longitude(bbox.getEastBoundLongitude())));
                }
                if (times.length >= 2) {
                    table.nextColumn();
                    table.append(times[1].toString());
                }
                table.nextLine();
                table.setCellAlignment(TableAppender.ALIGN_LEFT);
                vocabulary.appendLabel(Vocabulary.Keys.LowerBound, table);
                table.setCellAlignment(TableAppender.ALIGN_RIGHT);
                if (bbox != null) {
                    nf.setRoundingMode(RoundingMode.FLOOR);
                    table.nextColumn(); table.append(nf.format(new  Latitude(bbox.getSouthBoundLatitude())));
                    table.nextColumn(); table.append(nf.format(new Longitude(bbox.getWestBoundLongitude())));
                }
                if (times.length >= 1) {
                    table.nextColumn();
                    table.append(times[0].toString());
                }
                table.flush();
                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));
                    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());
                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);
            }
        }
    }
}
