/*
 * 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;

import java.util.Map;
import java.util.List;
import java.util.ArrayList;
import java.util.Collections;
import java.util.logging.Filter;
import java.util.logging.LogRecord;
import org.opengis.util.FactoryException;
import org.opengis.geometry.Envelope;
import org.opengis.referencing.NoSuchAuthorityCodeException;
import org.opengis.referencing.IdentifiedObject;
import org.opengis.referencing.cs.CartesianCS;
import org.opengis.referencing.cs.EllipsoidalCS;
import org.opengis.referencing.cs.AxisDirection;
import org.opengis.referencing.cs.CoordinateSystem;
import org.opengis.referencing.cs.CoordinateSystemAxis;
import org.opengis.referencing.crs.CRSFactory;
import org.opengis.referencing.crs.SingleCRS;
import org.opengis.referencing.crs.CompoundCRS;
import org.opengis.referencing.crs.CoordinateReferenceSystem;
import org.opengis.referencing.crs.CRSAuthorityFactory;
import org.opengis.referencing.crs.GeodeticCRS;
import org.opengis.referencing.crs.GeographicCRS;
import org.opengis.referencing.crs.GeneralDerivedCRS;
import org.opengis.referencing.crs.ProjectedCRS;
import org.opengis.referencing.crs.TemporalCRS;
import org.opengis.referencing.crs.VerticalCRS;
import org.opengis.referencing.crs.EngineeringCRS;
import org.opengis.referencing.datum.Datum;
import org.opengis.referencing.datum.GeodeticDatum;
import org.opengis.referencing.operation.Conversion;
import org.opengis.referencing.operation.OperationNotFoundException;
import org.opengis.metadata.citation.Citation;
import org.opengis.metadata.extent.Extent;
import org.opengis.metadata.extent.GeographicBoundingBox;
import org.opengis.referencing.operation.CoordinateOperation;
import org.opengis.referencing.operation.TransformException;
import org.apache.sis.measure.Units;
import org.apache.sis.geometry.Envelopes;
import org.apache.sis.geometry.GeneralEnvelope;
import org.apache.sis.internal.referencing.AxisDirections;
import org.apache.sis.internal.referencing.EllipsoidalHeightCombiner;
import org.apache.sis.internal.referencing.PositionalAccuracyConstant;
import org.apache.sis.internal.referencing.CoordinateOperations;
import org.apache.sis.internal.referencing.ReferencingUtilities;
import org.apache.sis.internal.referencing.DefinitionVerifier;
import org.apache.sis.internal.referencing.Resources;
import org.apache.sis.internal.system.DefaultFactories;
import org.apache.sis.internal.system.Modules;
import org.apache.sis.internal.system.Loggers;
import org.apache.sis.internal.util.Numerics;
import org.apache.sis.referencing.cs.AxisFilter;
import org.apache.sis.referencing.cs.CoordinateSystems;
import org.apache.sis.referencing.cs.DefaultVerticalCS;
import org.apache.sis.referencing.crs.DefaultGeographicCRS;
import org.apache.sis.referencing.crs.DefaultProjectedCRS;
import org.apache.sis.referencing.crs.DefaultVerticalCRS;
import org.apache.sis.referencing.crs.DefaultCompoundCRS;
import org.apache.sis.referencing.crs.DefaultEngineeringCRS;
import org.apache.sis.referencing.operation.AbstractCoordinateOperation;
import org.apache.sis.referencing.operation.CoordinateOperationContext;
import org.apache.sis.referencing.operation.DefaultCoordinateOperationFactory;
import org.apache.sis.referencing.operation.DefaultConversion;
import org.apache.sis.referencing.factory.GeodeticObjectFactory;
import org.apache.sis.referencing.factory.UnavailableFactoryException;
import org.apache.sis.metadata.iso.extent.DefaultGeographicBoundingBox;
import org.apache.sis.metadata.iso.extent.Extents;
import org.apache.sis.util.resources.Errors;
import org.apache.sis.util.logging.Logging;
import org.apache.sis.util.ArgumentChecks;
import org.apache.sis.util.Utilities;
import org.apache.sis.util.Static;


/**
 * Static methods working on {@linkplain CoordinateReferenceSystem Coordinate Reference Systems}.
 * The methods defined in this class can be grouped in three categories:
 *
 * <ul>
 *   <li>Factory methods, the most notable one being {@link #forCode(String)}.</li>
 *   <li>Methods providing information, like {@link #isHorizontalCRS(CoordinateReferenceSystem)}.</li>
 *   <li>Finding coordinate operations between a source and a target CRS.</li>
 * </ul>
 *
 * <div class="section">Usage example</div>
 * The most frequently used methods in this class are {@link #forCode forCode(…)}, {@link #fromWKT fromWKT(…)}
 * and {@link #findOperation findOperation(…)}. An usage example is like below
 * (see the <a href="http://sis.apache.org/tables/CoordinateReferenceSystems.html">Apache SIS™ Coordinate
 * Reference System (CRS) codes</a> page for the complete list of EPSG codes):
 *
 * {@preformat java
 *   CoordinateReferenceSystem source = CRS.forCode("EPSG:4326");                   // WGS 84
 *   CoordinateReferenceSystem target = CRS.forCode("EPSG:3395");                   // WGS 84 / World Mercator
 *   CoordinateOperation operation = CRS.findOperation(source, target, null);
 *   if (CRS.getLinearAccuracy(operation) > 100) {
 *       // If the accuracy is coarser than 100 metres (or any other threshold at application choice)
 *       // maybe the operation is not suitable. Decide here what to do (throw an exception, etc).
 *   }
 *   MathTransform mt = operation.getMathTransform();
 *   DirectPosition position = new DirectPosition2D(20, 30);            // 20°N 30°E   (watch out axis order!)
 *   position = mt.transform(position, position);
 *   System.out.println(position);
 * }
 *
 * <div class="section">Note on kinds of CRS</div>
 * The {@link #getSingleComponents(CoordinateReferenceSystem)} method decomposes an arbitrary CRS into a flat
 * list of single components. In such flat list, vertical and temporal components can easily be identified by
 * {@code instanceof} checks. But identifying the horizontal component is not as easy. The list below suggests
 * ways to classify the components:
 *
 * <ul>
 *   <li><code>if (crs instanceof TemporalCRS)</code> determines if the CRS is for the temporal component.</li>
 *   <li><code>if (crs instanceof VerticalCRS)</code> determines if the CRS is for the vertical component.</li>
 *   <li><code>if (CRS.{@linkplain #isHorizontalCRS(CoordinateReferenceSystem) isHorizontalCRS}(crs))</code>
 *       determines if the CRS is for the horizontal component.</li>
 * </ul>
 *
 * @author  Martin Desruisseaux (IRD, Geomatys)
 * @author  Alexis Manin (Geomatys)
 * @version 1.1
 * @since   0.3
 * @module
 */
public final class CRS extends Static {
    /**
     * Do not allow instantiation of this class.
     */
    private CRS() {
    }

    /**
     * Returns the Coordinate Reference System for the given authority code.
     * The set of available codes depends on the {@link CRSAuthorityFactory} instances available on the classpath.
     * There is many thousands of <a href="http://sis.apache.org/tables/CoordinateReferenceSystems.html">CRS
     * defined by EPSG authority or by other authorities</a>.
     * The following table lists a very small subset of codes which are guaranteed to be available
     * on any installation of Apache SIS:
     *
     * <blockquote><table class="sis">
     *   <caption>Minimal set of supported authority codes</caption>
     *   <tr><th>Code</th>      <th>Enum</th>                            <th>CRS Type</th>        <th>Description</th></tr>
     *   <tr><td>CRS:27</td>    <td>{@link CommonCRS#NAD27  NAD27}</td>  <td>Geographic</td>      <td>Like EPSG:4267 except for (<var>longitude</var>, <var>latitude</var>) axis order</td></tr>
     *   <tr><td>CRS:83</td>    <td>{@link CommonCRS#NAD83  NAD83}</td>  <td>Geographic</td>      <td>Like EPSG:4269 except for (<var>longitude</var>, <var>latitude</var>) axis order</td></tr>
     *   <tr><td>CRS:84</td>    <td>{@link CommonCRS#WGS84  WGS84}</td>  <td>Geographic</td>      <td>Like EPSG:4326 except for (<var>longitude</var>, <var>latitude</var>) axis order</td></tr>
     *   <tr><td>EPSG:4047</td> <td>{@link CommonCRS#SPHERE SPHERE}</td> <td>Geographic</td>      <td>GRS 1980 Authalic Sphere</td></tr>
     *   <tr><td>EPSG:4230</td> <td>{@link CommonCRS#ED50   ED50}</td>   <td>Geographic</td>      <td>European Datum 1950</td></tr>
     *   <tr><td>EPSG:4258</td> <td>{@link CommonCRS#ETRS89 ETRS89}</td> <td>Geographic</td>      <td>European Terrestrial Reference Frame 1989</td></tr>
     *   <tr><td>EPSG:4267</td> <td>{@link CommonCRS#NAD27  NAD27}</td>  <td>Geographic</td>      <td>North American Datum 1927</td></tr>
     *   <tr><td>EPSG:4269</td> <td>{@link CommonCRS#NAD83  NAD83}</td>  <td>Geographic</td>      <td>North American Datum 1983</td></tr>
     *   <tr><td>EPSG:4322</td> <td>{@link CommonCRS#WGS72  WGS72}</td>  <td>Geographic</td>      <td>World Geodetic System 1972</td></tr>
     *   <tr><td>EPSG:4326</td> <td>{@link CommonCRS#WGS84  WGS84}</td>  <td>Geographic</td>      <td>World Geodetic System 1984</td></tr>
     *   <tr><td>EPSG:4936</td> <td>{@link CommonCRS#ETRS89 ETRS89}</td> <td>Geocentric</td>      <td>European Terrestrial Reference Frame 1989</td></tr>
     *   <tr><td>EPSG:4937</td> <td>{@link CommonCRS#ETRS89 ETRS89}</td> <td>Geographic 3D</td>   <td>European Terrestrial Reference Frame 1989</td></tr>
     *   <tr><td>EPSG:4978</td> <td>{@link CommonCRS#WGS84  WGS84}</td>  <td>Geocentric</td>      <td>World Geodetic System 1984</td></tr>
     *   <tr><td>EPSG:4979</td> <td>{@link CommonCRS#WGS84  WGS84}</td>  <td>Geographic 3D</td>   <td>World Geodetic System 1984</td></tr>
     *   <tr><td>EPSG:4984</td> <td>{@link CommonCRS#WGS72  WGS72}</td>  <td>Geocentric</td>      <td>World Geodetic System 1972</td></tr>
     *   <tr><td>EPSG:4985</td> <td>{@link CommonCRS#WGS72  WGS72}</td>  <td>Geographic 3D</td>   <td>World Geodetic System 1972</td></tr>
     *   <tr><td>EPSG:5041</td> <td>{@link CommonCRS#WGS84  WGS84}</td>  <td>Projected</td>       <td>WGS 84 / UPS North (E,N)</td></tr>
     *   <tr><td>EPSG:5042</td> <td>{@link CommonCRS#WGS84  WGS84}</td>  <td>Projected</td>       <td>WGS 84 / UPS South (E,N)</td></tr>
     *   <tr><td>EPSG:322##</td><td>{@link CommonCRS#WGS72  WGS72}</td>  <td>Projected</td>       <td>WGS 72 / UTM zone ##N</td></tr>
     *   <tr><td>EPSG:323##</td><td>{@link CommonCRS#WGS72  WGS72}</td>  <td>Projected</td>       <td>WGS 72 / UTM zone ##S</td></tr>
     *   <tr><td>EPSG:326##</td><td>{@link CommonCRS#WGS84  WGS84}</td>  <td>Projected</td>       <td>WGS 84 / UTM zone ##N</td></tr>
     *   <tr><td>EPSG:327##</td><td>{@link CommonCRS#WGS84  WGS84}</td>  <td>Projected</td>       <td>WGS 84 / UTM zone ##S</td></tr>
     *   <tr><td>EPSG:5715</td> <td>{@link CommonCRS.Vertical#DEPTH DEPTH}</td> <td>Vertical</td> <td>Mean Sea Level depth</td></tr>
     *   <tr><td>EPSG:5714</td> <td>{@link CommonCRS.Vertical#MEAN_SEA_LEVEL MEAN_SEA_LEVEL}</td> <td>Vertical</td> <td>Mean Sea Level height</td></tr>
     * </table></blockquote>
     *
     * This method accepts also the URN and URL syntaxes.
     * For example the following codes are considered equivalent to {@code "EPSG:4326"}:
     * <ul>
     *   <li>{@code "EPSG::4326"}</li>
     *   <li>{@code "urn:ogc:def:crs:EPSG::4326"}</li>
     *   <li>{@code "http://www.opengis.net/def/crs/epsg/0/4326"}</li>
     *   <li>{@code "http://www.opengis.net/gml/srs/epsg.xml#4326"}</li>
     * </ul>
     *
     * URIs can be combined for creating larger objects. For example the following URIs combine a
     * two-dimensional WGS84 reference system (EPSG:4326) with a Mean Sea Level height (EPSG:5714).
     * The result is a three-dimensional {@linkplain org.apache.sis.referencing.crs.DefaultCompoundCRS
     * compound coordinate reference system}:
     *
     * <ul>
     *   <li>{@code "urn:ogc:def:crs,crs:EPSG::4326,crs:EPSG::5714"}</li>
     *   <li><code>"http://www.opengis.net/def/crs-compound?<br>
     *            1=http://www.opengis.net/def/crs/epsg/0/4326&amp;<br>
     *            2=http://www.opengis.net/def/crs/epsg/0/5714"</code></li>
     * </ul>
     *
     * <p>URNs (but not URLs) can also combine a
     * {@linkplain org.apache.sis.referencing.datum.DefaultGeodeticDatum geodetic datum} with an
     * {@linkplain org.apache.sis.referencing.cs.DefaultEllipsoidalCS ellipsoidal coordinate system} for creating a new
     * {@linkplain org.apache.sis.referencing.crs.DefaultGeographicCRS geographic CRS}, or a base geographic CRS with a
     * {@linkplain org.apache.sis.referencing.operation.DefaultConversion conversion} and a
     * {@linkplain org.apache.sis.referencing.cs.DefaultCartesianCS Cartesian coordinate system} for creating a new
     * {@linkplain org.apache.sis.referencing.crs.DefaultProjectedCRS projected coordinate reference system}.</p>
     *
     * Note that the {@link IdentifiedObjects#lookupURN(IdentifiedObject, Citation)}
     * method can be seen as a converse of this method.
     * More codes may also be supported depending on which extension modules are available.
     * See for example the {@linkplain org.apache.sis.storage.gdal bindings to Proj.4 library}.
     *
     * @param  code  the authority code.
     * @return the Coordinate Reference System for the given authority code.
     * @throws NoSuchAuthorityCodeException if there is no known CRS associated to the given code.
     * @throws FactoryException if the CRS creation failed for an other reason.
     *
     * @see #getAuthorityFactory(String)
     * @see org.apache.sis.referencing.factory.GeodeticAuthorityFactory
     * @see <a href="http://epsg-registry.org/">EPSG Geodetic Registry</a>
     *
     * @category factory
     */
    public static CoordinateReferenceSystem forCode(final String code)
            throws NoSuchAuthorityCodeException, FactoryException
    {
        ArgumentChecks.ensureNonNull("code", code);
        try {
            return AuthorityFactories.ALL.createCoordinateReferenceSystem(code);
        } catch (UnavailableFactoryException e) {
            return AuthorityFactories.fallback(e).createCoordinateReferenceSystem(code);
        }
    }

    /**
     * Creates a Coordinate Reference System object from a <cite>Well Known Text</cite> (WKT).
     * The default {@linkplain org.apache.sis.io.wkt Apache SIS parser} understands both
     * version 1 (a.k.a. OGC 01-009) and version 2 (a.k.a. ISO 19162) of the WKT format.
     *
     * <div class="note"><b>Example:</b> below is a slightly simplified WKT 2 string for a Mercator projection.
     * For making this example smaller, some optional {@code UNIT[…]} and {@code ORDER[…]} elements have been omitted.
     *
     * {@preformat wkt
     *   ProjectedCRS["SIRGAS 2000 / Brazil Mercator",
     *     BaseGeodCRS["SIRGAS 2000",
     *       Datum["Sistema de Referencia Geocentrico para las Americas 2000",
     *         Ellipsoid["GRS 1980", 6378137, 298.257222101]]],
     *     Conversion["Petrobras Mercator",
     *       Method["Mercator (variant B)", Id["EPSG",9805]],
     *       Parameter["Latitude of 1st standard parallel", -2],
     *       Parameter["Longitude of natural origin", -43],
     *       Parameter["False easting", 5000000],
     *       Parameter["False northing", 10000000]],
     *     CS[cartesian,2],
     *       Axis["easting (E)", east],
     *       Axis["northing (N)", north],
     *       LengthUnit["metre", 1],
     *     Id["EPSG",5641]]
     * }
     * </div>
     *
     * If the parsing produced warnings, they will be reported in a logger named {@code "org.apache.sis.io.wkt"}.
     * In particular, this method verifies if the description provided by the WKT matches the description provided
     * by the authority ({@code "EPSG:5641"} in above example) and reports discrepancies.
     * Note that this comparison between parsed CRS and authoritative CRS is specific to this convenience method;
     * other APIs documented in <cite>see also</cite> section do not perform this comparison automatically.
     * Should the WKT description and the authoritative description be in conflict, the WKT description prevails
     * as mandated by ISO 19162 standard (see {@link #fromAuthority fromAuthority(…)} if a different behavior is needed).
     *
     * <div class="section">Usage and performance considerations</div>
     * This convenience method delegates to
     * {@link org.apache.sis.referencing.factory.GeodeticObjectFactory#createFromWKT(String)}
     * using a default factory instance. This is okay for occasional use, but has the following limitations:
     *
     * <ul>
     *   <li>Performance may be sub-optimal in a multi-thread environment.</li>
     *   <li>No control on the WKT {@linkplain org.apache.sis.io.wkt.Convention conventions} in use.</li>
     *   <li>No control on the handling of {@linkplain org.apache.sis.io.wkt.Warnings warnings}.</li>
     * </ul>
     *
     * Applications which need to parse a large amount of WKT strings should consider to use
     * the {@link org.apache.sis.io.wkt.WKTFormat} class instead than this method.
     *
     * @param  text  coordinate system encoded in Well-Known Text format (version 1 or 2).
     * @return the parsed Coordinate Reference System.
     * @throws FactoryException if the given WKT can not be parsed.
     *
     * @see org.apache.sis.io.wkt.WKTFormat
     * @see org.apache.sis.referencing.factory.GeodeticObjectFactory#createFromWKT(String)
     * @see org.apache.sis.geometry.Envelopes#fromWKT(CharSequence)
     * @see <a href="http://docs.opengeospatial.org/is/12-063r5/12-063r5.html">WKT 2 specification</a>
     *
     * @since 0.6
     */
    public static CoordinateReferenceSystem fromWKT(final String text) throws FactoryException {
        ArgumentChecks.ensureNonNull("text", text);
        final CoordinateReferenceSystem crs = DefaultFactories.forBuildin(CRSFactory.class).createFromWKT(text);
        DefinitionVerifier.withAuthority(crs, Loggers.WKT, CRS.class, "fromWKT");
        return crs;
    }

    /**
     * Creates a coordinate reference system object from a XML string.
     * Note that the given argument is the XML document itself, <strong>not</strong> a URL to a XML document.
     * For reading XML documents from readers or input streams,
     * see static methods in the {@link org.apache.sis.xml.XML} class.
     *
     * <p>If the unmarshalling produced warnings, they will be reported in a logger named {@code "org.apache.sis.xml"}.
     * In particular, this method verifies if the description provided by the XML matches the description provided by
     * the authority code given in {@code <gml:identifier>} element, and reports discrepancies.
     * Note that this comparison between unmarshalled CRS and authoritative CRS is specific to this convenience method;
     * other APIs documented in <cite>see also</cite> section do not perform this comparison automatically.
     * Should the XML description and the authoritative description be in conflict, the XML description prevails
     * (see {@link #fromAuthority fromAuthority(…)} if a different behavior is needed).</p>
     *
     * @param  xml  coordinate reference system encoded in XML format.
     * @return the unmarshalled Coordinate Reference System.
     * @throws FactoryException if the object creation failed.
     *
     * @see org.apache.sis.referencing.factory.GeodeticObjectFactory#createFromXML(String)
     * @see org.apache.sis.xml.XML#unmarshal(String)
     *
     * @since 0.7
     */
    public static CoordinateReferenceSystem fromXML(final String xml) throws FactoryException {
        ArgumentChecks.ensureNonNull("text", xml);
        final CoordinateReferenceSystem crs = DefaultFactories.forBuildin(CRSFactory.class).createFromXML(xml);
        DefinitionVerifier.withAuthority(crs, Loggers.XML, CRS.class, "fromXML");
        return crs;
    }

    /**
     * Replaces the given coordinate reference system by an authoritative description, if one can be found.
     * This method can be invoked after constructing a CRS in a context where the EPSG (or other authority)
     * code is suspected more reliable than the rest of the description. A common case is a <cite>Well Known
     * Text</cite> (WKT) string declaring wrong projection method or parameter values for the EPSG code that
     * it pretends to describe. For example:
     *
     * <blockquote>
     *   {@code PROJCS["WGS 84 / Pseudo-Mercator",}<br>
     *   {@code   }(…base CRS omitted for brevity…)<br>
     *   {@code   PROJECTION["Mercator (variant A)"],} — <em><b>wrong:</b> shall be "Popular Visualisation Pseudo Mercator"</em><br>
     *   {@code   }(…parameters and axes omitted for brevity…)<br>
     *   {@code   AUTHORITY["EPSG", "3857"]]}
     * </blockquote>
     *
     * In such cases, Apache SIS behavior in {@link #fromWKT(String)}, {@link #fromXML(String)} and other methods is
     * conform to the <a href="http://docs.opengeospatial.org/is/12-063r5/12-063r5.html">ISO 19162 specification</a>:
     *
     * <blockquote><cite>"Should any attributes or values given in the cited identifier be in conflict with attributes
     * or values given explicitly in the WKT description, the WKT values shall prevail."</cite></blockquote>
     *
     * In situations where the opposite behavior is desired (i.e. to make the authority identifier prevails),
     * this method can be invoked. This method performs the following actions:
     *
     * <ul>
     *   <li>If the given CRS has an {@linkplain AbstractIdentifiedObject#getIdentifiers() identifier} and if the authority factory can
     *     {@linkplain org.apache.sis.referencing.factory.GeodeticAuthorityFactory#createCoordinateReferenceSystem(String) create a CRS}
     *     for that identifier, then:
     *     <ul>
     *       <li>If the CRS defined by the authority is {@linkplain Utilities#equalsIgnoreMetadata equal, ignoring metadata},
     *         to the given CRS, then this method returns silently the <em>authoritative</em> CRS.</li>
     *       <li>Otherwise if the CRS defined by the authority is equal, ignoring axis order and units, to the given CRS,
     *         then this method returns a <em>new</em> CRS derived from the authoritative one but with same
     *         {@linkplain org.apache.sis.referencing.cs.AxesConvention axes convention} than the given CRS.
     *         A warning is emitted.</li>
     *       <li>Otherwise this method discards the given CRS and returns the <em>authoritative</em> CRS.
     *         A warning is emitted with a message indicating where a difference has been found.</li>
     *     </ul>
     *   </li>
     *   <li>Otherwise if the given CRS does not have identifier, then this method
     *       {@linkplain org.apache.sis.referencing.factory.IdentifiedObjectFinder searches for an equivalent CRS}
     *       defined by the authority factory. If such CRS is found, then:
     *     <ul>
     *       <li>If the CRS defined by the authority is {@linkplain Utilities#equalsIgnoreMetadata equal, ignoring metadata},
     *         to the given CRS, then this method returns silently the <em>authoritative</em> CRS.</li>
     *       <li>Otherwise if the CRS defined by the authority is equal, ignoring axis order and units, to the given CRS,
     *         then this method returns silently a <em>new</em> CRS derived from the authoritative one but with same
     *         {@linkplain org.apache.sis.referencing.cs.AxesConvention axes convention} than the given CRS.</li>
     *     </ul>
     *   </li>
     *   <li>Otherwise this method silently returns the given CRS as-is.</li>
     * </ul>
     *
     * <div class="section">Avoiding warning redundancies</div>
     * The warnings logged by this method are redundant with warnings logged by other methods in this class,
     * in particular {@link #fromWKT(String)} and {@link #fromXML(String)} methods. For avoiding this annoyance,
     * a {@code null} value for the {@code warningFilter} argument means to shut off those redundant loggings.
     * A non-null {@code warningFilter} argument is more useful for CRS parsed by methods outside this class,
     * for example {@link org.apache.sis.io.wkt.WKTFormat} or {@link org.apache.sis.xml.XML#unmarshal(String)}.
     *
     * @param  crs            the CRS to replace by an authoritative CRS, or {@code null}.
     * @param  factory        the factory where to search for authoritative definitions, or {@code null} for the default.
     * @param  warningFilter  whether to log warnings, or {@code null} for the default behavior (which is to filter out
     *                        the warnings that are redundant with warnings emitted by other methods in this class).
     * @return the suggested CRS to use (may be the {@code crs} argument itself), or {@code null} if the given CRS was null.
     * @throws FactoryException if an error occurred while querying the authority factory.
     *
     * @since 1.0
     */
    public static CoordinateReferenceSystem fromAuthority(CoordinateReferenceSystem crs,
            final CRSAuthorityFactory factory, final Filter warningFilter) throws FactoryException
    {
        if (crs != null) {
            final DefinitionVerifier verification = DefinitionVerifier.withAuthority(crs, factory, true);
            if (verification != null) {
                crs = verification.authoritative;
                if (warningFilter != null) {
                    final LogRecord record = verification.warning(false);
                    if (record != null) {
                        record.setLoggerName(Modules.REFERENCING);
                        record.setSourceClassName(CRS.class.getName());
                        record.setSourceMethodName("fromAuthority");
                        if (warningFilter.isLoggable(record)) {
                            Logging.getLogger(Modules.REFERENCING).log(record);
                        }
                    }
                }
            }
        }
        return crs;
    }

    /**
     * Suggests a coordinate reference system which could be a common target for coordinate operations having the
     * given sources. This method compares the {@linkplain #getGeographicBoundingBox(CoordinateReferenceSystem)
     * domain of validity} of all given CRSs. If a CRS has a domain of validity that contains the domain of all other
     * CRS, then that CRS is returned. Otherwise this method verifies if a {@linkplain GeneralDerivedCRS#getBaseCRS()
     * base CRS} (usually a {@linkplain org.apache.sis.referencing.crs.DefaultGeographicCRS geographic CRS} instance)
     * would be suitable. If no suitable CRS is found, then this method returns {@code null}.
     *
     * <div class="note"><b>Use case:</b>
     * before to test if two arbitrary envelopes {@linkplain GeneralEnvelope#intersects(Envelope) intersect} each other,
     * they need to be {@linkplain Envelopes#transform(Envelope, CoordinateReferenceSystem) transformed} in the same CRS.
     * However if one CRS is a Transverse Mercator projection while the other CRS is a world-wide geographic CRS, then
     * attempts to use the Transverse Mercator projection as the common CRS is likely to fail since the geographic envelope
     * may span an area far outside the projection domain of validity. This {@code suggestCommonTarget(…)} method can used
     * for choosing a common CRS which is less likely to fail.</div>
     *
     * @param  regionOfInterest  the geographic area for which the coordinate operations will be applied,
     *                           or {@code null} if unknown. Will be intersected with CRS domains of validity.
     * @param  sourceCRS         the coordinate reference systems for which a common target CRS is desired.
     *                           May contain {@code null} elements, which are ignored.
     * @return a CRS that may be used as a common target for all the given source CRS in the given region of interest,
     *         or {@code null} if this method did not find a common target CRS. The returned CRS may be different than
     *         all given CRS.
     *
     * @since 0.8
     */
    public static CoordinateReferenceSystem suggestCommonTarget(GeographicBoundingBox regionOfInterest,
                                                                CoordinateReferenceSystem... sourceCRS)
    {
        CoordinateReferenceSystem bestCRS = null;
        /*
         * Compute the union of the domain of validity of all CRS. If a CRS does not specify a domain of validity,
         * then assume that the CRS is valid for the whole world if the CRS is geodetic (otherwise ignore that CRS).
         * Opportunistically remember the domain of validity of each CRS in this loop since we will need them later.
         */
        boolean worldwide = false;
        DefaultGeographicBoundingBox domain = null;
        final GeographicBoundingBox[] domains = new GeographicBoundingBox[sourceCRS.length];
        for (int i=0; i < sourceCRS.length; i++) {
            final CoordinateReferenceSystem crs = sourceCRS[i];
            final GeographicBoundingBox bbox = getGeographicBoundingBox(crs);
            if (bbox != null) {
                domains[i] = bbox;
                if (!worldwide) {
                    if (domain == null) {
                        domain = new DefaultGeographicBoundingBox(bbox);
                    } else {
                        domain.add(bbox);
                    }
                }
            } else if (crs instanceof GeodeticCRS) {
                /*
                 * Geodetic CRS (geographic or geocentric) can generally be presumed valid in a worldwide area.
                 * The 'worldwide' flag is a little optimization for remembering that we do not need to compute
                 * the union anymore, but we still need to continue the loop for fetching all bounding boxes.
                 */
                bestCRS = crs;                      // Fallback to be used if we don't find anything better.
                worldwide = true;
            }
        }
        /*
         * At this point we got the union of the domain of validity of all CRS. We are interested only in the
         * part that intersect the region of interest. If the union is whole world, we do not need to compute
         * the intersection; we can just leave the region of interest unchanged.
         */
        if (domain != null && !worldwide) {
            if (regionOfInterest != null) {
                domain.intersect(regionOfInterest);
            }
            regionOfInterest = domain;
            domain = null;
        }
        /*
         * Iterate again over the domain of validity of all CRS.  For each domain of validity, compute the area
         * which is inside the domain or interest and the area which is outside. The "best CRS" will be the one
         * which comply with the following rules, in preference order:
         *
         *   1) The CRS which is valid over the largest area of the region of interest.
         *   2) If two CRS are equally good according rule 1, then the CRS with the smallest "outside area".
         *
         * Example: given two source CRS, a geographic one and a projected one:
         *
         *   - If the projected CRS contains fully the region of interest, then it will be returned.
         *     The preference is given to the projected CRS because geometric operations are likely
         *     to be more accurate in that space. Furthermore forward conversions from geographic to
         *     projected CRS are usually faster than inverse conversions.
         *
         *   - Otherwise (i.e. if the region of interest is likely to be wider than the projected CRS
         *     domain of validity), then the geographic CRS will be returned.
         */
        final double roiArea  = Extents.area(regionOfInterest);   // NaN if 'regionOfInterest' is null.
        double maxInsideArea  = 0;
        double minOutsideArea = Double.POSITIVE_INFINITY;
        boolean tryDerivedCRS = false;
        do {
            for (int i=0; i < domains.length; i++) {
                final GeographicBoundingBox bbox = domains[i];
                if (bbox != null) {
                    double insideArea  = Extents.area(bbox);
                    double outsideArea = 0;
                    if (regionOfInterest != null) {
                        if (domain == null) {
                            domain = new DefaultGeographicBoundingBox(bbox);
                        } else {
                            domain.setBounds(bbox);
                        }
                        domain.intersect(regionOfInterest);
                        final double area = insideArea;
                        insideArea = Extents.area(domain);
                        outsideArea = area - insideArea;
                    }
                    if (insideArea > maxInsideArea || (insideArea == maxInsideArea && outsideArea < minOutsideArea)) {
                        maxInsideArea  = insideArea;
                        minOutsideArea = outsideArea;
                        bestCRS        = sourceCRS[i];
                    }
                }
            }
            /*
             * If the best CRS does not cover fully the region of interest, then we will redo the check again
             * but using base CRS instead. For example if the list of source CRS had some projected CRS, we
             * will try with the geographic CRS on which those projected CRS are based.
             */
            if (Double.isNaN(roiArea) || maxInsideArea < roiArea) {
                if (tryDerivedCRS) break;                                               // Do not try twice.
                final CoordinateReferenceSystem[] derivedCRS = new CoordinateReferenceSystem[sourceCRS.length];
                for (int i=0; i < derivedCRS.length; i++) {
                    GeographicBoundingBox bbox = null;
                    final CoordinateReferenceSystem crs = sourceCRS[i];
                    if (crs instanceof GeneralDerivedCRS) {
                        final CoordinateReferenceSystem baseCRS = ((GeneralDerivedCRS) crs).getBaseCRS();
                        bbox = getGeographicBoundingBox(baseCRS);
                        if (bbox == null && bestCRS == null && baseCRS instanceof GeodeticCRS) {
                            bestCRS = baseCRS;      // Fallback to be used if we don't find anything better.
                        }
                        tryDerivedCRS = true;
                        derivedCRS[i] = baseCRS;
                    }
                    domains[i] = bbox;
                }
                sourceCRS = derivedCRS;
            } else {
                break;
            }
        } while (tryDerivedCRS);
        return bestCRS;
    }

    /**
     * Finds a mathematical operation that transforms or converts coordinates from the given source to the
     * given target coordinate reference system. If an estimation of the geographic area containing the points
     * to transform is known, it can be specified for helping this method to find a better suited operation.
     * If no area of interest is specified, then the current default is the widest
     * {@linkplain AbstractCoordinateOperation#getDomainOfValidity() domain of validity}.
     * A future Apache SIS version may also take the country of current locale in account.
     *
     * <div class="note"><b>Note:</b>
     * the area of interest is just one aspect that may affect the coordinate operation.
     * Other aspects are the time of interest (because some coordinate operations take in account the
     * plate tectonics movement) or the desired accuracy. For more control on the coordinate operation
     * to create, see {@link CoordinateOperationContext}.</div>
     *
     * After the caller received a {@code CoordinateOperation} instance, the following methods can be invoked
     * for checking if the operation suits the caller's needs:
     *
     * <ul>
     *   <li>{@link #getGeographicBoundingBox(CoordinateOperation)}
     *       for checking if the operation is valid in the caller's area of interest.</li>
     *   <li>{@link #getLinearAccuracy(CoordinateOperation)}
     *       for checking if the operation has sufficient accuracy for caller's purpose.</li>
     * </ul>
     *
     * If the source and target CRS are equivalent, then this method returns an operation backed by an
     * {@linkplain org.apache.sis.referencing.operation.transform.AbstractMathTransform#isIdentity() identity}
     * transform. If there is no known operation between the given pair of CRS, then this method throws an
     * {@link OperationNotFoundException}.
     *
     * @param  sourceCRS       the CRS of source coordinates.
     * @param  targetCRS       the CRS of target coordinates.
     * @param  areaOfInterest  the area of interest, or {@code null} if none.
     * @return the mathematical operation from {@code sourceCRS} to {@code targetCRS}.
     * @throws OperationNotFoundException if no operation was found between the given pair of CRS.
     * @throws FactoryException if the operation can not be created for another reason.
     *
     * @see Envelopes#findOperation(Envelope, Envelope)
     * @see DefaultCoordinateOperationFactory#createOperation(CoordinateReferenceSystem, CoordinateReferenceSystem, CoordinateOperationContext)
     *
     * @since 0.7
     */
    public static CoordinateOperation findOperation(final CoordinateReferenceSystem sourceCRS,
                                                    final CoordinateReferenceSystem targetCRS,
                                                    final GeographicBoundingBox areaOfInterest)
            throws FactoryException
    {
        ArgumentChecks.ensureNonNull("sourceCRS", sourceCRS);
        ArgumentChecks.ensureNonNull("targetCRS", targetCRS);
        final CoordinateOperationContext context = CoordinateOperationContext.fromBoundingBox(areaOfInterest);
        /*
         * In principle following code should just delegate to factory.createOperation(…). However that operation
         * may fail if a connection to the EPSG database has been found, but the EPSG tables do not yet exist in
         * that database and we do not have the SQL scripts for creating them.
         */
        final DefaultCoordinateOperationFactory factory = CoordinateOperations.factory();
        try {
            return factory.createOperation(sourceCRS, targetCRS, context);
        } catch (UnavailableFactoryException e) {
            if (AuthorityFactories.failure(e)) {
                throw e;
            } else try {
                // Above method call replaced the EPSG factory by a fallback. Try again.
                return factory.createOperation(sourceCRS, targetCRS, context);
            } catch (FactoryException ex) {
                ex.addSuppressed(e);
                throw ex;
            }
        }
    }

    /**
     * Finds mathematical operations that transform or convert coordinates from the given source to the
     * given target coordinate reference system. If at least one operation exists, they are returned in
     * preference order: the operation having the widest intersection between its
     * {@linkplain AbstractCoordinateOperation#getDomainOfValidity() domain of validity}
     * and the given area of interest are returned first.
     *
     * @param  sourceCRS       the CRS of source coordinates.
     * @param  targetCRS       the CRS of target coordinates.
     * @param  areaOfInterest  the area of interest, or {@code null} if none.
     * @return mathematical operations from {@code sourceCRS} to {@code targetCRS}.
     * @throws OperationNotFoundException if no operation was found between the given pair of CRS.
     * @throws FactoryException if the operation can not be created for another reason.
     *
     * @see DefaultCoordinateOperationFactory#createOperations(CoordinateReferenceSystem, CoordinateReferenceSystem, CoordinateOperationContext)
     *
     * @since 1.0
     */
    public static List<CoordinateOperation> findOperations(final CoordinateReferenceSystem sourceCRS,
                                                           final CoordinateReferenceSystem targetCRS,
                                                           final GeographicBoundingBox areaOfInterest)
            throws FactoryException
    {
        ArgumentChecks.ensureNonNull("sourceCRS", sourceCRS);
        ArgumentChecks.ensureNonNull("targetCRS", targetCRS);
        final CoordinateOperationContext context = CoordinateOperationContext.fromBoundingBox(areaOfInterest);
        final DefaultCoordinateOperationFactory factory = CoordinateOperations.factory();
        try {
            return factory.createOperations(sourceCRS, targetCRS, context);
        } catch (UnavailableFactoryException e) {
            if (AuthorityFactories.failure(e)) {
                throw e;
            } else try {
                return Collections.singletonList(factory.createOperation(sourceCRS, targetCRS, context));
            } catch (FactoryException ex) {
                ex.addSuppressed(e);
                throw ex;
            }
        }
    }

    /**
     * Returns a positional accuracy estimation in metres for the given operation, or {@code NaN} if unknown.
     * This method applies the following heuristics:
     *
     * <ul>
     *   <li>If the given operation is an instance of {@link AbstractCoordinateOperation}, then delegate to the
     *       operation {@link AbstractCoordinateOperation#getLinearAccuracy() getLinearAccuracy()} method.</li>
     *
     *   <li>Otherwise if at least one {@linkplain org.apache.sis.metadata.iso.quality.DefaultQuantitativeResult
     *       quantitative result} is found with a linear unit, then return the largest value converted to metres.</li>
     *
     *   <li>Otherwise if the operation is a {@linkplain org.apache.sis.referencing.operation.DefaultConversion
     *       conversion}, then returns 0 since a conversion is by definition accurate up to rounding errors.</li>
     *
     *   <li>Otherwise if the operation is a {@linkplain org.apache.sis.referencing.operation.DefaultTransformation
     *       transformation}, then the returned value depends on whether the datum shift were applied with the help
     *       of Bursa-Wolf parameters of not.</li>
     * </ul>
     *
     * See {@link AbstractCoordinateOperation#getLinearAccuracy()} for more details on the above heuristic rules.
     *
     * @param  operation  the coordinate operation for which to get the accuracy estimation, or {@code null}.
     * @return the accuracy estimation (always in meters), or NaN if unknown.
     *
     * @see #findOperation(CoordinateReferenceSystem, CoordinateReferenceSystem, GeographicBoundingBox)
     *
     * @since 0.7
     */
    public static double getLinearAccuracy(final CoordinateOperation operation) {
        if (operation == null) {
            return Double.NaN;
        } else if (operation instanceof AbstractCoordinateOperation) {
            return ((AbstractCoordinateOperation) operation).getLinearAccuracy();
        } else {
            return PositionalAccuracyConstant.getLinearAccuracy(operation);
        }
    }

    /**
     * Returns the valid geographic area for the given coordinate operation, or {@code null} if unknown.
     * This method explores the {@linkplain AbstractCoordinateOperation#getDomainOfValidity() domain of validity}
     * associated with the given operation. If more than one geographic bounding box is found, then this method
     * computes their {@linkplain DefaultGeographicBoundingBox#add(GeographicBoundingBox) union}.
     *
     * <p><b>Fallback:</b> if the given operation does not declare explicitly a domain of validity, then this
     * method computes the intersection of the domain of validity declared by source and target CRS. If no CRS
     * declare a domain of validity, then this method returns {@code null}.</p>
     *
     * @param  operation  the coordinate operation for which to get the domain of validity, or {@code null}.
     * @return the geographic area where the operation is valid, or {@code null} if unspecified.
     *
     * @see #findOperation(CoordinateReferenceSystem, CoordinateReferenceSystem, GeographicBoundingBox)
     * @see Extents#getGeographicBoundingBox(Extent)
     *
     * @category information
     *
     * @since 0.7
     */
    public static GeographicBoundingBox getGeographicBoundingBox(final CoordinateOperation operation) {
        if (operation == null) {
            return null;
        }
        GeographicBoundingBox bbox = Extents.getGeographicBoundingBox(operation.getDomainOfValidity());
        if (bbox == null) {
            bbox = Extents.intersection(getGeographicBoundingBox(operation.getSourceCRS()),
                                        getGeographicBoundingBox(operation.getTargetCRS()));
        }
        return bbox;
    }

    /**
     * Returns the valid geographic area for the given coordinate reference system, or {@code null} if unknown.
     * This method explores the {@linkplain org.apache.sis.referencing.crs.AbstractCRS#getDomainOfValidity() domain of
     * validity} associated with the given CRS. If more than one geographic bounding box is found, then this method
     * computes their {@linkplain DefaultGeographicBoundingBox#add(GeographicBoundingBox) union}.
     * together.
     *
     * @param  crs  the coordinate reference system for which to get the domain of validity, or {@code null}.
     * @return the geographic area where the coordinate reference system is valid, or {@code null} if unspecified.
     *
     * @see #getDomainOfValidity(CoordinateReferenceSystem)
     * @see Extents#getGeographicBoundingBox(Extent)
     *
     * @category information
     */
    public static GeographicBoundingBox getGeographicBoundingBox(final CoordinateReferenceSystem crs) {
        return (crs != null) ? Extents.getGeographicBoundingBox(crs.getDomainOfValidity()) : null;
    }

    /**
     * Returns the domain of validity of the specified coordinate reference system, or {@code null} if unknown.
     * If non-null, then the returned envelope will use the same coordinate reference system them the given CRS
     * argument.
     *
     * @param  crs  the coordinate reference system, or {@code null}.
     * @return the envelope with coordinates in the given CRS, or {@code null} if none.
     *
     * @see #getGeographicBoundingBox(CoordinateReferenceSystem)
     *
     * @category information
     * @since 0.8
     */
    public static Envelope getDomainOfValidity(final CoordinateReferenceSystem crs) {
        Envelope envelope = null;
        GeneralEnvelope merged = null;
        /* if (envelope == null) */ {   // Condition needed on other branches but not on trunk.
            final GeographicBoundingBox bounds = getGeographicBoundingBox(crs);
            if (bounds != null && !Boolean.FALSE.equals(bounds.getInclusion())) {
                /*
                 * We do not assign WGS84 unconditionally to the geographic bounding box, because
                 * it is not defined to be on a particular datum; it is only approximated bounds.
                 * We try to get the GeographicCRS from the user-supplied CRS in order to reduce
                 * the amount of transformation needed.
                 */
                final SingleCRS targetCRS = getHorizontalComponent(crs);
                final GeographicCRS sourceCRS = ReferencingUtilities.toNormalizedGeographicCRS(targetCRS, false, false);
                if (sourceCRS != null) {
                    envelope = merged = new GeneralEnvelope(bounds);
                    merged.translate(-getGreenwichLongitude(sourceCRS), 0);
                    merged.setCoordinateReferenceSystem(sourceCRS);
                    try {
                        envelope = Envelopes.transform(envelope, targetCRS);
                    } catch (TransformException exception) {
                        /*
                         * The envelope is probably outside the range of validity for this CRS.
                         * It should not occurs, since the envelope is supposed to describe the
                         * CRS area of validity. Logs a warning and returns null, since it is a
                         * legal return value according this method contract.
                         */
                        unexpectedException("getEnvelope", exception);
                        envelope = null;
                    }
                }
            }
        }
        return envelope;
    }

    /**
     * Creates a compound coordinate reference system from an ordered list of CRS components.
     * A CRS is inferred from the given components and the domain of validity is set to the
     * {@linkplain org.apache.sis.metadata.iso.extent.DefaultExtent#intersect intersection}
     * of the domain of validity of all components.
     *
     * <div class="section">Ellipsoidal height</div>
     * If a two-dimensional geographic or projected CRS if followed or preceded by a vertical CRS with ellipsoidal
     * {@linkplain org.apache.sis.referencing.datum.DefaultVerticalDatum#getVerticalDatumType() datum type}, then
     * this method combines them in a single three-dimensional geographic or projected CRS.  Note that standalone
     * ellipsoidal heights are not allowed according ISO 19111. But if such situation is nevertheless found, then
     * the action described here fixes the issue. This is the reverse of <code>{@linkplain #getVerticalComponent
     * getVerticalComponent}(crs, true)</code>.
     *
     * <div class="section">Components order</div>
     * Apache SIS is permissive on the order of components that can be used in a compound CRS.
     * However for better inter-operability, users are encouraged to follow the order mandated by ISO 19162:
     *
     * <ol>
     *   <li>A mandatory horizontal CRS (only one of two-dimensional {@code GeographicCRS} or {@code ProjectedCRS} or {@code EngineeringCRS}).</li>
     *   <li>Optionally followed by a {@code VerticalCRS} or a {@code ParametricCRS} (but not both).</li>
     *   <li>Optionally followed by a {@code TemporalCRS}.</li>
     * </ol>
     *
     * @param  components  the sequence of coordinate reference systems making the compound CRS.
     * @return the compound CRS, or {@code components[0]} if the given array contains only one component.
     * @throws IllegalArgumentException if the given array is empty or if the array contains incompatible components.
     * @throws FactoryException if the geodetic factory failed to create the compound CRS.
     *
     * @since 0.8
     *
     * @see org.apache.sis.referencing.crs.DefaultCompoundCRS
     * @see GeodeticObjectFactory#createCompoundCRS(Map, CoordinateReferenceSystem...)
     * @see org.apache.sis.geometry.Envelopes#compound(Envelope...)
     * @see org.apache.sis.referencing.operation.transform.MathTransforms#compound(MathTransform...)
     */
    public static CoordinateReferenceSystem compound(final CoordinateReferenceSystem... components) throws FactoryException {
        ArgumentChecks.ensureNonNull("components", components);
        switch (components.length) {
            case 0: {
                throw new IllegalArgumentException(Errors.format(Errors.Keys.EmptyArgument_1, "components"));
            }
            case 1: {
                final CoordinateReferenceSystem crs = components[0];
                if (crs != null) return crs;
                break;
            }
        }
        return new EllipsoidalHeightCombiner().createCompoundCRS(components);
    }

    /**
     * Gets or creates a coordinate reference system with a subset of the dimensions of the given CRS.
     * This method can be used for dimensionality reduction, but not for changing axis order.
     * The specified dimensions are used as if they were in strictly increasing order without duplicated values.
     *
     * <div class="section">Ellipsoidal height</div>
     * This method can transform a three-dimensional geographic CRS into a two-dimensional geographic CRS.
     * In this aspect, this method is the converse of {@link #compound(CoordinateReferenceSystem...)}.
     * This method can also extract the {@linkplain CommonCRS.Vertical#ELLIPSOIDAL ellipsoidal height}
     * from a three-dimensional geographic CRS, but this is generally not recommended since ellipsoidal
     * heights make little sense without their (<var>latitude</var>, <var>longitude</var>) locations.
     *
     * @param  crs         the CRS to reduce the dimensionality, or {@code null} if none.
     * @param  dimensions  the dimensions to retain. The dimensions will be taken in increasing order, ignoring duplicated values.
     * @return a coordinate reference system for the given dimensions. May be the given {@code crs}, which may be {@code null}.
     * @throws IllegalArgumentException if the given array is empty or if the array contains invalid indices.
     * @throws FactoryException if the geodetic factory failed to create a compound CRS.
     *
     * @see #getComponentAt(CoordinateReferenceSystem, int, int)
     * @see #compound(CoordinateReferenceSystem...)
     *
     * @since 1.0
     */
    public static CoordinateReferenceSystem reduce(final CoordinateReferenceSystem crs, final int... dimensions) throws FactoryException {
        ArgumentChecks.ensureNonNull("dimensions", dimensions);
        if (crs == null) {
            return null;
        }
        final int dimension = ReferencingUtilities.getDimension(crs);
        long selected = 0;
        for (final int d : dimensions) {
            if (d < 0 || d >= dimension) {
                throw new IndexOutOfBoundsException(Errors.format(Errors.Keys.IndexOutOfBounds_1, d));
            }
            if (d >= Long.SIZE) {
                throw new IllegalArgumentException(Errors.format(Errors.Keys.ExcessiveNumberOfDimensions_1, d));
            }
            selected |= (1L << d);
        }
        if (selected == 0) {
            throw new IllegalArgumentException(Errors.format(Errors.Keys.EmptyArgument_1, "dimensions"));
        }
        final List<CoordinateReferenceSystem> components = new ArrayList<>(Long.bitCount(selected));
        reduce(0, crs, dimension, selected, components);
        return compound(components.toArray(new CoordinateReferenceSystem[components.size()]));
    }

    /**
     * Adds the components of reduced CRS into the given list.
     * This method may invoke itself recursively for walking through compound CRS.
     *
     * @param  previous    number of dimensions of previous CRS.
     * @param  crs         the CRS for which to select components.
     * @param  dimension   number of dimensions of {@code crs}.
     * @param  selected    bitmask of dimensions to select.
     * @param  addTo       where to add CRS components.
     * @return new bitmask after removal of dimensions of the components added to {@code addTo}.
     */
    private static long reduce(int previous, final CoordinateReferenceSystem crs, int dimension, long selected,
            final List<CoordinateReferenceSystem> addTo) throws FactoryException
    {
        final long current = (Numerics.bitmask(dimension) - 1) << previous;
        final long intersect = selected & current;
        if (intersect != 0) {
            if (intersect == current) {
                addTo.add(crs);
                selected &= ~current;
            } else if (crs instanceof CompoundCRS) {
                for (final CoordinateReferenceSystem component : ((CompoundCRS) crs).getComponents()) {
                    dimension = ReferencingUtilities.getDimension(component);
                    selected = reduce(previous, component, dimension, selected, addTo);
                    if ((selected & current) == 0) break;           // Stop if it would be useless to continue.
                    previous += dimension;
                }
            } else if (dimension == 3 && crs instanceof SingleCRS) {
                final Datum datum = ((SingleCRS) crs).getDatum();
                if (datum instanceof GeodeticDatum) {
                    final boolean isVertical = Long.bitCount(intersect) == 1;               // Presumed for now, verified later.
                    final int verticalDimension = Long.numberOfTrailingZeros((isVertical ? intersect : ~intersect) >>> previous);
                    final CoordinateSystemAxis verticalAxis = crs.getCoordinateSystem().getAxis(verticalDimension);
                    if (AxisDirections.isVertical(verticalAxis.getDirection())) try {
                        addTo.add(new EllipsoidalHeightSeparator((GeodeticDatum) datum).separate((SingleCRS) crs, isVertical));
                        selected &= ~current;
                    } catch (IllegalArgumentException | ClassCastException e) {
                        throw new FactoryException(Resources.format(Resources.Keys.CanNotSeparateCRS_1, crs.getName()));
                    }
                }
            }
        }
        if ((selected & current) != 0) {
            throw new FactoryException(Resources.format(Resources.Keys.CanNotSeparateCRS_1, crs.getName()));
        }
        return selected;
    }

    /**
     * Returns {@code true} if the given CRS is horizontal. The current implementation considers a
     * CRS as horizontal if it is two-dimensional and comply with one of the following conditions:
     *
     * <ul>
     *   <li>is an instance of {@link GeographicCRS} (or an equivalent {@link GeodeticCRS}), or</li>
     *   <li>is an instance of {@link ProjectedCRS}, or</li>
     *   <li>is an instance of {@link EngineeringCRS} (following
     *     <a href="http://docs.opengeospatial.org/is/12-063r5/12-063r5.html#111">ISO 19162 §16.1</a>
     *     definition of {@literal <horizontal crs>}).</li>
     * </ul>
     *
     * In case of doubt, this method conservatively returns {@code false}.
     *
     * @todo Future SIS implementation may extend the above conditions list. For example a radar station could
     *       use a polar coordinate system in a <code>DerivedCRS</code> instance based on a projected CRS.
     *       Conversely, a future SIS versions may impose more conditions on <code>EngineeringCRS</code>.
     *       See <a href="http://issues.apache.org/jira/browse/SIS-161">SIS-161</a>.
     *
     * @param  crs  the coordinate reference system, or {@code null}.
     * @return {@code true} if the given CRS is non-null and likely horizontal, or {@code false} otherwise.
     *
     * @see #getHorizontalComponent(CoordinateReferenceSystem)
     *
     * @category information
     */
    public static boolean isHorizontalCRS(final CoordinateReferenceSystem crs) {
        return horizontalCode(crs) == 2;
    }

    /**
     * If the given CRS would qualify as horizontal except for its number of dimensions, returns that number.
     * Otherwise returns 0. The number of dimensions can only be 2 or 3.
     */
    private static int horizontalCode(final CoordinateReferenceSystem crs) {
        /*
         * In order to determine if the CRS is geographic, checking the CoordinateSystem type is more reliable
         * then checking if the CRS implements the GeographicCRS interface.  This is because the GeographicCRS
         * type did not existed in ISO 19111:2007, so a CRS could be standard-compliant without implementing
         * the GeographicCRS interface.
         */
        boolean isEngineering = false;
        final boolean isGeodetic = (crs instanceof GeodeticCRS);
        if (isGeodetic || crs instanceof ProjectedCRS || (isEngineering = (crs instanceof EngineeringCRS))) {
            final CoordinateSystem cs = crs.getCoordinateSystem();
            final int dim = cs.getDimension();
            if ((dim & ~1) == 2 && (!isGeodetic || (cs instanceof EllipsoidalCS))) {
                if (isEngineering) {
                    int n = 0;
                    for (int i=0; i<dim; i++) {
                        if (AxisDirections.isCompass(cs.getAxis(i).getDirection())) n++;
                    }
                    // If we don't have exactly 2 east, north, etc. directions, consider as non-horizontal.
                    if (n != 2) return 0;
                }
                return dim;
            }
        }
        return 0;
    }

    /**
     * Returns the first horizontal coordinate reference system found in the given CRS, or {@code null} if there is
     * none. If the given CRS is already horizontal according {@link #isHorizontalCRS(CoordinateReferenceSystem)},
     * then this method returns it as-is. Otherwise if the given CRS is compound, then this method searches for the
     * first horizontal component in the order of the {@linkplain #getSingleComponents(CoordinateReferenceSystem)
     * single components list}.
     *
     * <p>In the special case where a three-dimensional geographic or projected CRS is found, this method
     * will create a two-dimensional geographic or projected CRS without the vertical axis.</p>
     *
     * @param  crs  the coordinate reference system, or {@code null}.
     * @return the first horizontal CRS, or {@code null} if none.
     *
     * @category information
     */
    public static SingleCRS getHorizontalComponent(final CoordinateReferenceSystem crs) {
        switch (horizontalCode(crs)) {
            /*
             * If the CRS is already two-dimensional and horizontal, return as-is.
             * We don't need to check if crs is an instance of SingleCRS since all
             * CRS accepted by horizontalCode(…) are SingleCRS.
             */
            case 2: {
                return (SingleCRS) crs;
            }
            case 3: {
                /*
                 * The CRS would be horizontal if we can remove the vertical axis. CoordinateSystems.replaceAxes(…)
                 * will do this task for us. We can verify if the operation has been successful by checking that
                 * the number of dimensions has been reduced by 1 (from 3 to 2).
                 */
                final CoordinateSystem cs = CoordinateSystems.replaceAxes(crs.getCoordinateSystem(), new AxisFilter() {
                    @Override public boolean accept(final CoordinateSystemAxis axis) {
                        return !AxisDirections.isVertical(axis.getDirection());
                    }
                });
                if (cs.getDimension() != 2) break;
                /*
                 * Most of the time, the CRS to rebuild will be geodetic. In such case we known that the
                 * coordinate system is ellipsoidal because (i.e. the CRS is geographic) because it was
                 * a condition verified by horizontalCode(…). A ClassCastException would be a bug.
                 */
                final Map<String, ?> properties = ReferencingUtilities.getPropertiesForModifiedCRS(crs);
                if (crs instanceof GeodeticCRS) {
                    return new DefaultGeographicCRS(properties, ((GeodeticCRS) crs).getDatum(), (EllipsoidalCS) cs);
                }
                /*
                 * In Apache SIS implementation, the Conversion contains the source and target CRS together with
                 * a MathTransform.   We need to recreate the same conversion, but without CRS and MathTransform
                 * for letting SIS create or associate new ones, which will be two-dimensional now.
                 */
                if (crs instanceof ProjectedCRS) {
                    final ProjectedCRS  proj = (ProjectedCRS) crs;
                    final GeographicCRS base = (GeographicCRS) getHorizontalComponent(proj.getBaseCRS());
                    Conversion fromBase = proj.getConversionFromBase();
                    fromBase = new DefaultConversion(IdentifiedObjects.getProperties(fromBase),
                            fromBase.getMethod(), null, fromBase.getParameterValues());
                    return new DefaultProjectedCRS(properties, base, fromBase, (CartesianCS) cs);
                }
                /*
                 * If the CRS is neither geographic or projected, then it is engineering.
                 */
                return new DefaultEngineeringCRS(properties, ((EngineeringCRS) crs).getDatum(), cs);
            }
        }
        if (crs instanceof CompoundCRS) {
            final CompoundCRS cp = (CompoundCRS) crs;
            for (final CoordinateReferenceSystem c : cp.getComponents()) {
                final SingleCRS candidate = getHorizontalComponent(c);
                if (candidate != null) {
                    return candidate;
                }
            }
        }
        return null;
    }

    /**
     * Returns the first vertical coordinate reference system found in the given CRS, or {@code null} if there is none.
     * If the given CRS is already an instance of {@code VerticalCRS}, then this method returns it as-is.
     * Otherwise if the given CRS is compound, then this method searches for the first vertical component
     * in the order of the {@linkplain #getSingleComponents(CoordinateReferenceSystem) single components list}.
     *
     * <div class="section">Height in a three-dimensional geographic CRS</div>
     * In ISO 19111 model, ellipsoidal heights are indissociable from geographic CRS because such heights
     * without their (<var>latitude</var>, <var>longitude</var>) locations make little sense. Consequently
     * a standard-conformant library should return {@code null} when asked for the {@code VerticalCRS}
     * component of a geographic CRS. This is what {@code getVerticalComponent(…)} does when the
     * {@code allowCreateEllipsoidal} argument is {@code false}.
     *
     * <p>However in some exceptional cases, handling ellipsoidal heights like any other kind of heights
     * may simplify the task. For example when computing <em>difference</em> between heights above the
     * same datum, the impact of ignoring locations may be smaller (but not necessarily canceled).
     * Orphan {@code VerticalCRS} may also be useful for information purpose like labeling a plot axis.
     * If the caller feels confident that ellipsoidal heights are safe for his task, he can set the
     * {@code allowCreateEllipsoidal} argument to {@code true}. In such case, this {@code getVerticalComponent(…)}
     * method will create a temporary {@code VerticalCRS} from the first three-dimensional {@code GeographicCRS}
     * <em>in last resort</em>, only if it can not find an existing {@code VerticalCRS} instance.
     * <strong>Note that this is not a valid CRS according ISO 19111</strong> — use with care.</p>
     *
     * @param  crs  the coordinate reference system, or {@code null}.
     * @param  allowCreateEllipsoidal {@code true} for allowing the creation of orphan CRS for ellipsoidal heights.
     *         The recommended value is {@code false}.
     * @return the first vertical CRS, or {@code null} if none.
     *
     * @see #compound(CoordinateReferenceSystem...)
     *
     * @category information
     */
    public static VerticalCRS getVerticalComponent(final CoordinateReferenceSystem crs, final boolean allowCreateEllipsoidal) {
        if (crs instanceof VerticalCRS) {
            return (VerticalCRS) crs;
        }
        if (crs instanceof CompoundCRS) {
            final CompoundCRS cp = (CompoundCRS) crs;
            boolean a = false;
            do { // Executed at most twice.
                for (final CoordinateReferenceSystem c : cp.getComponents()) {
                    final VerticalCRS candidate = getVerticalComponent(c, a);
                    if (candidate != null) {
                        return candidate;
                    }
                }
            } while ((a = !a) == allowCreateEllipsoidal);
        }
        if (allowCreateEllipsoidal && horizontalCode(crs) == 3) {
            final CoordinateSystem cs = crs.getCoordinateSystem();
            final int i = AxisDirections.indexOfColinear(cs, AxisDirection.UP);
            if (i >= 0) {
                final CoordinateSystemAxis axis = cs.getAxis(i);
                VerticalCRS c = CommonCRS.Vertical.ELLIPSOIDAL.crs();
                if (!c.getCoordinateSystem().getAxis(0).equals(axis)) {
                    final Map<String,?> properties = IdentifiedObjects.getProperties(c);
                    c = new DefaultVerticalCRS(properties, c.getDatum(), new DefaultVerticalCS(properties, axis));
                }
                return c;
            }
        }
        return null;
    }

    /**
     * Returns the first temporal coordinate reference system found in the given CRS, or {@code null} if there is none.
     * If the given CRS is already an instance of {@code TemporalCRS}, then this method returns it as-is.
     * Otherwise if the given CRS is compound, then this method searches for the first temporal component
     * in the order of the {@linkplain #getSingleComponents(CoordinateReferenceSystem) single components list}.
     *
     * @param  crs  the coordinate reference system, or {@code null}.
     * @return the first temporal CRS, or {@code null} if none.
     *
     * @category information
     */
    public static TemporalCRS getTemporalComponent(final CoordinateReferenceSystem crs) {
        if (crs instanceof TemporalCRS) {
            return (TemporalCRS) crs;
        }
        if (crs instanceof CompoundCRS) {
            final CompoundCRS cp = (CompoundCRS) crs;
            for (final CoordinateReferenceSystem c : cp.getComponents()) {
                final TemporalCRS candidate = getTemporalComponent(c);
                if (candidate != null) {
                    return candidate;
                }
            }
        }
        return null;
    }

    /**
     * Returns the ordered list of single coordinate reference systems for the specified CRS.
     * This method performs the following choices:
     *
     * <ul>
     *   <li>If the given CRS is null, returns an empty list.</li>
     *   <li>If the given CRS is an instance of {@link SingleCRS}, returns that instance in a singleton list.</li>
     *   <li>If the given CRS is an instance of {@link CompoundCRS}, returns a flattened list of its
     *       {@linkplain DefaultCompoundCRS#getComponents() components}. Some components may themselves be
     *       other {@code CompoundCRS} instances, in which case those compound CRS are also flattened in their
     *       list of {@code SingleCRS} components.</li>
     *   <li>Otherwise throws a {@code ClassCastException}.</li>
     * </ul>
     *
     * <div class="note"><b>Example:</b>
     * Apache SIS allows 4-dimensional (<var>x</var>,<var>y</var>,<var>z</var>,<var>t</var>)
     * coordinate reference system to be built in two different ways as shown below:
     *
     * <div class="horizontal-flow">
     * <div><p><b>Hierarchical structure</b></p>
     * <blockquote>
     *   <code>CompoundCRS</code> — (<var>x</var>, <var>y</var>, <var>z</var>, <var>t</var>)<br>
     *   <code>  ├─CompoundCRS</code> — (<var>x</var>, <var>y</var>, <var>z</var>)<br>
     *   <code>  │   ├─ProjectedCRS</code> — (<var>x</var>, <var>y</var>)<br>
     *   <code>  │   └─VerticalCRS</code> — (<var>z</var>)<br>
     *   <code>  └─TemporalCRS</code> — (<var>t</var>)
     * </blockquote></div>
     * <div><p><b>Flat list</b></p>
     * <blockquote>
     *   <code>CompoundCRS</code> — (<var>x</var>, <var>y</var>, <var>z</var>, <var>t</var>)<br>
     *   <code>  ├─ProjectedCRS</code> — (<var>x</var>, <var>y</var>)<br>
     *   <code>  ├─VerticalCRS</code> — (<var>z</var>)<br>
     *   <code>  └─TemporalCRS</code> — (<var>t</var>)
     * </blockquote>
     * </div></div>
     *
     * This method guaranteed that the returned list is a flat one as shown on the right side.
     * Note that such flat lists are the only one allowed by ISO/OGC standards for compound CRS.
     * The hierarchical structure is an Apache SIS flexibility.</div>
     *
     * @param  crs  the coordinate reference system, or {@code null}.
     * @return the single coordinate reference systems, or an empty list if the given CRS is {@code null}.
     * @throws ClassCastException if a CRS is neither a {@link SingleCRS} or a {@link CompoundCRS}.
     *
     * @see DefaultCompoundCRS#getSingleComponents()
     */
    public static List<SingleCRS> getSingleComponents(final CoordinateReferenceSystem crs) {
        final List<SingleCRS> singles;
        if (crs == null) {
            singles = Collections.emptyList();
        } else if (crs instanceof CompoundCRS) {
            if (crs instanceof DefaultCompoundCRS) {
                singles = ((DefaultCompoundCRS) crs).getSingleComponents();
            } else {
                final List<CoordinateReferenceSystem> elements = ((CompoundCRS) crs).getComponents();
                singles = new ArrayList<>(elements.size());
                ReferencingUtilities.getSingleComponents(elements, singles);
            }
        } else {
            // Intentional CassCastException here if the crs is not a SingleCRS.
            singles = Collections.singletonList((SingleCRS) crs);
        }
        return singles;
    }

    /**
     * Returns the coordinate reference system in the given range of dimension indices.
     * This method processes as below:
     *
     * <ul>
     *   <li>If the given {@code crs} is {@code null}, then this method returns {@code null}.</li>
     *   <li>Otherwise if {@code lower} is 0 and {@code upper} is the number of CRS dimensions,
     *       then this method returns the given CRS unchanged.</li>
     *   <li>Otherwise if the given CRS is an instance of {@link CompoundCRS}, then this method
     *       searches for a {@linkplain CompoundCRS#getComponents() component} where:
     *       <ul>
     *         <li>The {@linkplain org.apache.sis.referencing.cs.AbstractCS#getDimension() number of dimensions}
     *             is equals to {@code upper - lower};</li>
     *         <li>The sum of the number of dimensions of all previous CRS is equals to {@code lower}.</li>
     *       </ul>
     *       If such component is found, then it is returned.</li>
     *   <li>Otherwise (i.e. no component match), this method returns {@code null}.</li>
     * </ul>
     *
     * This method does <strong>not</strong> build new CRS from the components. For example this method does not
     * create a {@link CompoundCRS} or a three-dimensional CRS if the given range spans more than one component.
     *
     * @param  crs    the coordinate reference system to decompose, or {@code null}.
     * @param  lower  the first dimension to keep, inclusive.
     * @param  upper  the last  dimension to keep, exclusive.
     * @return the sub-coordinate system, or {@code null} if the given {@code crs} was {@code null}
     *         or can not be decomposed for dimensions in the [{@code lower} … {@code upper}] range.
     * @throws IndexOutOfBoundsException if the given index are out of bounds.
     *
     * @see #reduce(CoordinateReferenceSystem, int...)
     * @see org.apache.sis.geometry.GeneralEnvelope#subEnvelope(int, int)
     *
     * @since 0.5
     */
    public static CoordinateReferenceSystem getComponentAt(CoordinateReferenceSystem crs, int lower, int upper) {
        int dimension = ReferencingUtilities.getDimension(crs);
        ArgumentChecks.ensureValidIndexRange(dimension, lower, upper);
check:  while (lower != 0 || upper != dimension) {
            if (crs instanceof CompoundCRS) {
                // We need nested CompoundCRS (if any) below, not a flattened list of SingleCRS.
                final List<CoordinateReferenceSystem> components = ((CompoundCRS) crs).getComponents();
                final int size = components.size();
                for (int i=0; i<size; i++) {
                    crs = components.get(i);
                    dimension = crs.getCoordinateSystem().getDimension();
                    if (lower < dimension) {
                        /*
                         * The requested dimensions may intersect the dimension of this CRS.
                         * The outer loop will perform the verification, and eventually go
                         * down again in the tree of sub-components.
                         */
                        continue check;
                    }
                    lower -= dimension;
                    upper -= dimension;
                }
            }
            return null;
        }
        return crs;
    }

    /**
     * Returns the Greenwich longitude of the prime meridian of the given CRS in degrees.
     * If the prime meridian uses an other unit than degrees, then the value will be converted.
     *
     * @param  crs  the coordinate reference system from which to get the prime meridian.
     * @return the Greenwich longitude (in degrees) of the prime meridian of the given CRS.
     *
     * @see org.apache.sis.referencing.datum.DefaultPrimeMeridian#getGreenwichLongitude(Unit)
     *
     * @since 0.5
     */
    public static double getGreenwichLongitude(final GeodeticCRS crs) {
        ArgumentChecks.ensureNonNull("crs", crs);
        return ReferencingUtilities.getGreenwichLongitude(crs.getDatum().getPrimeMeridian(), Units.DEGREE);
    }

    /**
     * Returns the system-wide authority factory used by {@link #forCode(String)} and other SIS methods.
     * If the given authority is non-null, then this method returns a factory specifically for that authority.
     * Otherwise, this method returns the {@link org.apache.sis.referencing.factory.MultiAuthoritiesFactory}
     * instance that manages all other factories.
     *
     * <p>The {@code authority} argument can be {@code "EPSG"}, {@code "OGC"} or any other authority found
     * on the classpath. In the {@code "EPSG"} case, whether the full set of EPSG codes is supported or not
     * depends on whether a {@linkplain org.apache.sis.referencing.factory.sql connection to the database}
     * can be established. If no connection can be established, then this method returns a small embedded
     * EPSG factory containing at least the CRS defined in the {@link #forCode(String)} method javadoc.</p>
     *
     * <p>User-defined authorities can be added to the SIS environment by creating a {@code CRSAuthorityFactory}
     * implementation with a public no-argument constructor, and declaring the fully-qualified name of that class
     * in a file at the following location:</p>
     *
     * {@preformat text
     *     META-INF/services/org.opengis.referencing.crs.CRSAuthorityFactory
     * }
     *
     * @param  authority  the authority of the desired factory (typically {@code "EPSG"} or {@code "OGC"}),
     *         or {@code null} for the {@link org.apache.sis.referencing.factory.MultiAuthoritiesFactory}
     *         instance that manage all factories.
     * @return the system-wide authority factory used by SIS for the given authority.
     * @throws FactoryException if no factory can be returned for the given authority.
     *
     * @see #forCode(String)
     * @see org.apache.sis.referencing.factory.MultiAuthoritiesFactory
     *
     * @since 0.7
     */
    public static CRSAuthorityFactory getAuthorityFactory(final String authority) throws FactoryException {
        if (authority == null) {
            return AuthorityFactories.ALL;
        }
        return AuthorityFactories.ALL.getAuthorityFactory(CRSAuthorityFactory.class, authority, null);
    }

    /**
     * Invoked when an unexpected exception occurred. Those exceptions must be non-fatal, i.e. the caller
     * <strong>must</strong> have a reasonable fallback (otherwise it should propagate the exception).
     */
    private static void unexpectedException(final String methodName, final Exception exception) {
        Logging.unexpectedException(Logging.getLogger(Modules.REFERENCING), CRS.class, methodName, exception);
    }
}