/*
 * 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.referencing.gazetteer;

import java.util.Arrays;
import java.util.Map;
import java.util.IdentityHashMap;
import java.util.ConcurrentModificationException;
import java.util.Iterator;
import java.util.Spliterator;
import java.util.Spliterators;
import java.util.stream.Stream;
import java.util.function.Consumer;
import java.util.stream.StreamSupport;
import java.awt.geom.Rectangle2D;
import javax.xml.bind.annotation.XmlTransient;
import org.opengis.util.FactoryException;
import org.opengis.geometry.Envelope;
import org.opengis.geometry.DirectPosition;
import org.opengis.referencing.crs.SingleCRS;
import org.opengis.referencing.crs.ProjectedCRS;
import org.opengis.referencing.crs.CoordinateReferenceSystem;
import org.opengis.referencing.operation.MathTransform;
import org.opengis.referencing.operation.MathTransform2D;
import org.opengis.referencing.operation.OperationMethod;
import org.opengis.referencing.operation.CoordinateOperation;
import org.opengis.referencing.operation.Projection;
import org.opengis.referencing.operation.TransformException;
import org.apache.sis.internal.referencing.provider.TransverseMercator;
import org.apache.sis.internal.referencing.provider.PolarStereographicA;
import org.apache.sis.internal.gazetteer.Resources;
import org.apache.sis.referencing.CRS;
import org.apache.sis.referencing.CommonCRS;
import org.apache.sis.referencing.NamedIdentifier;
import org.apache.sis.referencing.IdentifiedObjects;
import org.apache.sis.referencing.cs.AxesConvention;
import org.apache.sis.referencing.crs.DefaultProjectedCRS;
import org.apache.sis.internal.referencing.j2d.IntervalRectangle;
import org.apache.sis.math.MathFunctions;
import org.apache.sis.util.CharSequences;
import org.apache.sis.util.ArgumentChecks;
import org.apache.sis.util.StringBuilders;
import org.apache.sis.util.Utilities;
import org.apache.sis.util.Workaround;
import org.apache.sis.util.resources.Errors;
import org.apache.sis.util.resources.Vocabulary;
import org.apache.sis.geometry.Shapes2D;
import org.apache.sis.geometry.Envelopes;
import org.apache.sis.geometry.Envelope2D;
import org.apache.sis.geometry.DirectPosition2D;
import org.apache.sis.internal.util.Strings;
import org.apache.sis.math.DecimalFunctions;
import org.apache.sis.measure.Longitude;
import org.apache.sis.measure.Latitude;

// Branch-dependent imports
import org.apache.sis.metadata.iso.citation.AbstractParty;


/**
 * The Military Grid Reference System (MGRS).
 * The MGRS is the geocoordinate standard used by NATO militaries for locating points on the earth.
 * It is based on the Universal Transverse Mercator (UTM) and the Universal Polar Stereographic (UPS) projections.
 * Despite its name, MGRS is used not only for military purposes; it is used also for organizing Earth Observation
 * data in directory trees for example.
 *
 * <p>MGRS references are sequences of digits and letters like “4Q FJ 12345 67890” (a reference with 1 metre accuracy),
 * optionally written with reduced resolution as in “4Q FJ 123 678” (a reference with 100 metres accuracy).
 * Those references form a hierarchy of 3 {@linkplain ModifiableLocationType location types}:</p>
 *
 * <blockquote>
 *   <b>Grid zone designator</b>         (example: “4Q”)<br>
 *     └─<b>100 km square identifier</b> (example: “FJ”)<br>
 *         └─<b>Grid coordinate</b>      (example: “12345 67890”)<br>
 * </blockquote>
 *
 * <p>Conversions between MGRS references and spatial coordinates can be performed by the {@link Coder Coder} inner class.
 * The result of decoding a MGRS reference is an envelope rather than a point, but a representative point can be obtained.
 * The encoding and decoding processes take in account Norway and Svalbard special cases (they have wider UTM zones for
 * historical reasons).</p>
 *
 * <div class="note"><b>Example:</b>
 * the following code:
 *
 * {@preformat java
 *     MilitaryGridReferenceSystem system = new MilitaryGridReferenceSystem();
 *     MilitaryGridReferenceSystem.Coder coder = system.createCoder();
 *     Location loc = coder.decode("32TNL83");
 *     System.out.println(loc);
 * }
 *
 * should display (locale may vary):
 *
 * {@preformat text
 *     ┌─────────────────────────────────────────────────────────────────┐
 *     │ Location type:               Grid coordinate                    │
 *     │ Geographic identifier:       32TNL83                            │
 *     │ West bound:                    580,000 m    —     9°57′00″E     │
 *     │ Representative value:          585,000 m    —    10°00′36″E     │
 *     │ East bound:                    590,000 m    —    10°04′13″E     │
 *     │ South bound:                 4,530,000 m    —    40°54′58″N     │
 *     │ Representative value:        4,535,000 m    —    40°57′42″N     │
 *     │ North bound:                 4,540,000 m    —    41°00′27″N     │
 *     │ Coordinate reference system: WGS 84 / UTM zone 32N              │
 *     │ Administrator:               North Atlantic Treaty Organization │
 *     └─────────────────────────────────────────────────────────────────┘
 * }
 * </div>
 *
 * <h2>Immutability and thread safety</h2>
 * This class is immutable and thus thread-safe.
 * However the {@link Coder Coder} instances performing conversions between references and coordinates
 * are not thread-safe; it is recommended to create a new {@code Coder} instance for each thread.
 *
 * @author  Martin Desruisseaux (Geomatys)
 * @version 0.8
 *
 * @see CommonCRS#universal(double, double)
 * @see <a href="https://en.wikipedia.org/wiki/Military_Grid_Reference_System">Military Grid Reference System on Wikipedia</a>
 *
 * @since 0.8
 * @module
 */
@XmlTransient
public class MilitaryGridReferenceSystem extends ReferencingByIdentifiers {
    /**
     * For cross-version compatibility.
     */
    private static final long serialVersionUID = 8337394374656125471L;

    /**
     * Height of latitude bands, in degrees.
     * Those bands are labeled from {@code 'C'} to {@code 'X'} inclusive, excluding {@code 'I'} and {@code 'O'}.
     */
    static final double LATITUDE_BAND_HEIGHT = 8;

    /**
     * Size of the 100 kilometres squares, in metres.
     */
    static final double GRID_SQUARE_SIZE = 100_000;

    /**
     * Number of letters available for grid rows. Those letters are "ABCDEFGHJKLMNPQRSTUV" (starting at letter
     * F for zones of even number), repeated in a cycle. Each row is {@value #GRID_SQUARE_SIZE} metres height.
     */
    static final int GRID_ROW_COUNT = 20;

    /**
     * The number of digits in a one-meter precision when formatting MGRS references.
     *
     * <p><b>Invariant:</b> the following relationship must hold:
     * {@code GRID_SQUARE_SIZE == Math.pow(10, METRE_PRECISION_DIGITS)}
     */
    static final int METRE_PRECISION_DIGITS = 5;

    /**
     * The first of the two letters ({@code 'I'} and {@code 'O'}) excluded in MGRS notation.
     * This letter and all following letters shall be shifted by one character. Example:
     *
     * {@preformat java
     *     char band = ...;
     *     if (band >= EXCLUDE_I) {
     *         band++;
     *         if (band >= EXCLUDE_O) band++;
     *     }
     * }
     *
     * or equivalently:
     *
     * {@preformat java
     *     char band = ...;
     *     if (band >= EXCLUDE_I && ++band >= EXCLUDE_O) band++;
     * }
     */
    private static final char EXCLUDE_I = 'I';

    /**
     * The second of the two letters ({@code 'I'} and {@code 'O'}) excluded in MGRS notation.
     */
    private static final char EXCLUDE_O = 'O';

    /**
     * The column letters used in Polar Stereographic Projections.
     * They are letters A to Z but omitting I, O, D, E, M, N, V, W.
     */
    private static final byte[] POLAR_COLUMNS = {
        'A','B','C','F','G','H','J','K','L','P','Q','R','S','T','U','X','Y','Z'
    };

    /**
     * The object to use for computing zone number and central meridian. This is a static final
     * field in current Apache SIS version, but could become configurable in a future version.
     */
    private static final TransverseMercator.Zoner ZONER = TransverseMercator.Zoner.UTM;

    /**
     * The datum to which to transform the coordinate before formatting the MGRS reference.
     * Only the datums enumerated in {@link CommonCRS} are currently supported.
     */
    final CommonCRS datum;

    /**
     * Whether {@link MilitaryGridReferenceSystem.Encoder} should infer the datum from the given coordinates
     * instead than using {@link #datum}.
     */
    final boolean avoidDatumChange;

    /**
     * Value to add to the row number in order to have the "A" letter on the northernmost value on Greenwich meridian of
     * the Universal Polar Stereographic (UPS) South projection. Value is initially zero and computed when first needed.
     * This is derived from the bottom of the 100 kilometres square labeled "A" in Grid Zone Designations A and B.
     */
    private transient short southOffset;

    /**
     * Value to add to the row number in order to have the "A" letter on the southernmost value on Greenwich meridian of
     * the Universal Polar Stereographic (UPS) North projection. Value is initially zero and computed when first needed.
     * This is derived from the bottom of the 100 kilometres square labeled "A" in Grid Zone Designations Y and Z.
     */
    private transient short northOffset;

    /**
     * Creates a new Military Grid Reference System (MGRS) using the default datum.
     * The current Apache SIS version uses the {@linkplain CommonCRS#WGS84 WGS84} datum,
     * but this choice may change in the future if there is a need to adapt to new MGRS specifications.
     */
    public MilitaryGridReferenceSystem() {
        super(properties(), types());
        datum = CommonCRS.WGS84;
        avoidDatumChange = false;
    }

    /**
     * Creates a new Military Grid Reference System (MGRS) using the specified datum.
     * Only the datums enumerated in {@link CommonCRS} are currently supported.
     *
     * @param  properties  the properties to be given to the reference system.
     * @param  datum       the datum to which to transform coordinates before formatting the MGRS references,
     *                     or {@code null} for inferring the datum from the CRS associated to each coordinate.
     */
    public MilitaryGridReferenceSystem(final Map<String,?> properties, final CommonCRS datum) {
        super(properties, types());
        this.datum = (datum != null) ? datum : CommonCRS.WGS84;
        avoidDatumChange = (datum == null);
    }

    /**
     * Work around for RFE #4093999 in Sun's bug database
     * ("Relax constraint on placement of this()/super() call in constructors").
     */
    @Workaround(library="JDK", version="1.8")
    private static Map<String,?> properties() {
        AbstractParty party = new AbstractParty("North Atlantic Treaty Organization", null);
        return properties(new NamedIdentifier(null, "NATO", Resources.formatInternational(Resources.Keys.MGRS), null, null), party);
    }

    /**
     * Work around for RFE #4093999 in Sun's bug database
     * ("Relax constraint on placement of this()/super() call in constructors").
     */
    @Workaround(library="JDK", version="1.8")
    private static ModifiableLocationType[] types() {
        final ModifiableLocationType gzd    = new ModifiableLocationType(Resources.formatInternational(Resources.Keys.GridZoneDesignator));
        final ModifiableLocationType square = new ModifiableLocationType(Resources.formatInternational(Resources.Keys.SquareIdentifier100));
        final ModifiableLocationType coord  = new ModifiableLocationType(Resources.formatInternational(Resources.Keys.GridCoordinates));
        gzd   .addIdentification(Vocabulary.formatInternational(Vocabulary.Keys.Code));
        coord .addIdentification(Vocabulary.formatInternational(Vocabulary.Keys.Coordinate));
        square.addParent(gzd);
        coord .addParent(square);
        return new ModifiableLocationType[] {gzd};
    }

    /**
     * Returns the value to add to the row number in order to have the "A" letter on the southernmost or
     * northernmost value on Greenwich meridian of the Universal Polar Stereographic (UPS) projection.
     * If {@code south} is {@code true}, then this is computed from the northernmost value of UPS South;
     * otherwise this is computed from the southernmost value of UPS North. This value is derived from
     * the bottom of the 100 kilometres square labeled "A" in Grid Zone Designations A, B, Y and Z.
     */
    final int polarOffset(final boolean south) throws TransformException {
        // No need to synchronized; not a big deal if computed twice.
        short origin = south ? southOffset : northOffset;
        if (origin == 0) {
            final DirectPosition2D position = new DirectPosition2D(
                    south ? TransverseMercator.Zoner.SOUTH_BOUNDS
                          : TransverseMercator.Zoner.NORTH_BOUNDS, 0);
            double northing = datum.universal(position.x * 1.01, position.y).getConversionFromBase()
                                   .getMathTransform().transform(position, position).getOrdinate(1);
            if (south) {
                northing = 2*PolarStereographicA.UPS_SHIFT - northing;
            }
            northing = Math.floor(northing / GRID_SQUARE_SIZE);
            origin = (short) northing;
            if (origin != northing) {                       // Paranoiac check (should never happen).
                throw new GazetteerException();
            }
            if (south) {
                southOffset = origin;
            } else {
                northOffset = origin;
            }
        }
        return origin;
    }

    /**
     * Returns a new object performing conversions between {@code DirectPosition} and MGRS references.
     * The returned object is <strong>not</strong> thread-safe; a new instance must be created for
     * each thread, or synchronization must be applied by the caller.
     *
     * @return a new object performing conversions between {@link DirectPosition} and MGRS references.
     */
    public Coder createCoder() {
        return new Coder();
    }

    /**
     * Conversions between direct positions and references in the Military Grid Reference System (MGRS).
     * Each {@code Coder} instance can read references at arbitrary precision, but formats at the
     * {@linkplain #setPrecision specified precision}.
     * The same {@code Coder} instance can be reused for reading or writing many MGRS references.
     *
     * <p>See the {@link MilitaryGridReferenceSystem} enclosing class for usage example.</p>
     *
     * <h2>Immutability and thread safety</h2>
     * This class is <strong>not</strong> thread-safe. A new instance must be created for each thread,
     * or synchronization must be applied by the caller.
     *
     * @author  Martin Desruisseaux (Geomatys)
     * @version 0.8
     * @since   0.8
     * @module
     */
    public class Coder {
        /**
         * Number of digits to use for formatting the numerical part of a MGRS reference.
         *
         * @see #getPrecision()
         */
        private byte digits;

        /**
         * The separator to insert between each component of the MGRS identifier, or an empty string if none.
         */
        private String separator;

        /**
         * Same separator, but without leading and trailing spaces.
         */
        private String trimmedSeparator;

        /**
         * Whether the decoded locations should be clipped to the valid area of MGRS cell.
         *
         * @see #getClipToValidArea()
         */
        private boolean clipToValidArea;

        /**
         * Cached information needed for building a MGRS reference from a direct position in the given CRS.
         */
        private final Map<CoordinateReferenceSystem,Encoder> encoders;

        /**
         * Temporary positions used for encoding. References are kept for avoiding
         * to recreate those temporary objects for every reference to parse or format.
         */
        transient DirectPosition normalized, geographic;

        /**
         * A buffer where to create reference, to be reused for each new reference.
         */
        final StringBuilder buffer = new StringBuilder(18);      // Length of "4 Q FJ 12345 67890" sample value.

        /**
         * Creates a new coder initialized to the default precision and separator.
         */
        protected Coder() {
            digits    = METRE_PRECISION_DIGITS;     // 1 metre precision.
            separator = trimmedSeparator = "";
            encoders  = new IdentityHashMap<>();
            clipToValidArea = true;
        }

        /**
         * Creates a new coder initialized to the same setting than the given separator.
         * The new instance will share the same {@link #encoders} map than the original instance.
         * This is okay only if all calls to {@link #encoder(CoordinateReferenceSystem)} are done
         * in the same thread before any call to {@link IteratorAllZones#trySplit()}.
         */
        Coder(final Coder other) {
            digits           = other.digits;
            separator        = other.separator;
            trimmedSeparator = other.trimmedSeparator;
            clipToValidArea  = other.clipToValidArea;
            encoders         = other.encoders;
        }

        /**
         * Returns the reference system for which MGRS references will be encoded or decoded.
         *
         * @return the enclosing reference system.
         */
        final MilitaryGridReferenceSystem getReferenceSystem() {
            return MilitaryGridReferenceSystem.this;
        }

        /**
         * Returns the precision of the references formatted by this coder.
         * This method returns one of the following values:
         *
         * <table class="sis">
         *   <caption>MGRS reference precisions</caption>
         *   <tr><th>Precision (m)</th>             <th>Reference example</th></tr>
         *   <tr><td style="text-align:right">1</td> <td>4 Q FJ 12345 67890</td></tr>
         *   <tr><td style="text-align:right">10</td> <td>4 Q FJ 1234 6789</td></tr>
         *   <tr><td style="text-align:right">100</td> <td>4 Q FJ 123 678</td></tr>
         *   <tr><td style="text-align:right">1000</td> <td>4 Q FJ 12 67</td></tr>
         *   <tr><td style="text-align:right">10 000</td> <td>4 Q FJ 1 6</td></tr>
         *   <tr><td style="text-align:right">100 000</td> <td>4 Q FJ</td></tr>
         *   <tr><td style="text-align:right">(approximated) 1 000 000</td> <td>4 Q</td></tr>
         * </table>
         *
         * Values smaller than 1 (e.g. 0.01 for a centimetre precision) may also be returned
         * if that value has been {@linkplain #setPrecision(double) explicitly set},
         * but sub-metric precision are usually not used with MGRS.
         *
         * @return precision of formatted references in metres.
         */
        public double getPrecision() {
            return MathFunctions.pow10(METRE_PRECISION_DIGITS - digits);
        }

        /**
         * Sets the desired precision of the references formatted by this coder.
         * This method rounds the given precision to one of the power of 10
         * documented in the {@link #getPrecision()} method.
         *
         * @param  precision  the desired precision in metres.
         */
        public void setPrecision(final double precision) {
            final int p;
            try {
                p = DecimalFunctions.floorLog10(precision);
            } catch (ArithmeticException e) {
                throw new IllegalArgumentException(Errors.format(Errors.Keys.IllegalArgumentValue_2, "precision", precision), e);
            }
            // The -3 is an arbitrary limit to millimetre precision.
            int n = Math.max(-3, Math.min(METRE_PRECISION_DIGITS + 1, p));
            digits = (byte) (METRE_PRECISION_DIGITS - n);
        }

        /**
         * For internal use by other internal classes in this Java source file.
         */
        final int digits() {
            return digits;
        }

        /**
         * Returns the separator to insert between each component of the MGRS identifier.
         * Components are zone number, latitude band, 100 000-metres square identifier and numerical values.
         * By default the separator is an empty string, which produce references like "4QFJ12345678".
         *
         * @return the separator to insert between each component of the MGRS identifier, or an empty string if none.
         */
        public String getSeparator() {
            return separator;
        }

        /**
         * Sets the separator to insert between each component of the MGRS identifier.
         * Components are zone number, latitude band, 100 000-metres square identifier and numerical values.
         * By default the separator is an empty string, which produce references like "4QFJ12345678".
         * If the separator is set to a space, then the references will be formatted like "4 Q FJ 1234 5678".
         *
         * <p>Note that a MGRS reference is normally written as an entity without spaces, parentheses, dashes,
         * or decimal points. Invoking this method with a non-empty separator produces non-conform MGRS, but
         * is sometime convenient for readability or for use in file systems (with the {@code '/'} separator).</p>
         *
         * @param  separator  the separator to insert between each component of the MGRS identifier.
         */
        public void setSeparator(final String separator) {
            ArgumentChecks.ensureNonNull("separator", separator);
            this.separator = separator;
            trimmedSeparator = CharSequences.trimWhitespaces(separator);
        }

        /**
         * Returns whether the decoded locations should be clipped to the valid area.
         * The default value is {@code true}.
         *
         * @return {@code true} if decoded locations are clipped to the valid area.
         */
        public boolean getClipToValidArea() {
            return clipToValidArea;
        }

        /**
         * Sets whether the decoded locations should be clipped to the valid area.
         * MGRS 100 km squares can actually be smaller than 100 km when the square overlaps two UTM zones or
         * two latitude bands. We may have half of a square in a zone and the other half in the other zone.
         * By default, the {@link #decode(CharSequence)} method clips the square to the zone where it belongs.
         * Invoking this method with the {@code false} value disables this behavior.
         *
         * @param clip  whether the decoded locations should be clipped to the valid area.
         */
        public void setClipToValidArea(final boolean clip) {
            clipToValidArea = clip;
        }

        /**
         * Bridge to {@link MilitaryGridReferenceSystem#datum}
         * for the {@link Encoder} and {@link Decoder} classes.
         */
        final ProjectedCRS projection(final double latitude, final double longitude) {
            return MilitaryGridReferenceSystem.this.datum.universal(latitude, longitude);
        }

        /**
         * Returns the encoder for the given coordinate reference system.
         * All calls to this method must be done in the same thread.
         *
         * @throws IllegalArgumentException if the given CRS do not use one of the supported datums.
         * @throws FactoryException if the creation of a coordinate operation failed.
         * @throws TransformException if the creation of an inverse operation failed.
         */
        final Encoder encoder(final CoordinateReferenceSystem crs) throws FactoryException, TransformException {
            if (crs == null) {
                throw new GazetteerException(Errors.format(Errors.Keys.UnspecifiedCRS));
            }
            // We can not use encoders.computeIfAbsent(crs, ...) because of checked exceptions.
            Encoder encoder = encoders.get(crs);
            if (encoder == null) {
                encoder = new Encoder(avoidDatumChange ? null : datum, crs);
                if (encoders.put(crs, encoder) != null) {
                    throw new ConcurrentModificationException();            // Opportunistic check.
                }
            }
            return encoder;
        }

        /**
         * Encodes the given position into a MGRS reference.
         * The given position must have a Coordinate Reference System (CRS) associated to it.
         *
         * @param  position  the coordinate to encode.
         * @return MGRS encoding of the given position.
         * @throws TransformException if an error occurred while transforming the given coordinate to a MGRS reference.
         */
        public String encode(final DirectPosition position) throws TransformException {
            ArgumentChecks.ensureNonNull("position", position);
            try {
                return encoder(position.getCoordinateReferenceSystem()).encode(this, position, true, getSeparator(), digits());
            } catch (IllegalArgumentException | FactoryException e) {
                throw new GazetteerException(e.getLocalizedMessage(), e);
            }
        }

        /**
         * Returns an iterator over all MGRS references that intersect the given envelope.
         * The given envelope must have a Coordinate Reference System (CRS) associated to it.
         * If the CRS is geographic, the envelope is allowed to span the anti-meridian.
         * The MGRS references may be returned in any iteration order.
         *
         * <div class="note"><b>Possible evolution:</b>
         * current implementation does not clip the cells to UPS/UTM valid areas before to test for intersection
         * with {@code areaOfInterest}. Consequently the iterator may return slightly more cells than expected.
         * A future version may filter the cells more accurately. If an application needs the same set of cells
         * than what current the implementation returns, it can invoke <code>{@linkplain #setClipToValidArea
         * setClipToValidArea}(false)</code> for preserving current behavior in future Apache SIS versions.</div>
         *
         * @param  areaOfInterest  envelope of desired MGRS references.
         * @return an iterator over MGRS references intersecting the given area of interest.
         * @throws TransformException if an error occurred while transforming the area of interest.
         */
        public Iterator<String> encode(final Envelope areaOfInterest) throws TransformException {
            ArgumentChecks.ensureNonNull("areaOfInterest", areaOfInterest);
            try {
                return Spliterators.iterator(new IteratorAllZones(areaOfInterest).simplify());
            } catch (IllegalArgumentException | FactoryException e) {
                throw new GazetteerException(e);
            }
        }

        /**
         * Returns a stream of all MGRS references that intersect the given envelope.
         * The given envelope must have a Coordinate Reference System (CRS) associated to it.
         * If the CRS is geographic, the envelope is allowed to span the anti-meridian.
         * The MGRS references may be returned in any order.
         *
         * <div class="note"><b>Possible evolution:</b>
         * current implementation does not clip the cells to UPS/UTM valid areas before to test for intersection
         * with {@code areaOfInterest}. Consequently the iterator may return slightly more cells than expected.
         * A future version may filter the cells more accurately. If an application needs the same set of cells
         * than what current the implementation returns, it can invoke <code>{@linkplain #setClipToValidArea
         * setClipToValidArea}(false)</code> for preserving current behavior in future Apache SIS versions.</div>
         *
         * @param  areaOfInterest  envelope of desired MGRS references.
         * @param  parallel        {@code true} for a parallel stream, or {@code false} for a sequential stream.
         * @return a stream of MGRS references intersecting the given area of interest.
         * @throws TransformException if an error occurred while transforming the area of interest.
         */
        public Stream<String> encode(final Envelope areaOfInterest, final boolean parallel) throws TransformException {
            ArgumentChecks.ensureNonNull("areaOfInterest", areaOfInterest);
            try {
                return StreamSupport.stream(new IteratorAllZones(areaOfInterest).simplify(), parallel);
            } catch (IllegalArgumentException | FactoryException e) {
                throw new GazetteerException(e);
            }
        }

        /**
         * Decodes the given MGRS reference into a position and an envelope.
         * The Coordinate Reference System (CRS) associated to the returned position depends on the given reference.
         *
         * <div class="warning"><b>Upcoming API change — generalization</b><br>
         * in a future SIS version, the type of returned element may be generalized
         * to the {@code org.opengis.referencing.gazetteer.Location} interface.
         * This change is pending GeoAPI revision.</div>
         *
         * @param  reference  MGRS string to decode.
         * @return a new position with the longitude at coordinate 0 and latitude at coordinate 1.
         * @throws TransformException if an error occurred while parsing the given string.
         */
        public AbstractLocation decode(final CharSequence reference) throws TransformException {
            ArgumentChecks.ensureNonEmpty("reference", reference);
            return new Decoder(this, reference);
        }


        /**
         * Iterator over the cells inside all UPS and UTM zones inside a given area of interest.
         * Each UPS or UTM zone is processed by a separated iterator, each of them with its own
         * {@link Encoder} instance.
         *
         * @see IteratorOneZone
         */
        private final class IteratorAllZones implements Spliterator<String> {
            /**
             * The iterators over a single UTM zone.
             */
            private final Spliterator<String>[] iterators;

            /**
             * Index of the current iterator.
             */
            private int index;

            /**
             * Index after the last iterator to return.
             */
            private int upper;

            /**
             * Creates a new iterator over MGRS cells in the given area of interest.
             * The borders of the given envelope are considered <strong>exclusive</strong>.
             */
            @SuppressWarnings({"unchecked", "rawtypes"})        // Generic array creation.
            IteratorAllZones(final Envelope areaOfInterest) throws FactoryException, TransformException {
                final SingleCRS sourceCRS = CRS.getHorizontalComponent(areaOfInterest.getCoordinateReferenceSystem());
                if (sourceCRS == null) {
                    throw new GazetteerException(org.apache.sis.internal.referencing.Resources.format(
                            org.apache.sis.internal.referencing.Resources.Keys.NonHorizontalCRS_1, "areaOfInterest"));
                }
                final int precision = (int) getPrecision();
                if (precision <= 0 || precision > (int) GRID_SQUARE_SIZE) {
                    throw new GazetteerException(Errors.format(Errors.Keys.ValueOutOfRange_4,
                            "precision", 1, (int) GRID_SQUARE_SIZE, precision));
                }
                /*
                 * Convert area of interest (AOI) from an envelope to a Rectangle2D for use with
                 * Envelope2D.intersect(Rectangle2D) during IteratorOneZone.advance(…) execution.
                 * We need to use the constructor expecting the two corners in order to preserve
                 * envelope spanning the anti-meridian.
                 */
                final IntervalRectangle aoi = new IntervalRectangle(areaOfInterest.getLowerCorner(),
                                                                    areaOfInterest.getUpperCorner());
                /*
                 * Project the area of interest (AOI) to normalized geographic coordinates for computing UTM zones
                 * and for IteratorOneZone construction. For computing UTM zones, envelopes that cross the anti-
                 * meridian should have a negative width. For IteratorOneZone construction, it does not matter.
                 */
                final Envelope geographicArea = Envelopes.transform(areaOfInterest, datum.normalizedGeographic());
                final double φmin   = geographicArea.getMinimum(1);
                final double φmax   = geographicArea.getMaximum(1);
                boolean   southPole = (φmin < TransverseMercator.Zoner.SOUTH_BOUNDS);
                boolean   northPole = (φmax > TransverseMercator.Zoner.NORTH_BOUNDS);
                boolean   southUTM  = false;
                boolean   northUTM  = false;                    // Default value for UPS, to be modified later if UTM.
                int       zoneStart = 0;                        // Zero-based (i.e. standard zone number minus 1).
                int       zoneEnd   = 0;                        // Exclusive.
                final int zoneCount = ZONER.zoneCount();
                if (φmax >= TransverseMercator.Zoner.SOUTH_BOUNDS && φmin < TransverseMercator.Zoner.NORTH_BOUNDS) {
                    southUTM = (φmin <  0);
                    northUTM = (φmax >= 0);
                    if (geographicArea.getSpan(0) > (Longitude.MAX_VALUE - Longitude.MIN_VALUE) - ZONER.width) {
                        zoneEnd = zoneCount;
                    } else {
                        /*
                         * Use of lower and upper corners below are not the same than calls to Envelope.getMinimum(0)
                         * or Envelope.getMaximum(0) if the envelope crosses the anti-meridian.
                         */
                        zoneStart = ZONER.zone(0, geographicArea.getLowerCorner().getOrdinate(0)) - 1;  // Inclusive.
                        zoneEnd   = ZONER.zone(0, geographicArea.getUpperCorner().getOrdinate(0));      // Exclusive.
                        if (zoneEnd < zoneStart) {
                            zoneEnd += zoneCount;                              // Envelope crosses the anti-meridian.
                        }
                    }
                }
                /*
                 * At this point we finished collecting the information that we will need (whether we will
                 * iterate over the north or south pole, the range of UTM zones, etc). For each UPS or UTM
                 * zone, we pick a representative φ value for the purpose of CommonCRS.universal(…) method
                 * call and we prepare an iterator on which we will delegate computations.
                 */
                upper = zoneEnd - zoneStart;
                if (southUTM & northUTM) upper *= 2;
                if (southPole)           upper += 2;
                if (northPole)           upper += 2;
                iterators = new Spliterator[upper];
                int zone = zoneStart;
                upper = 0;
                do {
                    final double λ = ZONER.centralMeridian((zone % zoneCount) + 1);
                    final double φ;
                    if (southPole) {
                        φ = Latitude.MIN_VALUE;
                        southPole = false;              // For next iteration.
                    } else if (southUTM) {
                        φ = -1;
                        if (++zone >= zoneEnd) {
                            zone = zoneStart;
                            southUTM = false;           // For next iteration.
                        }
                    } else if (northUTM) {
                        φ = +1;
                        if (++zone >= zoneEnd) {
                            northUTM = false;           // For next iteration.
                        }
                    } else if (northPole) {
                        φ = Latitude.MAX_VALUE;
                        northPole = false;              // For loop termination.
                    } else {
                        throw new AssertionError();
                    }
                    final ProjectedCRS targetCRS = datum.universal(φ, λ);
                    Spliterator<String> iter = new IteratorOneZone(Coder.this, aoi, geographicArea, sourceCRS, targetCRS, precision);
                    do iterators[upper++] = iter;
                    while ((iter = iter.trySplit()) != null);
                } while (southPole | northPole | southUTM | northUTM);
            }

            /**
             * Creates an iterator over the first half of the zones covered by the given iterator.
             * After construction, the given iterator will cover the second half. This constructor
             * is for {@link #trySplit()} method only.
             */
            private IteratorAllZones(final IteratorAllZones half) {
                iterators  =  half.iterators;
                index      =  half.index;
                upper      = (half.upper + index) / 2;
                half.index = upper;
            }

            /**
             * If this iterator can be partitioned, returns an iterator covering approximately
             * the first half of MGRS references and update this iterator for covering the other
             * half. Each iterator will use a disjoint set of projected CRS.
             */
            @Override
            public Spliterator<String> trySplit() {
                return (upper - index >= 2) ? new IteratorAllZones(this).simplify() : null;
            }

            /**
             * If this iterator is backed by only one worker iterator, returns that worker iterator.
             * Otherwise returns {@code this}. This method should be invoked after construction.
             */
            final Spliterator<String> simplify() {
                return (upper - index == 1) ? iterators[index] : this;
            }

            /**
             * Guess the number of elements to be returned. The value returned by this method is very rough,
             * and likely greater than the real amount of elements that will actually be returned.
             */
            @Override
            public long estimateSize() {
                long n = 0;
                for (int i=index; i<upper; i++) {
                    n += iterators[i].estimateSize();
                }
                return n;
            }

            /**
             * Performs the given action on the remaining MGRS reference, if any.
             */
            @Override
            public boolean tryAdvance(final Consumer<? super String> action) {
                while (index < upper) {
                    if (iterators[index].tryAdvance(action)) {
                        return true;
                    }
                    index++;
                }
                return false;
            }

            /**
             * Performs the given action on all remaining MGRS references.
             */
            @Override
            public void forEachRemaining(final Consumer<? super String> action) {
                while (index < upper) {
                    iterators[index++].forEachRemaining(action);
                }
            }

            /**
             * Specifies that the list of elements is immutable, that all elements will be distinct
             * and that this iterator never return {@code null} element.
             */
            @Override
            public int characteristics() {
                return IMMUTABLE | DISTINCT | NONNULL;
            }
        }
    }

    /**
     * Iterator over the cells in a single UTM or UPS zone. There is exactly one {@code IteratorOneZone} instance
     * for each Universal Polar Stereographic (UPS) or Universal Transverse Mercator (UTM) projection covered by
     * the area of interest. A given {@code IteratorOneZone} instance use the same projection for all cells.
     *
     * <p>This class extends {@link Coder} in order to freeze the configuration (separator, precision, <i>etc</i>)
     * to the values they have at iterator creation time, and because if we parallelize the iteration, each iterator
     * will need its own {@link Coder#buffer}, {@link Coder#normalized} and {@link Coder#geographic} cache.</p>
     *
     * @see Coder.IteratorAllZones
     */
    private final class IteratorOneZone extends Coder implements Spliterator<String> {
        /**
         * The region for which to return MGRS codes. This envelope can be in any CRS.
         * This shape shall not be modified since the same instance will be shared by
         * many {@code IteratorOneZone}s.
         */
        private final Rectangle2D areaOfInterest;

        /**
         * The transform from the CRS of {@link #encoder} to the CRS of {@link #areaOfInterest}.
         */
        private final MathTransform2D gridToAOI;

        /**
         * The encoder to use for creating MGRS codes.
         */
        private final Encoder encoder;

        /**
         * The easting value of the projection natural origin. This information is used for determining
         * whether {@link #gridX} is on the left side or on right side of the center of UPS or UTM zone.
         * We use this information for choosing the cell corner which is closest to map projection origin.
         */
        private final int xCenter;

        /**
         * The {@link #gridX} value where to stop iteration, exclusive. This value is always greater than
         * {@code gridX} (unless the iteration finished), but is not necessarily greater than {@link #xCenter}.
         */
        private final int xEnd;

        /**
         * The first <var>northing</var> value to use in iteration.
         * The {@link #gridY} value will need to be reset to this value for each new column.
         */
        private int yStart;

        /**
         * The {@link #gridY} values where to stop iteration, exclusive. Invariants:
         * <ul>
         *   <li>{@code gridY < yEnd} shall be true in the North hemisphere, and
         *       {@code gridY > yEnd} shall be true in the South hemisphere.</li>
         * </ul>
         *
         * The reason of {@code gridY} relationship dependency to the hemisphere is because
         * we try to iterate from equator to the pole.
         */
        private int yEnd;

        /**
         * Position of the next MGRS reference to encode. Position is composed of the latitude band number,
         * the row and column indices of current 100 km square, and finally the grid coordinates inside the
         * current 100 km square. All those components are inferred from the (easting, northing) values.
         */
        private int gridX, gridY;

        /**
         * The amount of metres to add to {@code gridX} and to add or subtract to {@code gridY} during iteration.
         * The sign to use when updating the {@code gridY} value depends on whether we are in the North or South
         * hemisphere.
         */
        private final int step;

        /**
         * Whether this iterator should iterates downward over rows. If {@code true}, then {@link #step} shall be
         * subtracted to {@link #gridY} instead of added. We iterate downward in UTM south zones in order to go
         * from equator to pole. This direction allows some optimizations.
         */
        private final boolean downward;

        /**
         * Whether this iterator is allowed to skip some cells when testing for inclusion in the area of interest.
         * Since {@code IteratorOneZone} iterates in UTM zone from equator to pole, the range of longitude values
         * will only decrease (the minimal longitude increase and the maximal longitude decrease). Consequently
         * if we found a longitude out of range, we don't need to test that longitude again in next row.
         *
         * <p>This optimization is allowed only under the following conditions:</p>
         * <ul>
         *   <li>{@link #areaOfInterest} is a rectangle in geographic coordinates.</li>
         *   <li>{@link #areaOfInterest} does not intersect Norway and Svalbard special cases.</li>
         *   <li>We iterate in a UTM zone, or in a UPS zone contained fully in the upper half or fully in lower half.</li>
         * </ul>
         */
        private final boolean optimize;

        /**
         * Latitude band of the last MGRS reference encoded by the iterator.
         * This information is used for detecting when we moved to a new latitude band.
         */
        private char latitudeBand;

        /**
         * A MGRS reference which was pending return by {@link #tryAdvance(Consumer)} before to continue iteration.
         * This field may be non-null immediately after a change of latitude band, and should be null otherwise.
         */
        private String pending;

        /**
         * Temporary rectangle for computation purpose. Needs to be an implementation from the
         * {@link org.apache.sis.geometry} in order to support AOI spanning the anti-meridian.
         */
        private final Envelope2D cell = new Envelope2D();

        /**
         * Returns a new iterator for creating MGRS codes in a single UTM or UPS zone.
         * The borders of the {@code areaOfInterest} rectangle are considered <strong>exclusive</strong>.
         *
         * <p>For envelopes that cross the anti-meridian, it does not matter if {@code geographicArea} uses
         * the negative width convention or is expanded to the [-180 … 180]° of longitude range, because it
         * will be clipped to the projection domain of validity anyway. However the {@code areaOfInterest}
         * should use the negative width convention.</p>
         *
         * @param areaOfInterest  the envelope for which to return MGRS codes. This envelope can be in any CRS.
         * @param geographicArea  the area of interest transformed into a normalized geographic CRS.
         * @param sourceCRS       the horizontal part of the {@code areaOfInterest} CRS.
         * @param targetCRS       the UTM or UPS projected CRS of the zone for which to create MGRS references.
         * @param step            the amount of metres to add or subtract to grid coordinates during iteration.
         */
        IteratorOneZone(final Coder coder, Rectangle2D areaOfInterest, final Envelope geographicArea,
                        final SingleCRS sourceCRS, final ProjectedCRS targetCRS, final int step)
                        throws FactoryException, TransformException
        {
            super(coder);
            this.areaOfInterest = areaOfInterest;
            this.encoder        = encoder(targetCRS);
            this.step           = step;
            /*
             * Compute the geographic bounds of the UPS or UTM zone of validity, together with a representative point
             * (φ,λ₀). We will need to clip the area of interest to those bounds before to project that area, because
             * the UPS and UTM projections can not cover the whole world.
             */
            double λmin, λmax, φmin, φmax;
            final int zone = Math.abs(encoder.crsZone);
            if (zone == Encoder.POLE) {
                xCenter = PolarStereographicA.UPS_SHIFT;
                if (encoder.crsZone < 0) {
                    φmin = Latitude.MIN_VALUE;
                    φmax = TransverseMercator.Zoner.SOUTH_BOUNDS;
                } else {
                    φmin = TransverseMercator.Zoner.NORTH_BOUNDS;
                    φmax = Latitude.MAX_VALUE;
                }
                λmin = Longitude.MIN_VALUE;
                λmax = Longitude.MAX_VALUE;
            } else {
                xCenter = (int) ZONER.easting;
                if (encoder.crsZone < 0) {
                    φmin = TransverseMercator.Zoner.SOUTH_BOUNDS;
                    φmax = 0;
                } else {
                    φmin = 0;
                    φmax = TransverseMercator.Zoner.NORTH_BOUNDS;
                }
                final double λ0 = ZONER.centralMeridian(zone);
                λmin = λ0 - ZONER.width / 2;
                λmax = λ0 + ZONER.width / 2;
            }
            double t;
            boolean clip = false;
            if ((t = geographicArea.getMinimum(1)) >= φmin) φmin = t; else clip = true;
            if ((t = geographicArea.getMaximum(1)) <= φmax) φmax = t; else clip = true;
            if ((t = geographicArea.getMinimum(0)) >= λmin) λmin = t; else clip = true;
            if ((t = geographicArea.getMaximum(0)) <= λmax) λmax = t; else clip = true;
            boolean isSpecialCase = ZONER.isSpecialCase(φmin, φmax, λmin, λmax);
            if (clip) {
                /*
                 * If we detected that the given area of interest is larger than UPS or UTM domain of validity
                 * (in which case above code clipped the geographic bounds), project the CRS domain of validity
                 * to the envelope CRS so we can clip it. Here, "domain of validity" is relative to MGRS grid,
                 * not necessarily to the ISO 19111 domain of validity.
                 */
                final IntervalRectangle r = new IntervalRectangle(λmin, φmin, λmax, φmax);
                r.setRect(Shapes2D.transform(CRS.findOperation(geographicArea.getCoordinateReferenceSystem(), sourceCRS, null), r, r));
                r.intersect(areaOfInterest);
                /*
                 * We need an envelope that do not cross the anti-meridian. If the specified area of interest (AOI)
                 * crosses the anti-meridian (which should happen only with geographic envelopes), then the call to
                 * r.intersect(areaOfInterest) may have overwritten our computation with the AOI longitude range.
                 * Overwrite again that range with the domain of validity of current UTM zone,
                 * based on the following assumptions:
                 *
                 *   1) 'sourceCRS' is geographic (otherwise the AOI should not cross the anti-meridian).
                 *   2) the "normalized" AOI longitude range is [-180 … 180]°, in which case intersection
                 *      with [λmin … λmax] is [λmin … λmax] itself.
                 */
                if (r.xmax < r.xmin) {
                    r.xmin = λmin;
                    r.xmax = λmax;
                }
                areaOfInterest = r;             // Never cross the anti-meridian, even if the original AOI did.
            }
            /*
             * At this point, the area of interest has been clipped to the UPS or UTM domain of validity.
             * Now we can project that area to the CRS managed by this iterator. All values after projection
             * should be positive (because UPS and UTM are designed that way).
             */
            final CoordinateOperation op = CRS.findOperation(sourceCRS, targetCRS, null);
            final Rectangle2D bounds = Shapes2D.transform(op, areaOfInterest, null);
            gridX = (((int)          (bounds.getMinX() / step))) * step;                    // Inclusive
            gridY = (((int)          (bounds.getMinY() / step))) * step;
            xEnd  = (((int) Math.ceil(bounds.getMaxX() / step))) * step;                    // Exclusive
            yEnd  = (((int) Math.ceil(bounds.getMaxY() / step))) * step;
            /*
             * Determine if we should iterate on rows upward or downward. The intent is to iterate from equator to pole
             * in UTM zones, or from projection center to projection border in UPS cases.  Those directions enable some
             * optimizations.
             */
            if (zone != Encoder.POLE) {
                downward = (encoder.crsZone < 0);           // Upward in UTM North zones, downward in UTM South zones.
            } else {
                downward = yEnd <= PolarStereographicA.UPS_SHIFT;  // Downward only if AOI is fully in the lower half.
                /*
                 * In the polar case, we can not apply the shortcut documented in 'optimize' if there is a hole
                 * in the UPS projection center. There is a hole if the latitude of the area of interest does not
                 * reach the pole, or if the longitude range does not make a full circle around the Earth.
                 */
                isSpecialCase        |=        λmin != Longitude.MIN_VALUE ||
                                               λmax != Longitude.MAX_VALUE ||
                        (encoder.crsZone < 0 ? φmin != Latitude .MIN_VALUE
                                             : φmax != Latitude .MAX_VALUE);
            }
            if (downward) {
                final int y = gridY;
                gridY = yEnd - step;
                yEnd  = y    - step;
            }
            yStart    = gridY;
            gridToAOI = (MathTransform2D) op.getMathTransform().inverse();
            /*
             * To be strict, we should also test that the region of interest does not intersect both the upper half
             * and lower half of Universal Polar Stereographic (UPS) projection domain. We do not check that because
             * we need the 'optimize' flag to have the value that we get after 'trySplit()' execution. This implies
             * that invoking 'trySplit()' is mandatory before using this iterator, but IteratorAllZones ensures that.
             */
            optimize = !isSpecialCase && Utilities.equalsIgnoreMetadata(geographicArea.getCoordinateReferenceSystem(), sourceCRS);
        }

        /**
         * Creates an iterator for the lower half of a Universal Polar Stereographic (UPS) projection,
         * and modifies the given iterator for restricting it to the upper half of UPS projection.
         * This method is for {@link #trySplit()} usage only.
         */
        private IteratorOneZone(final IteratorOneZone other) {
            super(other);
            areaOfInterest = other.areaOfInterest;
            gridToAOI      = other.gridToAOI;
            encoder        = other.encoder;
            optimize       = other.optimize;
            step           = other.step;
            gridX          = other.gridX;
            xCenter        = other.xCenter;
            xEnd           = other.xEnd;
            yEnd           = other.yStart - step;                   // Bottom of the zone to iterate, exclusive.
            yStart         = PolarStereographicA.UPS_SHIFT - step;  // Top of the zone to iterate, inclusive.
            other.yStart   = PolarStereographicA.UPS_SHIFT;         // Other iterator will be for the upper half.
            other.gridY    = other.yStart;
            gridY          = yStart;
            downward       = true;
            assert !other.downward;                                 // Fail if the other iterator is not for upper half.
        }

        /**
         * If this iterator intersects both the upper and lower half on UPS domain, returns an iterator for the
         * lower half and modifies this iterator for the upper half. This method <strong>must</strong> be invoked
         * before {@code IteratorOneZone} can be used.
         */
        @Override
        public Spliterator<String> trySplit() {
            if (!downward && Math.abs(encoder.crsZone) == Encoder.POLE && gridY < PolarStereographicA.UPS_SHIFT) {
                return new IteratorOneZone(this);
            }
            return null;
        }

        /**
         * Returns an estimation of the number of cells in the area covered by this iterator. The returned value
         * may be greater than the real amount since we do not take in account the fact that the number of cells
         * in a row become lower as we approach poles.
         */
        @Override
        public long estimateSize() {
            return (xEnd - (long) gridX) * Math.abs(yEnd - (long) yStart) / (step * (long) step);
        }

        /**
         * Computes the next cell reference, if any. This method computes the bounding box in UPS or UTM
         * coordinates, verifies if that box intersects the area of interest, and (if it intersects)
         * delegates to {@link Encoder} for creating the MGRS reference.
         */
        @Override
        public boolean tryAdvance(final Consumer<? super String> action) {
            return advance(action, false);
        }

        /**
         * Performs the given action for each remaining MGRS codes.
         */
        @Override
        public void forEachRemaining(final Consumer<? super String> action) {
            advance(action, true);
            if (pending != null) {
                action.accept(pending);
                pending = null;
            }
        }

        /**
         * Implementation of {@link #tryAdvance(Consumer)} and {@link #forEachRemaining(Consumer)}.
         * The {@code all} argument specifies whether this method is invoked for a single element
         * or for all remaining ones.
         */
        private boolean advance(final Consumer<? super String> action, final boolean all) {
            final int digits = digits();
            final String separator = getSeparator();
            if (normalized == null) {
                normalized = new DirectPosition2D();
            }
            boolean found = false;
            try {
                do {
                    if (pending != null) {
                        action.accept(pending);
                        pending = null;
                        found = true;
                        continue;
                    }
                    /*
                     * Verifies if the current cell should be accepted.
                     * To be accepted, the cell must complies with two conditions:
                     *
                     *   1) It must be inside the area of interest (AOI). Note that the AOI may be in any CRS.
                     *   2) It must be inside the area of validity (AOV). Since the constructor clipped the AOI
                     *      to the area of validity, this test #2 is redundant with test #1 if both AOI and AOV
                     *      use the same CRS. However if AOI and AOV do not use the same CRS, then condition #1
                     *      does not automatically implies condition #2, so we test both.
                     *
                     * Condition #1 is verified by the call to 'areaOfInterest.intersects(…)' below.
                     * Condition #2 is verified indirectly by the call to 'encoder.encode(…)', which return null
                     * if the given point is outside the area of validity.
                     *
                     * Since MGRS references are created by calls to Encoder.encode(…, DirectPosition, …), we need
                     * to select the corner closest to projection center. Otherwise Encoder.encode(…) may consider
                     * that a position is outside the domain of validity while another corner of the same cell would
                     * have been ok. We resolve this issue by shifting 'gridX' toward projection center if 'gridX'
                     * is on the left side.
                     */
                    cell.setRect(gridX, gridY, step, step);
                    cell.setRect(Shapes2D.transform(gridToAOI, cell, cell));
                    if (cell.intersects(areaOfInterest)) {          // Must be invoked on Envelope2D implementation.
                        int x = gridX;
                        int y = gridY;
                        if (x < xCenter) x += step - 1;
                        if (downward)    y += step - 1;
                        normalized.setOrdinate(0, x);
                        normalized.setOrdinate(1, y);
                        String ref = encoder.encode(this, normalized, false, separator, digits);
                        if (ref != null) {
                            /*
                             * If there is a change of latitude band, we may have missed a cell before this one.
                             * We verify that by encoding the cell for the position just before current cell and
                             * comparing the latitude band of the result with our previous latitude band.
                             */
                            char previous = latitudeBand;
                            latitudeBand = encoder.latitudeBand;
                            if (latitudeBand != previous && previous != 0) {
                                pending = ref;
                                normalized.setOrdinate(1, y + (downward ? +1 : -1));
                                ref = encoder.encode(this, normalized, false, separator, digits);
                                if (ref == null || encoder.latitudeBand == previous) {
                                    ref = pending;  // No result or same result than previous iteration - cancel.
                                    pending = null;
                                }
                            }
                            action.accept(ref);
                            found = true;
                        }
                    } else if (optimize) {
                        /*
                         * If this cell is not in the area of interest (AOI) when iterating away from projection
                         * origin (from equator to pole in UTM case), then all other cells after this one in the
                         * same column will not be inside the AOI neither. Set 'gridY' to the limit value so the
                         * condition below will stop the iteration in this column and move to the next column.
                         */
                        gridY = yEnd;
                    }
                    /*
                     * Move to the next cell. We need to do that regardless if the previous block found a cell or not.
                     * We move from equator to pole (UTM case) or projection center to projection border (UPS case)
                     * on the same column, than move one column on the right side when we have reached the last row.
                     */
                    final boolean end = downward ? ((gridY -= step) <= yEnd)        // UTM South or lower part of UPS.
                                                 : ((gridY += step) >= yEnd);       // UTM North or upper part of UPS.
                    if (end) {
                        gridY = yStart;
                        latitudeBand = 0;
                        if ((gridX += step) >= xEnd) {
                            break;
                        }
                    }
                } while (all || !found);
            } catch (FactoryException | TransformException e) {
                // Should never happen since we clipped the area of interest to UPS or UTM domain of valididty.
                throw (ArithmeticException) new ArithmeticException(Errors.format(Errors.Keys.OutsideDomainOfValidity)).initCause(e);
            }
            return found;
        }

        /**
         * Specifies that the list of elements is immutable, that all elements will be distinct
         * and that this iterator never return {@code null} element.
         */
        @Override
        public int characteristics() {
            return IMMUTABLE | DISTINCT | NONNULL;
        }

        /**
         * Returns a string representation of this iterator for debugging purpose.
         */
        @Override
        public String toString() {
            return Strings.toString(getClass(), "zone", encoder.crsZone,
                    "downward", downward, "yStart", yStart, "yEnd", yEnd, "gridX", gridX, "xEnd", xEnd);
        }
    }


    /**
     * Conversions from direct positions to Military Grid Reference System (MGRS) references.
     * Each {@code Encoder} instance is configured for one {@code DirectPosition} CRS.
     * If a position is given in another CRS, another {@code Encoder} instance must be created.
     *
     * <h4>Immutability and thread safety</h4>
     * This class is <strong>not</strong> thread-safe. A new instance must be created for each thread,
     * or synchronization must be applied by the caller.
     *
     * @author  Martin Desruisseaux (Geomatys)
     * @version 0.8
     *
     * @see <a href="https://en.wikipedia.org/wiki/Military_Grid_Reference_System">Military Grid Reference System on Wikipedia</a>
     *
     * @since 0.8
     * @module
     */
    static final class Encoder {
        /**
         * Special {@link #crsZone} value for the UPS (Universal Polar Stereographic) projection.
         * Positive value is used for North pole and negative value for South pole.
         */
        private static final int POLE = 100;

        /**
         * The datum to which to transform the coordinate before formatting the MGRS reference.
         * Only the datums enumerated in {@link CommonCRS} are currently supported.
         */
        private final CommonCRS datum;

        /**
         * UTM zone of position CRS (negative for South hemisphere), or {@value #POLE} (negative of positive)
         * if the CRS is a Universal Polar Stereographic projection, or 0 if the CRS is not a recognized projection.
         * Note that this is not necessarily the same zone than the one to use for formatting any given coordinate
         * in that projected CRS, since the {@link TransverseMercator.Zoner#zone(double, double)} method has special
         * rules for some latitudes.
         */
        final int crsZone;

        /**
         * Coordinate conversion from the position CRS to a CRS of the same type but with normalized axes,
         * or {@code null} if not needed. After conversion, axis directions will be (East, North) and axis
         * units will be metres.
         *
         * <p>This transform should perform only simple operations like swapping axis order and unit conversions.
         * It should not perform more complex operations that would require to go back to geographic coordinates.</p>
         */
        private final MathTransform toNormalized;

        /**
         * Coordinate conversion or transformation from the <em>normalized</em> CRS to geographic CRS.
         * Axis directions are (North, East) as in EPSG geodetic dataset and axis units are degrees.
         * This transform is never {@code null}.
         *
         * <p>This transform may contain datum change from the position datum to the target {@link #datum}.</p>
         */
        private final MathTransform toGeographic;

        /**
         * A transform from a position in the CRS given at construction time to a position in the CRS identified by
         * {@link #actualZone}. This field is updated only when a given position is not located in the zone of the
         * CRS given at construction time.
         */
        private MathTransform toActualZone;

        /**
         * The actual zone where the position to encode is located. Legal values are the same than {@link #crsZone}.
         * If non-zero, then this is the zone of the {@link #toActualZone} transform. This field is updated only when
         * a given position is not located in the zone of the CRS given at construction time.
         */
        private int actualZone;

        /**
         * The latitude band of the last encoded reference.
         * This information is provided for {@link IteratorOneZone} purpose.
         */
        char latitudeBand;

        /**
         * Creates a new converter from direct positions to MGRS references.
         *
         * @param  datum  the datum to which to transform the coordinate before formatting the MGRS reference,
         *                or {@code null} for inferring the datum from the given {@code crs}.
         * @param  crs    the coordinate reference system of the coordinates for which to create MGRS references.
         * @throws IllegalArgumentException if the given CRS do not use one of the supported datums.
         * @throws FactoryException if the creation of a coordinate operation failed.
         * @throws TransformException if the creation of an inverse operation failed.
         */
        Encoder(CommonCRS datum, CoordinateReferenceSystem crs) throws FactoryException, TransformException {
            CoordinateReferenceSystem horizontal = CRS.getHorizontalComponent(crs);
            if (horizontal == null) {
                horizontal = crs;
            }
            if (datum == null) {
                datum = CommonCRS.forDatum(horizontal);
            }
            this.datum = datum;
            if (horizontal instanceof ProjectedCRS) {
                ProjectedCRS  projCRS = (ProjectedCRS) horizontal;
                Projection projection = projCRS.getConversionFromBase();
                final OperationMethod method = projection.getMethod();
                if (IdentifiedObjects.isHeuristicMatchForName(method, TransverseMercator.NAME)) {
                    crsZone = ZONER.zone(projection.getParameterValues());
                } else if (IdentifiedObjects.isHeuristicMatchForName(method, PolarStereographicA.NAME)) {
                    crsZone = POLE * PolarStereographicA.isUPS(projection.getParameterValues());
                } else {
                    crsZone = 0;                                    // Neither UTM or UPS projection.
                }
                if (crsZone != 0) {
                    /*
                     * Usually, the projected CRS already has (E,N) axis orientations with metres units,
                     * so we let 'toNormalized' to null. In the rarer cases where the CRS axes do not
                     * have the expected orientations and units, then we build a normalized version of
                     * that CRS and compute the transformation to that CRS.
                     */
                    final DefaultProjectedCRS userAxisOrder = DefaultProjectedCRS.castOrCopy(projCRS);
                    projCRS = userAxisOrder.forConvention(AxesConvention.NORMALIZED);
                    if (crs != horizontal || projCRS != userAxisOrder) {
                        toNormalized = CRS.findOperation(crs, projCRS, null).getMathTransform();
                        projection   = projCRS.getConversionFromBase();
                        horizontal   = projCRS;
                        crs          = projCRS;         // Next step in the chain of transformations.
                    } else {
                        toNormalized = null;            // ProjectedCRS (UTM or UPS) is already normalized.
                    }
                } else {
                    toNormalized = null;    // ProjectedCRS is neither UTM or UPS — will need full reprojection.
                }
                /*
                 * We will also need the transformation from the normalized projected CRS to latitude and
                 * longitude (in that order) in degrees. We can get this transform directly from the
                 * projected CRS if its base CRS already has the expected axis orientations and units.
                 */
                if (crs == horizontal && Utilities.equalsIgnoreMetadata(projCRS.getBaseCRS(), datum.geographic())) {
                    toGeographic = projection.getMathTransform().inverse();
                    return;
                }
            } else {
                crsZone      = 0;
                toNormalized = null;
            }
            toGeographic = CRS.findOperation(crs, datum.geographic(), null).getMathTransform();
        }

        /**
         * Returns the band letter for the given latitude. It is caller responsibility to ensure that the given latitude
         * is between {@value TransverseMercator.Zoner#SOUTH_BOUNDS} and {@value TransverseMercator.Zoner#NORTH_BOUNDS}
         * inclusive. The returned letter will be one of {@code "CDEFGHJKLMNPQRSTUVWX"} (note that I and O letters are
         * excluded). All bands are 8° height except the X band which is 12° height.
         *
         * @param  φ  the latitude in degrees for which to get the band letter.
         * @return the band letter for the given latitude.
         */
        static char latitudeBand(final double φ) {
            int band = 'C' + (int) ((φ - TransverseMercator.Zoner.SOUTH_BOUNDS) / LATITUDE_BAND_HEIGHT);
            if (band >= EXCLUDE_I && ++band >= EXCLUDE_O && ++band == 'Y') {
                band = 'X';         // Because the last latitude band ('X') is 12° height instead of 8°.
            }
            assert band >= 'C' && band <= 'X' : band;
            return (char) band;
        }

        /**
         * Encodes the given position into a MGRS reference. It is caller responsibility to ensure that
         * the position CRS is the same than the CRS specified at this {@code Encoder} creation time.
         *
         * @param  owner      the {@code Coder} which own this {@code Encoder}.
         * @param  position   the direct position to format as a MGRS reference.
         * @param  reproject  whether this method is allowed to reproject {@code position} when needed.
         * @param  separator  the separator to insert between each component of the MGRS identifier.
         * @param  digits     number of digits to use for formatting the numerical part of a MGRS reference.
         * @return the value of {@code buffer.toString()}, or {@code null} if a reprojection was necessary
         *         but {@code reproject} is {@code false}.
         */
        String encode(final Coder owner, DirectPosition position, final boolean reproject,
                final String separator, final int digits) throws FactoryException, TransformException
        {
            final StringBuilder buffer = owner.buffer;
            if (toNormalized != null) {
                owner.normalized = position = toNormalized.transform(position, owner.normalized);
            }
            final DirectPosition geographic = toGeographic.transform(position, owner.geographic);
            owner.geographic     = geographic;                  // For reuse in next method calls.
            final double  λ      = geographic.getOrdinate(1);
            final double  φ      = geographic.getOrdinate(0);
            final boolean isUTM  = φ >= TransverseMercator.Zoner.SOUTH_BOUNDS &&
                                   φ <  TransverseMercator.Zoner.NORTH_BOUNDS;
            final int zone       = isUTM ? ZONER.zone(φ, λ) : POLE;
            final int signedZone = MathFunctions.isNegative(φ) ? -zone : zone;
            if (signedZone == 0) {
                // Zero value at this point is the result of NaN of infinite coordinate value.
                throw new GazetteerException(Errors.format(Errors.Keys.NotANumber_1, "longitude"));
            }
            /*
             * If the DirectPosition given to this method is not in the expected Coordinate Reference System,
             * transform it now. This may happen because the UTM zone computed above is not the same UTM zone
             * than the coordinate one, or because the coordinate is geographic instead than projected.
             */
            if (signedZone != crsZone) {
                if (!reproject) {
                    return null;
                }
                if (signedZone != actualZone) {
                    actualZone   = 0;                           // In case an exception is thrown on the next line.
                    toActualZone = CRS.findOperation(datum.geographic(), datum.universal(φ, λ), null).getMathTransform();
                    actualZone   = signedZone;
                }
                geographic.setOrdinate(1, Longitude.normalize(λ));
                owner.normalized = position = toActualZone.transform(geographic, owner.normalized);
            }
            /*
             * Grid Zone Designator (GZD).
             */
            buffer.setLength(0);
            if (isUTM) {
                buffer.append(zone).append(separator);
                latitudeBand = latitudeBand(φ);
            } else {
                latitudeBand = (signedZone < 0) ? 'A' : 'Y';
                if (λ >= 0) latitudeBand++;
            }
            buffer.append(latitudeBand);
            /*
             * 100 kilometres square identification.
             */
            if (digits >= 0) {
                final double  x = position.getOrdinate(0);
                final double  y = position.getOrdinate(1);
                final double cx = Math.floor(x / GRID_SQUARE_SIZE);
                final double cy = Math.floor(y / GRID_SQUARE_SIZE);
                int col = (int) cx;
                int row = (int) cy;
                if (isUTM) {
                    /*
                     * Specification said that 100,000-meters columns are lettered from A through Z (omitting I and O)
                     * starting at the 180° meridian, proceeding easterly for 18°, and repeating for each 18° intervals.
                     * Since a UTM zone is 6° width, a 18° interval is exactly 3 standard UTM zones. Columns in zone 1
                     * are A-H, zone 2 are J-R (skipping O), zone 3 are S-Z, then repeating every 3 zones.
                     */
                    if (col < 1 || col > 8) {
                        /*
                         * UTM northing values at the equator range from 166021 to 833979 meters approximately
                         * (WGS84 ellipsoid). Consequently 'cx' ranges from approximately 1.66 to 8.34, so 'c'
                         * should range from 1 to 8 inclusive.
                         */
                        throw new GazetteerException(Errors.format(Errors.Keys.OutsideDomainOfValidity));
                    }
                    switch (zone % 3) {                          // First A-H sequence starts at zone number 1.
                        case 1: col += ('A' - 1); break;
                        case 2: col += ('J' - 1); if (col >= EXCLUDE_O) col++; break;
                        case 0: col += ('S' - 1); break;
                    }
                    /*
                     * Rows in odd  zones are ABCDEFGHJKLMNPQRSTUV
                     * Rows in even zones are FGHJKLMNPQRSTUVABCDE
                     * Those 20 letters are repeated in a cycle.
                     */
                    if ((zone & 1) == 0) {
                        row += ('F' - 'A');
                    }
                    row %= GRID_ROW_COUNT;
                    // Row calculation to be completed after the 'else' block.
                } else {
                    /*
                     * Universal Polar Stereographic (UPS) case. Row letters go from A to Z, omitting I and O.
                     * The column letters go from A to Z, omitting I, O, D, E, M, N, V, W. Rightmost column in
                     * grid zones A and Y has column letter Z, and the next column in grid zones B and Z starts
                     * over with column letter A.
                     */
                    final byte[] columns = POLAR_COLUMNS;
                    col -= PolarStereographicA.UPS_SHIFT / GRID_SQUARE_SIZE;
                    if (!(λ >= 0)) {                    // Same condition than in GZD block. Use of ! is for NaN.
                        col += columns.length;          // Letters Z to A from right to left.
                    }
                    if (col < 0 || col >= columns.length) {
                        throw new GazetteerException(Errors.format(Errors.Keys.OutsideDomainOfValidity));
                    }
                    col  = columns[col];
                    row -= owner.getReferenceSystem().polarOffset(signedZone < 0);
                }
                row += 'A';
                if (row >= EXCLUDE_I && ++row >= EXCLUDE_O) row++;
                buffer.append(separator).append(letter(col)).append(letter(row));
                /*
                 * Numerical location at the given precision.
                 * The specification requires us to truncate the number, not to round it.
                 */
                if (digits > 0) {
                    final double precision = MathFunctions.pow10(METRE_PRECISION_DIGITS - digits);
                    append(buffer.append(separator), (int) ((x - cx * GRID_SQUARE_SIZE) / precision), digits);
                    append(buffer.append(separator), (int) ((y - cy * GRID_SQUARE_SIZE) / precision), digits);
                }
            }
            return buffer.toString();
        }

        /**
         * Returns the given character as a {@code char} if it is a letter, or throws an exception otherwise.
         * The exception should never happen, unless the the encoder is used for a planet larger than Earth
         * for which we do not have enough letters.
         */
        private static char letter(final int c) throws GazetteerException {
            if (c >= 'A' && c <= 'Z') return (char) c;
            throw new GazetteerException(Errors.format(Errors.Keys.OutsideDomainOfValidity));
        }

        /**
         * Appends the given value in the given buffer, padding with zero digits in order to get
         * the specified total amount of digits.
         */
        private static void append(final StringBuilder buffer, final int value, int digits) throws GazetteerException {
            if (value >= 0) {
                final int p = buffer.length();
                digits -= (buffer.append(value).length() - p);
                if (digits >= 0) {
                    StringBuilders.repeat(buffer, p, '0', digits);
                    return;
                }
            }
            throw new GazetteerException(Errors.format(Errors.Keys.OutsideDomainOfValidity));
        }
    }


    /**
     * The result of decoding a MGRS reference.
     * The {@linkplain #getPosition() position} represents the centroid of the decoded MGRS reference.
     *
     * @author  Martin Desruisseaux (Geomatys)
     * @version 0.8
     * @since   0.8
     * @module
     */
    static final class Decoder extends SimpleLocation.Projected {
        /**
         * Number of bits reserved for storing the minimal northing value of latitude bands in the
         * {@link #ROW_RESOLVER} table.
         */
        static final int NORTHING_BITS_COUNT = 4;

        /**
         * Mask over the {@value #NORTHING_BITS_COUNT} lowest bits.
         */
        static final int NORTHING_BITS_MASK = (1 << NORTHING_BITS_COUNT) - 1;

        /**
         * Mapping from latitude bands to the minimal northing values, together with the set of valid rows.
         * There is 20 latitude bands identified by letters C to X inclusively, with letters I and O skipped.
         * First, the band letter must be converted to an index <var>i</var> in the [0 … 19] range.
         * Then, the {@code ROW_RESOLVER[i]} integer values provides the following information:
         *
         * <ul>
         *   <li>
         *     The lowest {@value #NORTHING_BITS_COUNT} bits gives the minimal northing value of all valid
         *     coordinates in the latitude band, as a multiple of the number of metres in a full cycle of
         *     {@value #GRID_ROW_COUNT} rows. That northing value can be computed in metre as below:
         *
         *     {@preformat java
         *         double northing = (ROW_RESOLVER[i] & NORTHING_BITS_MASK) * (GRID_SQUARE_SIZE * GRID_ROW_COUNT);
         *     }
         *   </li><li>
         *     Given a row number <var>row</var> in the [0 … 19] range, the following expression tells
         *     if that row can be inside the latitude band:
         *
         *     {@preformat java
         *         boolean isValidRow = (ROW_RESOLVER[i] & (1 << (row + NORTHING_BITS_COUNT))) != 0;
         *     }
         *
         *     Note that the same row may be valid in two consecutive latitude bands.
         *     The trailing {@code _0000} parts make room for {@value #NORTHING_BITS_COUNT} bits.
         *   </li>
         * </ul>
         *
         * The content of this table is verified by {@code MilitaryGridReferenceSystemTest.verifyDecoderTables()}.
         * The same method with minor edition can be used for generating this table.
         */
        @SuppressWarnings("PointlessBitwiseExpression")
        private static final int[] ROW_RESOLVER = {
            /*
             * First column is the minimal northing value as a multiple of 2000 km.
             * Second column enumerates the valid rows in that latitude band.
             * Trailing _0000 are for making room for NORTHING_BITS_COUNT bits.
             */
            /* Latitude band C (from 80°S) */   0  |  0b11111111100000000001_0000,
            /* Latitude band D (from 72°S) */   1  |  0b00000000001111111111_0000,
            /* Latitude band E (from 64°S) */   1  |  0b00111111111100000000_0000,
            /* Latitude band F (from 56°S) */   1  |  0b11100000000001111111_0000,
            /* Latitude band G (from 48°S) */   2  |  0b00001111111111000000_0000,
            /* Latitude band H (from 40°S) */   2  |  0b11111000000000011111_0000,
            /* Latitude band J (from 32°S) */   3  |  0b00000011111111110000_0000,
            /* Latitude band K (from 24°S) */   3  |  0b11111110000000000111_0000,
            /* Latitude band L (from 16°S) */   4  |  0b00000000111111111100_0000,
            /* Latitude band M (from  8°S) */   4  |  0b11111111100000000000_0000,
            /* Latitude band N (from  0° ) */   0  |  0b00000000000111111111_0000,
            /* Latitude band P (from  8°N) */   0  |  0b00111111111100000000_0000,
            /* Latitude band Q (from 16°N) */   0  |  0b11100000000001111111_0000,
            /* Latitude band R (from 24°N) */   1  |  0b00001111111111000000_0000,
            /* Latitude band S (from 32°N) */   1  |  0b11111000000000011111_0000,
            /* Latitude band T (from 40°N) */   2  |  0b00000011111111110000_0000,
            /* Latitude band U (from 48°N) */   2  |  0b11111110000000000111_0000,
            /* Latitude band V (from 56°N) */   3  |  0b00000000111111111100_0000,
            /* Latitude band W (from 64°N) */   3  |  0b11111111110000000000_0000,
            /* Latitude band X (from 72°N) */   3  |  0b10000011111111111111_0000
        };

        /**
         * The Coordinate Reference System of the decoded MGRS reference.
         * This is an Universal Transverse Mercator (UTM) or Universal Polar Stereographic (UPS) projection.
         *
         * @see #getCoordinateReferenceSystem()
         */
        private final ProjectedCRS crs;

        /**
         * Decodes the given MGRS reference.
         *
         * @param owner  the {@code Coder} which is creating this {@code Decoder}.
         */
        Decoder(final Coder owner, final CharSequence reference) throws TransformException {
            super(owner.getReferenceSystem().rootType(), reference);
            final int zone;                     // UTM zone, or 0 if UPS.
            boolean hasSquareIdentification;    // Whether a square identification is present.
            final double φs;                    // Southernmost bound of latitude band (UTM only).
            final double λ0;                    // Central meridian of UTM zone (ignoring Norway and Svalbard).
            boolean isValid = true;             // Whether the given reference passes consistency checks.

            final int end  = CharSequences.skipTrailingWhitespaces(reference, 0, reference.length());
            final int base = CharSequences.skipLeadingWhitespaces (reference, 0, end);
            int i = endOfDigits(reference, base, end);
            if (i == base && i < end) {
                /*
                 * Universal Polar Stereographic (UPS) case. The reference has no zone number
                 * and begins directly with 3 parts, where each part is exactly one letter:
                 *
                 *   part 0  —  A zone indicator: A or B for South pole, or Y or Z for North pole.
                 *   part 1  —  column letter:    A to Z, omitting I, O, D, E, M, N, V, W.
                 *   part 2  —  row letter:       A to Z, omitting I, O.
                 */
                zone          = 0;                  // Not used in UPS case.
                boolean south = false;              // True for A and B zones, false for Y and Z zones.
                boolean west  = false;              // True for A and Y zones, false for B and Z zones.
                int col = 0, row = 0;               // Parsed row and column indices.
                for (int part = 0; part <= 2; part++) {
                    int c = Character.codePointAt(reference, i);
                    final int ni = i + Character.charCount(c);
                    if (isLetter(c) || isLetter(c -= ('a' - 'A'))) {
parse:                  switch (part) {
                            case 0: {
                                switch (c) {
                                    case 'A': south = true; west = true; break;
                                    case 'B': south = true;              break;
                                    case 'Y':               west = true; break;
                                    case 'Z':                            break;
                                    default : break parse;                              // Invalid UPS zone.
                                }
                                i = nextComponent(owner, reference, base, ni, end);
                                continue;
                            }
                            case 1: {
                                col = Arrays.binarySearch(POLAR_COLUMNS, (byte) c);
                                if (col < 0) break;                                     // Invalid column letter.
                                if (west) col -= POLAR_COLUMNS.length;
                                col += PolarStereographicA.UPS_SHIFT / GRID_SQUARE_SIZE;
                                i = nextComponent(owner, reference, base, ni, end);
                                continue;
                            }
                            case 2: {
                                if (c >= EXCLUDE_O) c--;
                                if (c >= EXCLUDE_I) c--;
                                row = (c - 'A') + owner.getReferenceSystem().polarOffset(south);
                                i = ni;
                                continue;
                            }
                        }
                    }
                    /*
                     * We reach this point only if the current character is invalid (not a letter,
                     * or not one of the letters expected for the current part).
                     */
                    final short key;
                    final CharSequence token;
                    if (part == 0) {
                        key = Resources.Keys.IllegalUPSZone_1;
                        token = reference.subSequence(i, ni);
                    } else {
                        key = Resources.Keys.IllegalSquareIdentification_1;
                        token = CharSequences.token(reference, i);
                    }
                    throw new GazetteerException(Resources.format(key, token));
                }
                crs  = owner.projection(φs = (south ? Latitude.MIN_VALUE : Latitude.MAX_VALUE), 0);
                minX = col * GRID_SQUARE_SIZE;
                minY = row * GRID_SQUARE_SIZE;
                hasSquareIdentification = true;
                λ0 = 0;
            } else {
                /*
                 * Universal Transverse Mercator (UTM) case.
                 */
                zone = parseInt(reference, base, i, Resources.Keys.IllegalUTMZone_1);
                if (zone < 1 || zone > 60) {
                    throw new GazetteerException(Resources.format(Resources.Keys.IllegalUTMZone_1, zone));
                }
                /*
                 * Parse the sub-sequence made of letters. That sub-sequence can have one or three parts.
                 * The first part is mandatory and the two other parts are optional, but if the two last
                 * parts are omitted, then they must be omitted together.
                 *
                 *   part 0  —  latitude band: C-X (excluding I and O) for UTM. Other letters (A, B, Y, Z) are for UPS.
                 *   part 1  —  column letter: A-H in zone 1, J-R (skipping O) in zone 2, S-Z in zone 3, then repeat.
                 *   part 2  —  row letter:    ABCDEFGHJKLMNPQRSTUV in odd zones, FGHJKLMNPQRSTUVABCDE in even zones.
                 */
                int latitudeBand = -1;                              // The latitude band in [0 … 19] range.
                int col = 1, row = 0;
                hasSquareIdentification = true;
                for (int part = 0; part <= 2; part++) {
                    if (part == 1 && i >= end) {
                        hasSquareIdentification = false;
                        break;                                      // Allow to stop parsing only after part 1.
                    }
                    i = nextComponent(owner, reference, base, i, end);
                    int c = Character.codePointAt(reference, i);
                    final int ni = i + Character.charCount(c);
                    if (!isLetter(c) && !isLetter(c -= ('a' - 'A'))) {
                        final short key;
                        final CharSequence token;
                        if (part == 0) {
                            key = Resources.Keys.IllegalLatitudeBand_1;
                            token = reference.subSequence(i, ni);
                        } else {
                            key = Resources.Keys.IllegalSquareIdentification_1;
                            token = CharSequences.token(reference, i);
                        }
                        throw new GazetteerException(Resources.format(key, token));
                    }
                    /*
                     * At this point, 'c' is a valid letter. First, applies a correction for the fact that 'I' and 'O'
                     * letters were excluded. Next, the conversion to latitude or 100 000 meters grid indices depends
                     * on which part we are parsing. The formulas used below are about the same than in Encoder class,
                     * with terms moved on the other side of the equations.
                     */
                    if (c >= EXCLUDE_O) c--;
                    if (c >= EXCLUDE_I) c--;
                    switch (part) {
                        case 0: {
                            latitudeBand = (c - 'C');
                            break;
                        }
                        case 1: {
                            switch (zone % 3) {                         // First A-H sequence starts at zone number 1.
                                case 1: col = c - ('A' - 1); break;
                                case 2: col = c - ('J' - 2); break;     // -2 because 'I' has already been excluded.
                                case 0: col = c - ('S' - 3); break;     // -3 because 'I' and 'O' have been excluded.
                            }
                            break;
                        }
                        case 2: {
                            if ((zone & 1) != 0) {
                                row = c - 'A';
                            } else {
                                row = c - 'F';
                                if (row < 0) {
                                    row += GRID_ROW_COUNT;
                                }
                            }
                            break;
                        }
                    }
                    i = ni;
                }
                /*
                 * Create a UTM projection for exactly the zone specified in the MGRS reference, regardless the
                 * Norway and Svalbard special cases.  Then fetch the smallest northing value that the latitude
                 * band may contain in that UTM zone. This is the projection of the geographic coordinate on the
                 * central meridian in the North hemisphere, or on the UTM zone border in the South hemisphere
                 * (because the zone center is always closer to equator and the zone borders closer to poles).
                 * This estimation is needed because the 100 kilometres square identification is insufficient;
                 * we may need to add some multiple of 2000 kilometres (20 squares).
                 */
                φs = latitudeBand * LATITUDE_BAND_HEIGHT + TransverseMercator.Zoner.SOUTH_BOUNDS;
                if (latitudeBand < 0 || latitudeBand >= ROW_RESOLVER.length) {
                    throw new GazetteerException(Resources.format(Resources.Keys.IllegalLatitudeBand_1, Encoder.latitudeBand(φs)));
                }
                λ0  = ZONER.centralMeridian(zone);
                crs = owner.projection(Math.signum(φs), λ0);
                final int info = ROW_RESOLVER[latitudeBand];        // Contains the above-cited northing value.
                if (hasSquareIdentification) {
                    int rowBit = 1 << (row + NORTHING_BITS_COUNT);  // Bit mask of the row to check for existence.
                    isValid = (info & rowBit) != 0;                 // Whether the row exists in the latitude band.
                    if (isValid) {
                        /*
                         * Get the bit mask of the first invalid row north of current row. For example if 'info' has
                         * the following value and the current 'rowBit' mask has its bit set at the position of (1),
                         * then after this line of code the 'rowBit' bit will be set at the position of (2).
                         *
                         *     11111000000000011111
                         *         ^         ^   ^
                         *        (3)       (2) (1)
                         *
                         * This line works by forcing to 1 all bits on the right side of (1) then searching the lowest
                         * zero bit (actually inverting everything and searching the lowest one bit, since there is no
                         * "lowestZeroBit" method in java.lang.Integer).
                         *
                         * The purpose is to verify if there another, distinct, sequence of 1 values on the left of
                         * current sequence. We will check that (after this block) by verifying if there is any bit
                         * set to 1 on the left side of zero bits, like the bit at position (3) in above figure. If
                         * such distinct sequence exists, then that sequence is the current cycle of rows while the
                         * sequence on the right side is a new cycle of rows.
                         *
                         * Note that if there is no zero bit on the left side, then rowBit = 0.
                         * This implies ~(rowBit - 1) == 0, which is okay for next line of code.
                         * Note: ~(rowBit - 1)  ==  -rowBit
                         */
                        rowBit = Integer.lowestOneBit(~(info | (rowBit - 1)));
                    }
                    if ((info & -rowBit) != 0) {    // Test if there is valid rows on the left side of sequence of zero bits.
                        row += GRID_ROW_COUNT;      // Left bits were from previous cycle, which means that we started a new cycle.
                    }
                }
                row += (info & NORTHING_BITS_MASK) * GRID_ROW_COUNT;        // Add the pre-computed northing value.
                minX = col * GRID_SQUARE_SIZE;
                minY = row * GRID_SQUARE_SIZE;
            }
            /*
             * If we have not yet reached the end of string, parse the numerical location.
             * That location is normally encoded as a single number with an even number of digits.
             * The first half is the easting and the second half is the northing, both relative to the
             * 100 kilometer square. However some variants of MGRS use a separator, in which case we get
             * two distinct numbers. In both cases, the resolution is determined by the amount of digits.
             */
            final double sx, sy;    // Scale factors for converting MGRS values in to easting and northing in metres.
            if (i < end) {
                i = nextComponent(owner, reference, base, i, end);
                int s = endOfDigits(reference, i, end);
                final double x, y;
                if (s >= end) {
                    int length = s - i;
                    if ((length & 1) != 0) {
                        throw new GazetteerException(Resources.format(Resources.Keys.OddGridCoordinateLength_1,
                                reference.subSequence(i, s)));
                    }
                    final int h = i + (length >>>= 1);
                    sx = sy = MathFunctions.pow10(METRE_PRECISION_DIGITS - length);
                    x  = parseCoordinate(reference, i, h, sx);
                    y  = parseCoordinate(reference, h, s, sy);
                } else {
                    sx = MathFunctions.pow10(METRE_PRECISION_DIGITS - (s - i));
                    x  = parseCoordinate(reference, i, s, sx);
                    i  = nextComponent(owner, reference, base, s, end);
                    s  = endOfDigits(reference, i, end);
                    sy = MathFunctions.pow10(METRE_PRECISION_DIGITS - (s - i));
                    y = parseCoordinate(reference, i, s, sy);
                    if (s < end) {
                        throw new GazetteerException(Errors.format(Errors.Keys.UnexpectedCharactersAfter_2,
                                reference.subSequence(base, s), CharSequences.trimWhitespaces(reference, s, end)));
                    }
                }
                minX += x;
                minY += y;
            } else if (hasSquareIdentification) {
                sx = sy = GRID_SQUARE_SIZE;
            } else {
                /*
                 * Not used for scaling anymore. Choose value that will cause the point to be in zone zenter.
                 * The '- GRID_SQUARE_SIZE' in 'sx' is because the westernmost 100-km grid squares are in the
                 * column at index 1 (the column at index 0 is outside the UTM zone).
                 */
                sx = (ZONER.easting - GRID_SQUARE_SIZE) * 2;
                sy =  ZONER.northing;
            }
            maxX = minX + sx;
            maxY = minY + sy;
            /*
             * At this point the non-clipped projected envelope has been computed. Now compute the geographic envelope.
             * We need this information for clipping the projected envelope to the domain of validity of UTM zone.
             */
            if (!hasSquareIdentification) {
                if (zone != 0) {
                    if (φs < 0) {
                        southBoundLatitude = TransverseMercator.Zoner.SOUTH_BOUNDS;
                        northBoundLatitude = 0;
                    } else {
                        southBoundLatitude = 0;
                        northBoundLatitude = TransverseMercator.Zoner.NORTH_BOUNDS;
                    }
                    westBoundLongitude = λ0 - ZONER.width / 2;
                    eastBoundLongitude = λ0 + ZONER.width / 2;
                } else {
                    if (φs < 0) {
                        southBoundLatitude = Latitude.MIN_VALUE;
                        northBoundLatitude = TransverseMercator.Zoner.SOUTH_BOUNDS;
                    } else {
                        southBoundLatitude = TransverseMercator.Zoner.NORTH_BOUNDS;
                        northBoundLatitude = Latitude.MAX_VALUE;
                    }
                    westBoundLongitude = Longitude.MIN_VALUE;
                    eastBoundLongitude = Longitude.MAX_VALUE;
                }
            } else {
                final MathTransform projection = crs.getConversionFromBase().getMathTransform();
                computeGeographicBoundingBox(projection.inverse());
                /*
                 * Update the LocationType according the precision.
                 * This update is mostly for documentation purpose.
                 * There is three levels:
                 *
                 *   - Grid zone designator         (if hasSquareIdentification == false)
                 *   - 100 km square identifier     (if resolution == 100 km)
                 *   - Grid coordinate              (if resolution < 100 km)
                 */
                setTypeToChild();                   // Replace "Grid zone designator" by "100 km square identifier".
                if (sx < GRID_SQUARE_SIZE || sy < GRID_SQUARE_SIZE) {
                    setTypeToChild();               // Replace "100 km square identifier" by "Grid coordinate".
                }
                /*
                 * At this point we finished computing the position. Now perform error detection, by verifying
                 * if the given 100 kilometres square identification is consistent with grid zone designation.
                 * We verify both φ and λ, but the verification of φ is actually redundant with the check of
                 * 100 km square validity that we did previously with the help of ROW_RESOLVER bitmask.
                 * We check φ anyway in case of bug.
                 */
                if (isValid && zone != 0) {
                    isValid = (northBoundLatitude >= φs) && (southBoundLatitude < upperBound(φs));
                    if (isValid) {
                        /*
                         * Verification of UTM zone. We allow a tolerance for latitudes close to a pole because
                         * not all users may apply the UTM special rules for Norway and Svalbard. Anyway, using
                         * the neighbor zone at those high latitudes is less significant. For other latitudes,
                         * we allow a tolerance if the point is close to a line of zone change.
                         */
                        final double λ = (westBoundLongitude + eastBoundLongitude) / 2;
                        final double φ = (southBoundLatitude + northBoundLatitude) / 2;
                        int zoneError = ZONER.zone(φ, λ) - zone;
                        if (zoneError != 0) {
                            final int zc = ZONER.zoneCount();
                            if (zoneError < zc/-2) zoneError += zc;         // If change of zone crosses the anti-meridian.
                            if (zoneError > zc/+2) zoneError -= zc;
                            if (ZONER.isSpecialCase(zone, φ)) {
                                isValid = Math.abs(zoneError) <= 2;         // Tolerance in zone numbers for high latitudes.
                            } else {
                                final double rλ = Math.IEEEremainder(λ - ZONER.origin, ZONER.width);    // Distance to closest zone change, in degrees of longitude.
                                final double cv = (minX - ZONER.easting) / (λ - λ0);                    // Approximate conversion factor from degrees to metres.
                                isValid = (Math.abs(rλ) * cv <= sx);                                    // Be tolerant if distance in metres is less than resolution.
                                if (isValid) {
                                    isValid = (zoneError == (rλ < 0 ? -1 : +1));                        // Verify also that the error is on the side of the zone change.
                                }
                            }
                        }
                    }
                }
                /*
                 * At this point we finished verifying the cell validity using the coordinates specified by the
                 * MGRS reference. If the cell is valid, we can now check for cells that are on a zone border.
                 * Those cells will be clipped to the zone valid area.
                 */
                if (isValid && owner.getClipToValidArea()) {
                    final boolean changed;
                    if (zone != 0) {
                        double width = ZONER.width;
                        if (!ZONER.isSpecialCase(zone, φs)) width /= 2;       // Be strict only if not Norway or Svalbard.
                        changed = clipGeographicBoundingBox(λ0 - width, φs,
                                                            λ0 + width, upperBound(φs));
                    } else if (φs < 0) {
                        changed = clipGeographicBoundingBox(Longitude.MIN_VALUE, Latitude.MIN_VALUE,
                                                            Longitude.MAX_VALUE, TransverseMercator.Zoner.SOUTH_BOUNDS);
                    } else {
                        changed = clipGeographicBoundingBox(Longitude.MIN_VALUE, TransverseMercator.Zoner.NORTH_BOUNDS,
                                                            Longitude.MAX_VALUE, Latitude.MAX_VALUE);
                    }
                    if (changed) {
                        clipProjectedEnvelope(projection, sx / 100, sy / 100);
                    }
                }
            }
            if (!isValid) {
                final String gzd;
                try {
                    gzd = owner.encoder(crs).encode(owner, getDirectPosition(), true, "", 0);
                } catch (IllegalArgumentException | FactoryException e) {
                    throw new GazetteerException(e.getLocalizedMessage(), e);
                }
                final CharSequence ref = reference.subSequence(base, end);
                throw new ReferenceVerifyException(Resources.format(Resources.Keys.InconsistentWithGZD_2, ref, gzd));
            }
        }

        /**
         * Skips spaces, then the separator if present (optional).
         *
         * @param  reference  the reference to parse.
         * @param  base       index where the parsing began. Used for formatting error message only.
         * @param  start      current parsing position.
         * @param  end        where the parsing is expected to end.
         * @return position where to continue parsing (with spaces skipped).
         * @throws GazetteerException if this method unexpectedly reached the end of string.
         */
        private static int nextComponent(final Coder owner, final CharSequence reference,
                final int base, int start, final int end) throws GazetteerException
        {
            start = CharSequences.skipLeadingWhitespaces(reference, start, end);
            if (start < end) {
                if (!CharSequences.regionMatches(reference, start, owner.trimmedSeparator)) {
                    return start;               // Separator not found, but it was optional.
                }
                start += owner.trimmedSeparator.length();
                start = CharSequences.skipLeadingWhitespaces(reference, start, end);
                if (start < end) {
                    return start;
                }
            }
            throw new GazetteerException(Errors.format(Errors.Keys.UnexpectedEndOfString_1, reference.subSequence(base, end)));
        }

        /**
         * Returns {@code true} if the given character is a valid upper case ASCII letter, excluding I and O.
         */
        private static boolean isLetter(final int c) {
            return (c >= 'A' && c <= 'Z') && c != EXCLUDE_I && c != EXCLUDE_O;
        }

        /**
         * Returns the index after the last digit in a sequence of ASCII characters.
         * Leading whitespaces must have been skipped before to invoke this method.
         */
        private static int endOfDigits(final CharSequence reference, int i, final int end) {
            while (i < end) {
                final char c = reference.charAt(i);     // Code-point API not needed here because we restrict to ASCII.
                if (c < '0' || c > '9') {               // Do not use Character.isDigit(…) because we restrict to ASCII.
                    break;
                }
                i++;
            }
            return i;
        }

        /**
         * Parses part of the given character sequence as an integer.
         *
         * @param  reference  the MGRS reference to parse.
         * @param  start      index of the first character to parse as an integer.
         * @param  end        index after the last character to parse as an integer.
         * @param  errorKey   {@link Resources.Keys} value to use in case of error.
         *                    The error message string shall accept exactly one argument.
         * @return the parsed integer.
         * @throws GazetteerException if the string can not be parsed as an integer.
         */
        private static int parseInt(final CharSequence reference, final int start, final int end, final short errorKey)
                throws GazetteerException
        {
            NumberFormatException cause = null;
            final CharSequence part;
            if (start == end) {
                part = CharSequences.token(reference, start);
            } else {
                part = reference.subSequence(start, end);
                try {
                    return Integer.parseInt(part.toString());
                } catch (NumberFormatException e) {
                    cause = e;
                }
            }
            throw new GazetteerException(Resources.format(errorKey, part), cause);
        }

        /**
         * Parses part of the given character sequence as a grid coordinate.
         * The resolution is determined by the amount of digits.
         *
         * @param  reference  the MGRS reference to parse.
         * @param  start      index of the first character to parse as a grid coordinate.
         * @param  end        index after the last character to parse as a grid coordinate.
         * @param  scale      value of {@code MathFunctions.pow10(METRE_PRECISION_DIGITS - (end - start))}.
         * @return the parsed grid coordinate (also referred to as rectangular coordinates).
         * @throws GazetteerException if the string can not be parsed as a grid coordinate.
         */
        private static double parseCoordinate(final CharSequence reference,
                final int start, final int end, final double scale) throws GazetteerException
        {
            return parseInt(reference, start, end, Resources.Keys.IllegalGridCoordinate_1) * scale;
        }

        /**
         * Returns the upper bound of the latitude band specified by the given lower bound.
         */
        static double upperBound(final double φ) {
            return φ < TransverseMercator.Zoner.SVALBARD_BOUNDS ? φ + LATITUDE_BAND_HEIGHT
                     : TransverseMercator.Zoner.NORTH_BOUNDS;
        }

        /**
         * Returns the Coordinate Reference System of the decoded MGRS reference.
         * This is an Universal Transverse Mercator (UTM) or Universal Polar Stereographic (UPS) projection.
         */
        @Override
        public CoordinateReferenceSystem getCoordinateReferenceSystem() {
            return crs;
        }
    }
}