Given a source coordinate reference system (CRS) in which existing coordinate values are expressed, and a target coordinate reference system in which coordinate values are desired, Apache SIS can provide a coordinate operation performing the conversion or transformation work. The search for coordinate operations may use a third argument, optional but recommended, which is the geographic area of the data to transform. That later argument is recommended because coordinate operations are often valid only in a some geographic area (typically a particular country or state), and many transformations may exist for the same pair of source and target CRS but different domain of validity. Different coordinate operations may also be different compromises between accuracy and their domain of validity, and specifying a smaller area of interest may allow Apache SIS to select a more accurate operation.

Getting a coordinate operation

The easiest way to obtain a coordinate operation from above-cited information is to use the org.apache.sis.referencing.CRS convenience class:

CoordinateOperation cop = CRS.findOperation(sourceCRS, targetCRS, areaOfInterest);

Among the information provided by CoordinateOperation object, the following are of special interest:

• The domain of validity, either as a textual description (e.g. “Canada – onshore and offshore”) or with the coordinates of a geographic bounding box.
• The positional accuracy, which may be anything from 1 centimetre to a few kilometres.
• The coordinate operation subtype. Among them, two sub-types provide the same functionalities but with a significant conceptual difference:
  • Coordinate conversions are fully determined by mathematical formulas. Those conversions would have an infinite precision if it was not for the unavoidable rounding errors inherent to floating-point calculations. Map projections are in this category.
  • Coordinate transformations are defined empirically. They often have errors of a few metres which are not caused by limitation in computer accuracy. Those errors exist because transformations are only approximations of a more complex reality. Datum shifts from NAD27 to NAD83 are in this category.

If the coordinate operation is an instance of Transformation, then the instance selected by SIS may be one among many possibilities depending on the area of interest. Furthermore its accuracy is usually less than the centimetric accuracy that we can expect from a Conversion. Consequently verifying the domain of validity and the positional accuracy declared in the transformation metadata is of particular importance.

Executing an operation on coordinate values

The CoordinateOperation object introduced in above section provides high-level informations (source and target CRS, domain of validity, positional accuracy, operation parameters, etc). The actual mathematical work is performed by a separated object obtained by a call to CoordinateOperation.getMathTransform(). At the difference of CoordinateOperation instances, MathTransform instances do not carry any metadata. They are kind of black box which know nothing about the source and target CRS (actually the same MathTransform can be used for different pairs of CRS if the mathematical work is the same), domain or accuracy. Furthermore MathTransform may be implemented in a very different way than what CoordinateOperation said. In particular many conceptually different coordinate operations (e.g. longitude rotations, change of units of measurement, conversions between two Mercator projections on the same datum, etc.) are implemented by MathTransform as affine transforms and concatenated for efficiency, even if CoordinateOperation reports them as a chain of Mercator and other operations. The “conceptual versus real chain of coordinate operations” section explains the differences in more details.

The following Java code performs a map projection from geographic coordinates on the World Geodetic System 1984 (WGS84) datum coordinates in the WGS 84 / UTM zone 33N coordinate reference system. In order to make the example a little bit simpler, this code uses predefined constants given by the CommonCRS convenience class. But more advanced applications will typically use EPSG codes instead. Note that all geographic coordinates below express latitude before longitude.

import org.opengis.geometry.DirectPosition;
import org.opengis.referencing.crs.CoordinateReferenceSystem;
import org.opengis.referencing.operation.CoordinateOperation;
import org.opengis.referencing.operation.TransformException;
import org.opengis.util.FactoryException;
import org.apache.sis.referencing.CRS;
import org.apache.sis.referencing.CommonCRS;
import org.apache.sis.geometry.DirectPosition2D;

public class MyApp {
    public static void main(String[] args) throws FactoryException, TransformException {
        CoordinateReferenceSystem sourceCRS = CommonCRS.WGS84.geographic();
        CoordinateReferenceSystem targetCRS = CommonCRS.WGS84.universal(40, 14);  // Get whatever zone is valid for 14°E.
        CoordinateOperation operation = CRS.findOperation(sourceCRS, targetCRS, null);

        // The above lines are costly and should be performed only once before to project many points.
        // In this example, the operation that we got is valid for coordinates in geographic area from
        // 12°E to 18°E (UTM zone 33) and 0°N to 84°N.

        DirectPosition ptSrc = new DirectPosition2D(40, 14);           // 40°N 14°E
        DirectPosition ptDst = operation.getMathTransform().transform(ptSrc, null);

        System.out.println("Source: " + ptSrc);
        System.out.println("Target: " + ptDst);
    }
}

Partial derivatives of coordinate operations

Previous section shows how to project a coordinate from one reference system to another one. There is another, less known, operation which does not compute the projected coordinates of a given point, but instead the derivative of the projection function at that point. Let P be a map projection converting degrees of latitude and longitude (φ, λ) into projected coordinates (x, y) in metres. The formula below represents the map projection result as a column matrix (reason will become clearer soon):

DirectPosition geographic = new DirectPosition2D(φ, λ);
DirectPosition projected = P.transform(geographic, null);
double x = projected.getOrdinate(0);
double y = projected.getOrdinate(1);

The map projection partial derivate at this point can be represented by a Jacobian matrix:

DirectPosition geographic = new DirectPosition2D(φ, λ);
Matrix jacobian = P.derivative(geographic);
double dx_dλ = jacobian.getElement(0,1);
double dy_dφ = jacobian.getElement(1,0);

The first matrix column tells us that if we apply a displacement of 1° of latitude from the (φ, λ) position, — in other words if we move at the (φ + 1, λ) geographic position — then the projected coordinates would be displaced by (∂x, ∂λ) metres — in other words they would become (x + ∂x, y + ∂y) — if the map projection is approximated by an affine transform valid at the (φ, λ) position. Similarly the last matrix column gives us the displacement that happen on the projected coordinate if we apply a displacement of 1° of longitude on the source geographic coordinate under the same assumption. We can visualize such displacements in a figure like below. This figure shows the derivative at two points, P₁ and P₂, for emphasing that the result change for every points. In that figure, vectors U et V stand for the first and second column respectively in the Jacobian matrices.

where vectors are related to the matrix by:

$$\begin{array}{ccc}\overrightarrow{U}=\left[\begin{array}{c}\frac{\partial x}{\partial \phi }\\ \frac{\partial y}{\partial \phi }\end{array}\right]& \text{et}& \overrightarrow{V}=\left[\begin{array}{c}\frac{\partial x}{\partial \lambda }\\ \frac{\partial y}{\partial \lambda }\end{array}\right]\end{array}$$

Above figure shows one usage of map projection derivatives: they provide the direction of parallels and meridians at a given location in a map projection. One can use that information for determining if axes have been swapped or their direction reversed. But the usefulness of map projection derivatives goes further. The annex explains how derivatives are used by the Apache SIS implementation of envelope and raster reprojections.