The meaning of words sometime depend on the community using them. The Apache SIS library prefers as much as possible to use terms in the sense of OGC and ISO standards. Particular care must be taken with the interfaces between SIS and certain other external libraries. For example, the ISO 19123 standard see CV_Coverage as functions in which the domain is the set of spatio-temporal coordinates encompassed by the data, and the range is the set of values encompassed. However, UCAR’s netCDF library applies these terms instead to the function for converting pixel indices (its domain) to spatial-temporal coordinates (its range). Thus the UCAR library’s range may be the domain of ISO 19123.

Coverage

Feature that acts as a function to return values from its range for any direct position within its spatial, temporal or spatiotemporal domain. ISO 19123

Coordinate

One of a sequence of numbers designating the position of a point. ISO 19111

Coordinate operation

Process using a mathematical model, based on a one-to-one relationship, that changes coordinates in a source coordinate reference system to coordinates in a target coordinate reference system, or that changes coordinates at a source coordinate epoch to coordinates at a target coordinate epoch within the same coordinate reference system. ISO 19111

Coordinate reference system

Coordinate system that is related to an object by a datum. ISO 19111

Coordinate system

Set of mathematical rules for specifying how coordinates are to be assigned to points. ISO 19111

Coordinate tuple

Tuple composed of coordinates. The number of coordinates in the coordinate tuple equals the dimension of the coordinate system. The order of coordinates in the coordinate tuple is identical to the order of the axes of the coordinate system. ISO 19111

Datum

Parameter or set of parameters that realize the position of the origin, the scale, and the orientation of a coordinate system. ISO 19111

Domain

Well-defined set. ISO 19123

Range

Set of feature attribute values associated by a function with the elements of the domain of a coverage. ISO 19123

The elements defined in a computer language, such as classes and methods in Java or elements in an XML document, appear in monospaced font in this document. In order to facilitate an understanding of the relationships between Apache SIS and the standards, these elements are also represented using the following color codes:

• Elements in blue are defined in an ISO or OGC standard other than GeoAPI.
• Elements in green are Java element defined in GeoAPI.
• Elements in brown are defined in Apache SIS.
• Elements left in black are either defined elsewhere (for example the standard Java library), or there is no emphasis on that element for the discussion.

For example to represent a projected coordinate reference system (Mercator, Lambert, etc), SC_ProjectedCRS is an UML and XML element defined by the ISO 19111 standard. Then org.​opengis.​referencing.​crs.​ProjectedCRS is the implementation-neutral GeoAPI interface derived from that standard, and org.​apache.​sis.​referencing.​crs.​DefaultProjectedCRS is the implementation class provided by Apache SIS.

Apache SIS implements most GeoAPI interfaces by a classes of the same name than the interface but prefixed by Abstract, Default or General word. The General prefix is sometimes used instead of Default to indicate that alternative implementations are available for some specific cases. For example the Envelope interface is implemented by at least two Apache SIS classes: GeneralEnvelope and Envelope2D. The first implementation can represent envelopes with any number of dimensions while the second implementation is specialized for two-dimensional envelopes.

Apache SIS classes prefixed by Abstract should not (in principle) be instantiated. Users should instantiate a non-abstract subclass instead. However many SIS classes are only conceptually abstract, without abstract Java keyword in class definition. Such classes can be instantiated by a new AbstractXXX(…) statement despite being conceptually abstract. Such instantiations should be avoided, but are nevertheless permitted in last resort when it is not possible to determine the exact subtype.

The easiest way to use Apache SIS is to declare Maven dependencies in the application project. SIS is divided in about 20 modules, which allow applications to import a subset of the library. The Apache SIS downloads page lists the main modules. The pom.xml fragment below gives all dependencies needed by the code snippets in this document (ignoring core modules such as sis-referencing which are inherited by transitive dependencies). Note that the sis-epsg optional module is not under Apache license. Inclusion of that module is subject to acceptation of EPSG terms of use. It is optional but recommended; see How to use EPSG geodetic dataset page for more information.

<properties>
  <sis.version>1.4</sis.version>
</properties>
<dependencies>
  <dependency>
    <groupId>org.apache.sis.storage</groupId>
    <artifactId>sis-geotiff</artifactId>
    <version>${sis.version}</version>
  </dependency>
  <dependency>
    <groupId>org.apache.sis.storage</groupId>
    <artifactId>sis-netcdf</artifactId>
    <version>${sis.version}</version>
  </dependency>

  <!-- Specialization of GeoTIFF reader for Landsat data. -->
  <dependency>
    <groupId>org.apache.sis.storage</groupId>
    <artifactId>sis-earth-observation</artifactId>
    <version>${sis.version}</version>
  </dependency>

  <!-- The following dependency can be omitted if XML support is not desired. -->
  <dependency>
    <groupId>org.glassfish.jaxb</groupId>
    <artifactId>jaxb-runtime</artifactId>
    <version>4.0.4</version>
    <scope>runtime</scope>
  </dependency>

  <!-- This optional dependency requires agreement with EPSG terms of use. -->
  <dependency>
    <groupId>org.apache.sis.non-free</groupId>
    <artifactId>sis-epsg</artifactId>
    <version>${sis.version}</version>
    <scope>runtime</scope>
  </dependency>
</dependencies>

The sis-epsg optional module needs a directory where it will install the geodetic database. That directory can be anywhere on the local machine, it shall exist (but should be initially empty), and its location should be specified by the SIS_DATA environment variable. For example on a Unix system (replace user by the actual user name and some_directory by anything):

export SIS_DATA=/home/user/some_directory
mkdir $SIS_DATA

It is possible to avoid the need to setup SIS_DATA directory if the sis-epsg dependency is replaced by sis-embedded-data. However the latter is slower, and an SIS_DATA directory is still needed for other purposes such as the installation of datum shift grids.

It is possible to instantiate data structures programmatically in memory. But more often, data are read from files or other kinds of data stores. There is different ways to access those data, but an easy way is to use the DataStores​.open(Object) convenience method. The method argument can be a path to a data file (File, Path, URL, URI), a stream (Channel, DataInput, InputStream, Reader), a connection to a data base (DataSource, Connection) or other kinds of object specific to the data source. The DataStores​.open(Object) method detects data formats and returns a DataStore instance for that format.

DataStore functionalities depend on the kind of data (coverage, feature set, time series, etc.). But in all cases, there is always some metadata that can be obtained. Metadata allows to identify the phenomenon or features described by the data (temperature, land occupation, etc.), the geographic area or temporal period covered by the data, together with their resolution. Some rich data source provides also a data quality estimation, contact information for the responsible person or organization, legal or technical constraints on data usage, the history of processing apply on the data, expected updates schedule, etc. The metadata structures depends on the data formats, but Apache SIS translates all of them in a unique metadata model in order to hide this heterogeneity. This pivot model approach is often used by various libraries, with Dublin Core as a popular choice. For Apache SIS, the chosen pivot model is the ISO 19115 international standard. This model organizes metadata in a tree structure. For example if a data format can provides a geographic bounding box encompassing all data, then that information will always be accessible (regardless the data format) from the root Metadata object under the identification­Info node, then the extent sub-node, and finally the geographic­Element sub-node. For example, the following code read a metadata file from a Landsat-8 image and prints the declared geographic bounding box:

import org.opengis.metadata.Metadata;
import org.opengis.metadata.extent.GeographicBoundingBox;
import org.apache.sis.storage.DataStore;
import org.apache.sis.storage.DataStores;
import org.apache.sis.storage.DataStoreException;
import org.apache.sis.metadata.iso.extent.Extents;

void main() throws DataStoreException {
    try (DataStore store = DataStores.open(new File("LC81230522014071LGN00_MTL.txt"))) {
        Metadata overview = store.getMetadata();

        // Convenience method for fetching value at the "metadata/identification­Info/geographic­Element" path.
        GeographicBoundingBox bbox = Extents.getGeographicBoundingBox(overview);

        System.out.println("The geographic bounding box is:");
        System.out.println(bbox);
    }
}

Above example produces the following output (this area is located in Vietnam):

The geographic bounding box is:
Geographic Bounding Box
  ├─West bound longitude…………………………… 108°20′10.464″E
  ├─East bound longitude…………………………… 110°26′39.66″E
  ├─South bound latitude…………………………… 10°29′59.604″N
  └─North bound latitude…………………………… 12°37′25.716″N

Metadata are covered in more details in a latter chapter. Among metadata elements, there is one which will be the topic of a dedicated chapter: reference­System­Info. Its content is essential for accurate data positioning; without this element, even positions given by latitudes and longitudes are ambiguous. Reference systems have many characteristics that make them apart from other metadata: they are immutable, have a particular Well-Known Text representation and are associated to an engine performing coordinate transformation from one reference system to another.

The domain of a coverage is the set of valid input values. In Apache SIS, the domain of grid coverages is described by the Grid­Geometry class. This class contains the following information:

• A grid extent (a.k.a. grid envelope), often inferred from the image size in pixels.
• A grid to CRS conversion, typically as a scale followed by a translation.
• A georeferenced envelope, which can be inferred from the grid extent and the grid to CRS conversion.
• A Coordinate Reference System (CRS) which is the target of the grid to CRS conversion.
• An estimation of grid resolution along each CRS axes.
• An indication of whether conversion for some axes is linear or not.

One of the most important property listed above is the grid to CRS conversion, which defines how to map pixel coordinates to "real world" coordinates such as latitudes and longitudes. This relationship is often linear (an affine transform), but not necessarily; GridGeometry accepts non-linear conversions as well.

The range of a coverage is the set of valid output values. In Apache SIS, the distinction between ranges of numerical values and range of any types of values is represented by Number­Range and Range classes respectively. The Number­Range is used more often, and is also the one that most closely approaches the the common mathematical concept of an interval. This textual representation approaches the specifications of ISO 31-11 standard, except that the comma is replaced by the character “…” as the separator of minimal and maximal values. For example, “[0 … 256)” represents the range of values from 0 inclusive to 256 exclusive.

Range objects are only indirectly associated with coverages. In SIS, the values that can return coverages are described by objects of the Sample­Dimension type. It is these that contain instances of Range, as well as other information such as transfer function (described later).

The SampleDimension.Builder provides convenience methods for building the sample dimensions of a coverage. The usage pattern is to invoke the following methods:

• setName(…) for giving a name to a band.
• addQuantitative(…) for declaring a range of sample values to convert to units of measurement.
• addQualitative(…) for declaring "no data" values.
• setBackground(…) for declaring a "no data" value which can also be used for filling empty space.

Envelopes store minimal and maximal coordinate values of a geometry. Envelopes are not geometries themselves; they are not infinite sets of points (TransfiniteSet). There is no guarantee that all the positions contained within the limits of an envelope are geographically valid. Envelopes must be seen as information about extreme values that might take the coordinates of a geometry as if each dimension were independent of the others, nothing more. Nevertheless, we speak of envelopes as rectangles, cubes or hyper-cubes (depending on the number of dimensions) in order to facilitate discussion, while bearing in mind their non-geometric nature.

An envelope might be represented by two positions corresponding to two opposite corners of a rectangle, cube or hyper-cube. For the first corner, we often take the one whose ordinates all have the maximal value (upperCorner). When displayed using a conventional system of coordinates (with y axis values running upwards), these two positions appear respectively in the lower left corner and the upper right corner of a rectangle. Care must be taken with different coordinate systems, however, which may vary the positions of these corners on the screen. The expressions lower corner and upper corner should thus be understood in the mathematical rather than the visual sense.

3.1.1. Crossing the antimeridian

Minimums and maximums are the values most often assigned to lowerCorner and upperCorner. But the situation becomes complicated when an envelope crosses the antimeridian (−180° or 180° longitude). For example, an envelope 10° in size may begin at 175° longitude and end at −175°. In this case, the longitude value assigned to lowerCorner is greater than that assigned to upperCorner. Apache SIS therefore uses a slightly different definition of these two corners:

• lowerCorner: the starting point, if we move along the inside of the envelope in the direction of ascending values.
• upperCorner: the end-point, if we move along the inside of the envelope in the direction of ascending values.

If the envelope does not cross the antimeridian, these two definitions are equivalent to the selection of minimal and maximal values respectively. This is the case in the green rectangle in the figure below. When the envelope crosses the antimeridian, the lowerCorner and the upperCorner appear again at the bottom and top of the rectangle (assuming a standard system of coordinates), so their names remain appropriate from a visual standpoint. However, the left and right positions are switched. This case is illustrated by the red rectangle in the figure below.

The notions of inclusion and intersection, however, are interpreted slightly differently in these two cases. In the usual case where the envelope does not cross the antimeridian, the green rectangle covers a region of inclusion. The regions excluded from this rectangle continue on to infinity in all directions. In other words, the region of inclusion is not repeated every 360°. But in the case of the red rectangle, the information provided by the envelope actually covers a region of exclusion between the two edges of the rectangle. The region of inclusion extends to infinity to the left and right. We could stipulate that all longitudes below −180° or above 180° are considered excluded, but this would be an arbitrary decision that would not be an exact counterpart to the usual case (green rectangle). A developer may wish to use these values, for example, in a mosaic where the map of the world is repeated several times horizontally and each repetition is considered distinct. If developers wish to perform operations as though the regions of inclusion or exclusion were repeated every 360°, they themselves will have to bring the longitudinal values between −180° and 180° in advance. All the add(…), contains(…), intersect(…), etc. functions of all the envelopes defined in the org.apache.sis.geometry package perform their calculations according to this convention.

In order for functions such as add(…) to work correctly, all objects involved must use the same coordinate reference system, including the same range of values. Thus an envelope that expresses longitudes in the range [−180 … +180]° is not compatible with an envelope that expresses longitudes in the range [0 … 360]°. The conversions, if necessary, are up to the user (the Envelopes class provides convenience methods to do this). Moreover, the envelope’s coordinates must be included within the system of coordinates, unless the developer explicitly decides to consider (for example) 300° longitude as a position distinct from −60°. The GeneralEnvelope class provides a normalize() method to bring coordinates within the desired limits, sometimes at the cost of lower values being higher than upper values.

3.1.2. Transforming to another reference system

Geographic information systems often need to transform an envelope from one Coordinate Reference System (CRS) to another. But a naive approach transforming the 4 corners is not sufficient. The figure below shows an envelope before a map projection and the geometric shape that we would get if all points (not only the corners) were projected. The resulting geometric shape is more complex than a rectangle because of the curvature caused by the map projection. Computing the envelope that contains the 4 corners of that shape is not enough, because the area in the bottom of the projected shape is lower than the two bottom corners. That surface would be outside the envelope.

Envelope before projection

Geometric shape after projection

Sampling a larger number of points reduces the problem but does not resolve it. Map projection derivatives offer a more efficient way to resolve this problem (see the annex for more mathematical details). Another complication occurs if the envelope contains the North or South pole. For making a long story short, transforming an envelope is a lot more complicated than it looks like. Apache SIS contains a few utility methods for making this task easier. For transforming an envelope to another CRS (WGS 84 / World Mercator in this example):

CoordinateReferenceSystem targetCRS = CRS.forCode("EPSG:3395");
Envelope transformed = Envelopes.transform(envelope, targetCRS);

If envelopes are transformed in the goal of using a common CRS before to compute the union of many envelopes, an additional complication is that each envelope may use a CRS with a relatively small domain of validity. The union operation needs to find a CRS valid in a domain large enough for containing all envelopes. It may be a CRS different than all CRS used by the source envelopes. Apache SIS has an utility method for handling this additional complexity. This method accepts an arbitrary number of envelopes that may be in different CRS:

Envelope union = Envelopes.union(envelope1, envelope2, envelope3);

Spatial reference systems by coordinates provide necessary information for mapping numerical coordinate values to real-world locations. In Apache SIS, most information is contained (directly or indirectly) in classes with a name ending in CRS, the abbreviation of Coordinate Reference System. Those objects contain:

• A datum, which specifies among other things which ellipsoid to use as an Earth shape approximation.
• A description for each axis: name, direction, units of measurement, range of values.
• Sometimes a list of parameters, especially when using map projections.

Those systems are described by the ISO 19111 standard (Referencing by Coordinates), which replaces for most parts the older OGC 01-009 standard (Coordinate Transformation Services). Those standards are completed by two other standards defining exchange formats: ISO 19136 and 19162 respectively for the Geographic Markup Language (GML) — a XML format which is quite detailed but verbose — and the Well-Known Text (WKT) — a text format easier to read by humans.

4.1.1. Map projections

Map projections represent a curved surface (the Earth surface) on a plane surface (a map or a computer screen). Every rendering of geospatial data on a flat screen uses some kind of map projection, sometimes implicitly. Well-designed map projections provide some control over deformations: one can preserve the angles, another projection can preserve the areas, but none can preserve both in same time. The geometric properties to preserve depend on the feature to represent and the work to do on that feature. For example countries elongated along the East-West axis often use a Lambert projection, while countries elongated along the North-South axis prefer a Transverse Mercator projection.

There is thousands of projected CRS in use around the world. Many of them are published in the EPSG geodetic database. The easiest way to get a projected CRS with Apache SIS is to use its EPSG code. For example the following code gets the definition of the JGD2000 / UTM zone 54N CRS (for Japan from 138°E to 144°E):

CoordinateReferenceSystem crs = CRS.forCode("EPSG:3100");

Other ways to get a coordinate reference system will be given in a next section.

4.1.2. Geographic reference systems

All map projections are based on a geodetic (usually geographic) CRS. A geodetic CRS is a coordinate reference system with latitude, longitude and sometimes height axes. There is many kinds of latitudes and longitudes, but two common kinds supported by Apache SIS are geodetic and geocentric latitudes. Those two angles differ slightly in the way they intersect the ellipsoid surface. On Earth surface, the difference between those two kinds of latitude varies between 0 and about 20 km.

When peoples talk about latitudes and longitudes, they usually mean geodetic latitudes and longitudes. A coordinate reference system using such latitudes and longitudes is said geographic and is represented by the GeographicCRS interface. Systems using the other kinds of latitude are represented by other CRS interfaces.

Theoretically, data expressed in a geographic CRS can never be rendered directly on a flat screen (they could be rendered directly on a planetarium dome however). In practice we allow data rendering in a geographic CRS, but this process implicitly uses a Plate Carrée projection.

4.1.3. Vertical and temporal dimensions

TODO

4.1.4. Coordinate systems

A Coordinate System (CS) defines the set of axes that spans a given coordinate space. Each axis defines an approximative direction (north, south, east, west, up, down, port, starboard, past, future, etc.), units of measurement, minimal and maximal values, and what happen after reaching those extremum. For example in longitude case, after +180° the coordinate values continue at −180°. Axes having such behavior are flagged by the RangeMeaning.WRAPAROUND code.

Each Coordinate Reference System (CRS) is associated with exactly one Coordinate System (CS). Some properties that we can get from a coordinate system and its axes are shown below. Axes are numbered from 0 to cs.getDimension()-1 inclusive.

CoordinateSystem cs = crs.getCoordinateSystem();
CoordinateSystemAxis secondAxis = cs.getAxis(1);            // For a geographic CRS, this is usually geodetic longitude.
String        abbreviation = secondAxis.getAbbreviation();  // For a longitude axis, this is usually "λ", "L" or "lon".
AxisDirection direction    = secondAxis.getDirection();     // For a longitude axis, this is usually EAST. Another occasional value is WEST.
Unit<?>       units        = secondAxis.getUnit();          // For a longitude axis, this is usually Units.DEGREE.
double        minimum      = secondAxis.getMinimumValue();  // For a longitude axis, this is usually −180°. Another common value is 0°.
double        maximum      = secondAxis.getMaximumValue();  // For a longitude axis, this is usually +180°. Another common value is 360°.
RangeMeaning  atEnds       = secondAxis.getRangeMeaning();  // For a longitude axis, this is WRAPAROUND.

In addition to axis definitions, another important coordinate system characteristic is their type (CartesianCS, SphericalCS, etc.). The CS type implies the set of mathematical rules for calculating geometric quantities like angles, distances and surfaces. Usually the various CS subtypes do not define any new Java methods compared to the parent type, but are nevertheless important for type safety. For example many calculations or associations are legal only when all axes are perpendicular to each other. In such case the coordinate system type is restricted to CartesianCS in method signatures.

Coordinate systems are mathematical concepts; they do not contain any information about where on Earth is located the system origin. Consequently coordinate systems alone are not sufficient for describing a location; they must be combined with a datum (or reference frame). Those combinations form the coordinate reference systems described in previous sections.

4.1.5. Geodetic datum

Since the real topographic surface is difficult to represent mathematically, it is not used directly. A slightly more convenient surface is the geoid, a surface where the gravitational field has the same value everywhere (an equipotential surface). This surface is perpendicular to the direction of a plumb line at all points. The geoid surface would be equivalent to the mean sea level if all oceans where at rest, without winds or permanent currents like the Gulf Stream.

While much smoother than topographic surface, the geoid surface still have hollows and bumps caused by the uneven distribution of mass inside Earth. For more convenient mathematical operations, the geoid surface is approximated by an ellipsoid. This “figure of Earth” is represented in GeoAPI by the Ellipsoid interface, which is a fundamental component in coordinate reference systems of type GeographicCRS and ProjectedCRS. Tenth of ellipsoids are commonly used for datum definitions. Some of them provide a very good approximation for a particular geographic area at the expense of the rest of the world for which the datum was not designed. Other datums are compromises applicable to the whole world.

The main properties that we can get from an ellipsoid are given below. The semi-major axis length is sometimes called equatorial radius and the semi-minor axis length the polar radius. The inverse flattening factor is apparently superfluous since it can be derived from other quantities, but many ellipsoid definitions provide this factor instead of semi-minor axis length.

Unit<Length> units = ellipsoid.getAxisUnit();
double semiMajor   = ellipsoid.getSemiMajorAxis();          // In units of measurement given above.
double semiMinor   = ellipsoid.getSemiMinorAxis();          // In units of measurement given above.
double ivf         = ellipsoid.getInverseFlattening();      // = semiMajor / (semiMajor - semiMinor).

For defining a geodetic system in a country, a national authority selects an ellipsoid matching closely the country surface. Differences between that ellipsoid and the geoid’s hollows and bumps are usually less than 100 metres. Parameters that relate an Ellipsoid to the Earth surface (for example the position of ellipsoid center) are represented by instances of GeodeticDatum. Many GeodeticDatum definitions can use the same Ellipsoid, but with different orientations or center positions.

Before the satellite era, geodetic measurements were performed exclusively from Earth surface. Consequently, two islands or continents not in range of sight from each other were not geodetically related. So the North American Datum 1983 (NAD83) and the European Datum 1950 (ED50) are independent: their ellipsoids have different sizes and are centered at a different positions. The same geographic coordinate will map different locations on Earth depending on whether the coordinate uses one reference system or the other.

The GPS invention implied the creation of a world geodetic system named WGS84. The ellipsoid is then unique and centered at the Earth gravity center. GPS provides at any moment the receptor absolute position on that world geodetic system. But since WGS84 is a world-wide system, it may differs significantly from local systems. For example the difference between WGS84 and the European system ED50 is about 150 metres, and the average difference between WGS84 and the Réunion 1947 system is 1.5 kilometres. Consequently we shall not blindly use GPS coordinates on a map, as transformations to the local system may be required. Those transformations are represented in GeoAPI by instances of the Transformation interface.

The WGS84 ubiquity tends to reduce the need for Transformation operations with recent data, but does not eliminate it. The Earth moves under the effect of plate tectonic and new systems are defined every years for taking that fact in account. For example while NAD83 was originally defined as practically equivalent to WGS84, there is now (as of 2016) a 1.5 metres difference. The Japanese Geodetic Datum 2000 was also defined as practically equivalent to WGS84, but the Japanese Geodetic Datum 2011 now differs. Even the WGS84 datum, which was a terrestrial model realization at a specific time, got revisions because of improvements in instruments accuracy. Today, at least six WGS84 versions exist. Furthermore many borders were legally defined in legacy datums, for example NAD27 in USA. Updating data to the new datum would imply transforming some straight lines or simple geometric shapes into more irregular shapes, if the shapes are large enough.

Contrarily to other kinds of objects introduced in this section, there is not many useful information that we can get from a Datum instance except its name. It is difficult to translate in programming language how a datum is related to the Earth. Often, the most we can do is to consider that having two datums with different names implies that the same location on Earth has different coordinate values when using those different datums, even if the ellipsoid is identical in both cases. Coordinate transformations between datums require some kind of database.

TODO:

• Using CommonCRS
• Looking CRS defined by authorities with CRSAuthorityFactory
• Reading definitions in GML or WKT format
• Constructing programmatically using CRSFactory

4.2.1. Adding new CRS definitions

TODO

The axis order is specified by the authority (typically a national agency) defining the Coordinate Reference System (CRS). The order depends on the CRS type and the country defining the CRS. In the case of geographic CRS, the (latitude, longitude) axis order is widely used by geographers and pilots for centuries. However software developers tend to consistently use the (x, y) order for every kind of CRS. Those different practices resulted in contradictory definitions of axis order for almost every CRS of kind GeographicCRS, for some ProjectedCRS in the South hemisphere (South Africa, Australia, etc.) and for some polar projections among others.

Recent OGC standards mandate the use of axis order as defined by the authority. Oldest OGC standards used the (x, y) axis order instead, ignoring any authority specification. Many software products still use the old (x, y) axis order, maybe because such uniformization makes CRS implementation and usage apparently easier. Apache SIS supports both conventions with the following approach: by default, SIS creates CRS with axis order as defined by the authority. Those CRS are created by calls to the CRS.forCode(String) method and the actual axis order can be verified after the CRS creation with System.out​.println(crs). But if (x, y) axis order is wanted for compatibility with older OGC specifications or other software products, then CRS forced to longitude first axis order can be created by a call to the following method:

CoordinateReferenceSystem crs = …;               // CRS obtained by any means.
crs = AbstractCRS.castOrCopy(crs).forConvention(AxesConvention.RIGHT_HANDED)

Among the legacy OGC standards that used the non-conform axis order, an influent one is version 1 of the Well Known Text (WKT) format specification. According that widely-used format, WKT 1 definitions without explicit AXIS[…] elements shall default to (longitude, latitude) or (x, y) axis order. In version 2 of WKT format, AXIS[…] elements are no longer optional and should contain an explicit ORDER[…] sub-element for making the intended order yet more obvious. But if AXIS[…] elements are nevertheless missing in a WKT 2 definition, Apache SIS defaults to (latitude, longitude) order. So in summary:

• Default WKT 1 axis order of geographic CRS is (longitude, latitude) as mandated by OGC 01-009 specification.
• Default WKT 2 axis order of geographic CRS is (latitude, longitude), but this is SIS-specific as ISO 19162 does not mention default axis order.

To avoid ambiguities, users are encouraged to always provide explicit AXIS[…] elements in their WKT.

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.

4.4.1. 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.

4.4.2. 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);
    }
}

4.4.3. 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.

4.4.4. Chain of coordinate operation steps

Coordinate operations may include many steps, each with their own set of parameters. For example transformations from one datum (e.g. NAD27) to another datum (e.g. WGS84) can be approximated by an affine transform (translation, rotation and scale) applied on the geocentric coordinates. This implies that the coordinates must be converted from geographic to geocentric domain before the affine transform, then back to geographic domain after the affine transform. The result is a three-steps process illustrated in the “Conceptual chain of operations” column of the example below. However because that operation chain is very common, the EPSG geodetic dataset provides a shortcut named “Geocentric translation in geographic domain”. Using this operation, the conversion steps between geographic and geocentric CRS are implicit. Consequently the datum shifts as specified by EPSG appears as if it was a single operation, but the real operation executed by Apache SIS is divided in more steps.

All those operation chains can be viewed in Well Known Text (WKT) or pseudo-WKT format. The simplest operation chain, as specified by the authority, is given directly by the String representation of the CoordinateOperation instance. This WKT 2 representation contains not only a description of operations with their parameter values, but also additional information about the context in which the operation applies (the source and target CRS) together with some metadata like the accuracy and domain of validity. Some operation steps and parameters may be omitted if they can be inferred from the context.

CoordinateOperation["NAD27 to WGS 84 (3)",
  SourceCRS[full CRS definition required here but omitted for brevity],
  TargetCRS[full CRS definition required here but omitted for brevity],
  Method["Geocentric translations (geog2D domain)"],
    Parameter["X-axis translation", -10.0, Unit["metre", 1]],
    Parameter["Y-axis translation", 158.0, Unit["metre", 1]],
    Parameter["Z-axis translation", 187.0, Unit["metre", 1]],
  OperationAccuracy[20.0],
  Area["Canada - onshore and offshore"],
  BBox[40.04, -141.01, 86.46, -47.74],
  Id["EPSG", 1172, "8.9"]]

An operation chain closer to what Apache SIS really performs is given by the String representation of the MathTransform instance. In this WKT 1 representation, contextual information and metadata are lost; a MathTransform is like a mathematical function with no knowledge about the meaning of the coordinates on which it operates. Since contextual information are lost, implicit operations and parameters become explicit. This representation is useful for debugging since any axis swapping operation (for example) become visible. Apache SIS constructs this representation from the data structure in memory, but convert them in a more convenient form for human, for example by converting radians to degrees.

Concat_MT[
  Param_MT["Affine parametric transformation",
    Parameter[parameters performing axis swapping omitted for brevity]],
  Inverse_MT[Param_MT["Geographic3D to 2D conversion"]],
  Param_MT["Geographic/geocentric conversions",
    Parameter["semi_major", 6378206.4],
    Parameter["semi_minor", 6356583.8]],
  Param_MT["Geocentric translations (geocentric domain)",
    Parameter["X-axis translation", -10.0],
    Parameter["Y-axis translation", 158.0],
    Parameter["Z-axis translation", 187.0]],
  Param_MT["Geocentric_To_Ellipsoid",
    Parameter["semi_major", 6378137.0],
    Parameter["semi_minor", 6356752.314245179]],
  Param_MT["Geographic3D to 2D conversion"],
  Param_MT["Affine parametric transformation",
    Parameter[parameters performing axis swapping omitted for brevity]]]

The latter form is often useful for debugging. If a coordinate operation seems to produce wrong results, inspecting the Well Known Text like above should be the first thing to do. The Frequently Asked Questions page gives more tips about common causes of coordinate transformation errors.

The metadata modules provide support methods for handling the metadata objects through Java Reflection. This is an approach similar to Java Beans, in that users are encouraged to use directly the API of Plain Old Java objects every time their type is known at compile time, and fallback on the reflection technic when the type is known only at runtime. When using Java reflection, a metadata can be viewed in different ways:

• As key-value pairs in a Map (from java.util).
• As a TreeTable (from org.apache.sis.util.collection).
• As a table record in a database (using org.apache.sis.metadata.sql).
• As an XML document conforms to ISO standard schema.

The use of reflection is described below. The XML representation is described in a separated chapter.

5.1.1. Direct access via getter methods

All metadata classes provide getter, and sometime setter, methods for their properties. Some properties accept many values, in which case the property type is a collection. The following example prints the range of latitudes of all data descriptions found in a given root Metadata object:

import org.opengis.metadata.metadata;
import org.opengis.metadata.extent.Extent;
import org.opengis.metadata.extent.GeographicExtent;
import org.opengis.metadata.extent.GeographicBoundingBox;
import org.opengis.metadata.identification.Identification;
import org.opengis.metadata.identification.DataIdentification;

void main() {
    Metadata metadata = ...;    // For example, metadata read from a data store.
    for (Identification identification : metadata.getIdentificationInfo()) {
        if (identification instanceof DataIdentification data) {
            for (Extent extent : data.getExtents()) {
                // Extents may have horizontal, vertical and temporal components.
                for (GeographicExtent horizontal : extent.getGeographicElements()) {
                    if (horizontal instanceof GeographicBoundingBox bbox) {
                        double south = bbox.getSouthBoundLatitude();
                        double north = bbox.getNorthBoundLatitude();
                        System.out.println("Latitude range: " + south + " to " + north);
                    }
                }
            }
        }
    }
}

Because of ISO 19115 richness, interesting information may be buried deeply in the metadata tree, as in above example. For a few frequently-used elements, some convenience methods are provided. Those conveniences are generally defined as static methods in classes having a name in plural form. For example the Extents class defines static methods for fetching more easily some information from Extent metadata elements. For example the following method navigates through different branches where North, South, East and West data bounds may be found:

import org.opengis.metadata.metadata;
import org.opengis.metadata.extent.GeographicBoundingBox;
import org.apache.sis.metadata.iso.extent.Extents;

void main() {
    Metadata metadata = ...;    // For example, metadata read from a data store.
    GeographicBoundingBox bbox = Extents.getGeographicBoundingBox(extent);
}

Those conveniences are defined as static methods in order to allow their use with different metadata implementations. Some other classes providing static methods for specific interfaces are Citations, Envelopes, Matrices and MathTransforms.

5.1.2. View as key-value pairs

Above static methods explore fragments of metadata tree in search for requested information, but the searches are still targeted to elements whose types and at least part of their paths are known at compile-time. Sometimes the element to search is known only at runtime, or sometimes there is a need to iterate over all elements. In such cases, one can view the metadata as a java.util.Map like below:

import java.util.Map;
import org.apache.sis.metadata.MetadataStandard;
import org.apache.sis.metadata.KeyNamePolicy;
import org.apache.sis.metadata.ValueExistencePolicy;

void main() {
    Map<String,Object> elements = MetadataStandard.ISO_19115.asValueMap(
            metadata,                           // Any instance from the org.opengis.metadata package or a sub-package.
            null,                               // Used for resolving ambiguities. We can ignore for this example.
            KeyNamePolicy.JAVABEANS_PROPERTY,   // Keys in the map will be getter method names without "get" prefix.
            ValueExistencePolicy.NON_EMPTY);    // Entries with null or empty values will be omitted.

    // Print the names of all root metadata elements having a value.
    for (String name : elements.keySet()) {
        System.out.println(name);
    }
}

The Map object returned by asValueMap(…) is live: any change in the metadata instance will be immediately reflected in the view. Actually, each map.get("foo") call is forwarded to the corresponding metadata​.getFoo() method. Conversely, any map.put("foo", …) or map.remove("foo") operation applied on the view will be forwarded to the corresponding metadata​.setFoo(…) method, if that method exists. The view is lenient regarding keys given in arguments to Map methods: keys may be property names ("foo"), method names ("getFoo"), or names used in ISO 19115 standard UML diagrams (similar to property names but not always identical). Differences in upper cases and lower cases are ignored when this tolerance does not introduce ambiguities. For more information on metadata views, see org.apache.sis.metadata package javadoc.

5.1.3. View as tree table

A richer alternative to the view as a map is the view as a tree table. With this view, the metadata structure is visible as a tree, and each tree node is a table row with the following columns:

Tree table views are obtained in a way similar to map views, but using the asTreeTable(Object) method instead of asValueMap(Object).

Apache SIS allows any metadata property to be null (for values that are not collections) or an empty collection. However, ISO 19115 defines some properties as mandatory. For example, the title of a Citation is mandatory, while the edition is optional. Apache SIS will not raise any warning or error if a mandatory property is missing. But for strict ISO compliance, the reason why the property is missing should be provided. Predefined reasons are:

The transmission of this information requires the use of a non-null object, even when the value is missing. In such case, SIS will return an object that, besides implementing the desired GeoAPI interface, also implements the org.apache.sis.xml.NilObject interface. This interface flags the instances where all methods return an empty collection, an empty table, null, NaN, 0 or false, in this preference order, as permitted by the return types of the methods. Each instance that implements NilObject provides a getNilReason() method indicating why the object is nil.

For each metadata class, there is an XML type with the same name than in the abstract specification (for example, mdb:MD_Metadata and cit:CI_Citation). All of these types may be used as the root of an XML document. It is therefore possible to write a document representing a complete MD_Metadata object, or to write a document representing only a CI_Citation object.

ISO 19115-3 standard arranges the content of these objects in an unusual way: for each property whose value type is itself another class of ISO 19115, the value is wrapped in an element that represents its type, rather than being written directly. For example, in an object of the CI_Citation type, the value of the citedResponsibleParty property is incorporated into a CI_Responsibility element. This practice doubles the depth of the hierarchy, and introduces duplication at all levels for each value, as in the following example:

<MD_Metadata>
  <identificationInfo>
    <MD_DataIdentification>
      <citation>
        <CI_Citation>
          <citedResponsibleParty>
            <CI_Responsibility>
              <party>
                <CI_Party>
                  <contactInfo>
                    <CI_Contact>
                      <onlineResource>
                        <CI_OnlineResource>
                          <linkage>
                            <URL>https://www.ogc.org</URL>
                          </linkage>
                        </CI_OnlineResource>
                      </onlineResource>
                    </CI_Contact>
                  </contactInfo>
                </CI_Party>
              </party>
            </CI_Responsibility>
          </citedResponsibleParty>
        </CI_Citation>
      </citation>
    </MD_DataIdentification>
  </identificationInfo>
</MD_Metadata>

The preceding example, like all documents that conform to ISO 19115-3, consists of a systematic alternation of two types of XML elements:

1. It begins with the name of the property, which always begins with a lower-case letter (ignoring prefixes). In Java APIs, each property corresponds to a method in its enclosing class. In the example above, mdb:identificationInfo (a property of MD_Metadata class) corresponds to the Metadata​.getIdentificationInfo() method.

2. The value type is included under each property, unless it has been replaced with a reference (the following sub-section will elaborate on this subject). The value type is an XML element whose name always begins with an upper-case letter, ignoring prefixes. In the example above we had MD_DataIdentification, which corresponds to the DataIdentification Java interface. It is this XML element that contains the child properties.

In order to reduce the complexity of the libraries, GeoAPI and Apache SIS only expose publicly a single unified view of these two types of elements. The public API basically corresponds to the second group.

6.1.1. Links to previously-defined instances

The parent element may contain an id or uuid attribute. If one of these attributes is present, the parent attribute may be completely omitted; it will be replaced at the time of reading by the element referenced by the attribute. In the following example, the part on the left defines an element associated with the identifier “my_id,” while the part on the right references this element:

<MD_MetaData>
  <identificationInfo>
    <MD_DataIdentification id="my_id">
      <!-- insert child properties here -->
    </MD_DataIdentification>
  </identificationInfo>
</MD_MetaData>

<MD_MetaData>
  <identificationInfo xlink:href="#my_id"/>
</MD_MetaData>

The decision of which attribute to use depends on the scope of the referenced item:

• id is only valid in the same XML document that defines the object it references.
• uuid may be valid outside the XML document, but someone must maintain a database providing the objects for each given UUID.
• xlink:href may reference another XML document accessible on the Internet.

In the SIS library, all objects that can be identified in an XML document implements the org.apache.sis.xml.IdentifiedObject interface. Each instance of this interface provides a view of its identifiers in the form of a Map<Citation,String>, in which the Citation key indicates the type of identifier and the value is the identifier itself. Some constants representing different types of identifiers are listed in IdentifierSpace, including ID, UUID and HREF. Each of these keys may be associated with a different type of value (usually String, UUID or URI) depending on the key. For example, the following code defines a value for the uuid attribute:

import org.apache.sis.metadata.iso.DefaultMetadata;
import org.apache.sis.xml.IdentifierSpace;
import java.util.UUID;

public class MyClass {
    public void myMethod() {
        UUID identifier = UUID.randomUUID();
        DefaultMetadata metadata = new DefaultMetadata();
        metadata.getIdentifierMap().putSpecialized(IdentifierSpace.UUID, identifier);
    }
}

Although this mechanism has been defined in order to better support the representation of XML attributes of the gco:ObjectIdentification group, it also conveniently allows other types of identifiers to be manipulated. For example, the ISBN and ISSN attributes of Citation may be manipulated in this way. The methods of the IdentifiedObject interface therefore provides a specific location where all types of identifiers (not only XML) associated with an object may be manipulated.

6.1.2. Placeholders for missing values

If a mandatory property value is missing, a nil reason should be provided. The reason is specified in a nilReason attribute, which can take values such as unknown, missing, inapplicable, withheld, template, etc. In the following example, the left side shows a CI_Citation element containing a CI_Series element, while on the right side the series is unknown. If the CI_Series element had been completely omitted, then the Citation​.getSeries() method would return null in Java. But when a nilReason is present, the SIS implementation of getSeries() returns instead an object that implements both the Series and NilReason interfaces, and in which the getNilReason() method returns the constant UNKNOWN.

<CI_Citation>
  <series>
    <CI_Series>
      <!-- insert child properties here -->
    </CI_Series>
  </series>
</CI_Citation>

<CI_Citation>
  <series nilReason="unknown"/>
</CI_Citation>

There are various opinions on how to implement Java standard’s Object​.equals(Object) method. According to some, it should be possible to compare different implementations of the same interface or base class. But to follow this policy, each interface or base class’s javadoc must define the algorithms that all implementations shall use for their equals(Object) and hashCode() methods. This approach is common in java.util.Collection and its child interfaces. Transferring this approach to certain GeoAPI interfaces, however, would be a difficult task, and would probably not be followed in many implementations. Moreover, it comes at the expense of being able to take into account supplementary attributes in the child interfaces, unless this possibility has been specified in the parent interface. This constraint arises from the following points of the equals(Object) and hashCode() method contracts:

• A.equals(B) implies B.equals(A) (symmetry);
• A.equals(B) and B.equals(C) implies A.equals(C) (transitivity);
• A.equals(B) implies A.hashCode() == B.hashCode().

For example, these three constraints are violated if A (and eventually C) can contain attributes which B ignores. To bypass this problem, an alternative approach is to require that the objects compared by the Object​.equals(Object) method be of the same class; in other words, A.getClass() == B.getClass(). This approach is sometimes regarded as contrary to the principles of object oriented programming. In practice, for relatively complex applications, the important accorded to these principles depends on the context in which the objects are compared: if the objects are added to a HashSet or used as keys in a HashMap, we would need a stricter adherence to the equals(Object) and hashCode() contract. But if the developer is comparing the objects his- or herself, for example to check that the relevant information has been changed, then the constraints of symmetry, transitivity or coherence with the hash values may be of little interest. More permissive comparisons may be desirable, sometimes going so far as to tolerate minor discrepancies in numerical values.

In order to allow developers a certain amount of flexibility, many classes in the SIS library implement the org.apache.sis.util.LenientComparable interface, which defines a equals(Object, ComparisonMode) method. The principle modes of comparison are:

• STRICT — The objects compared must share the same class and have exactly equal attributes, including any possible public attributes specific to the implementation.

• BY_CONTRACT — The objects compared must implement the same GeoAPI (or other standard) interface, but need not be of the same implementation class. Only the attributes defined in the interface are compared; all other attributes specific to the implementation — even if they are public — are ignored.

• IGNORE_METADATA — Like BY_CONTRACT, but only compares attributes that influence the operations (numeric calculations or otherwise) performed by the object. For example, in a geodesic datum, the longitude (in relation to Greenwich) of the original meridian would be taken into account, while the name of the meridian would be ignored.

• APPROXIMATIVE — Like IGNORE_METADATA, but tolerates minor discrepancies in numerical values.

The default mode, used in all equals(Object) methods in SIS, is STRICT. This mode is chosen for a safe operation — particularly with HashMap — without the need to rigorously define equals(Object) and hashCode() operations in every interface. With this mode, the order of objects (A.equals(B) or B.equals(A)) is unimportant. It is, however, the only mode that offers this guarantee. In the expression A.equals(B), the BY_CONTRACT mode (and so by extension all other modes that depend on it) only compares the properties known to A, regardless of whether B knows more.

There is sometimes a need to convert instances from a source Java type to a target Java type while those types are unknown at compile time. Various projects (Apache Common Convert, Spring, etc.) have created their own interface for performing object conversions between types known only at runtime. Details vary, but such interfaces typically look like below:

interface ObjectConverter<S,T> {   // Some projects use only "Converter" as interface name.
    T apply(S object);             // Another method name commonly found in other projects is "convert".
}

Like other projects, Apache SIS also defines its own ObjectConverter interface. The main difference between SIS converter interface and the interfaces found in other projects is that SIS converters provide some information about their mathematical properties. An Apache SIS converter can have zero, one or many of the following properties:

Injective

A function is injective if no pair of S values can produce the same T value.

Example: the Integer → String conversion performed by Integer​.toString() is an injective function because if two Integer values are not equal, then it is guaranteed that their conversions will result in different String values. However the String → Integer conversion performed by Integer​.valueOf(String) is not an injective function because many distinct String values can be converted to the same Integer value. For example converting the "42", "+42" and "0042" character strings all result in the same 42 integer value.

Surjective

A function is surjective if each values of T can be created from at least one value of S.

Example: the String → Integer conversion performed by Integer​.valueOf(String) is a surjective function because every Integer value can be created from at least one String value. However the Integer → String conversion performed by Integer​.toString() is not a surjective function because it cannot produce all possible String values. For example there is no way to produce the "ABC" value with the Integer​.toString() method.

Bijective

A function is bijective if there is a one-to-one relationship between S and T values.

Note: the bijective property is defined here for clarity, but actually does not have an explicit item in Apache SIS FunctionProperty enumeration. It is not necessary since a function that is both injective and surjective is necessarily bijective.

Order preserving

A function is order preserving if any sequence of increasing S values is mapped to a sequence of increasing T values.

Example: conversion from Integer to Long preserve the natural ordering of elements. However conversions from Integer to String do not preserve natural ordering, because some sequences of increasing integer values are ordered differently when their string representations are sorted by lexicographic order. For example 1, 2, 10 become "1", "10", "2".

Order reversing

A function is order reversing if any sequence of increasing S values is mapped to a sequence of decreasing T values.

Example: a conversion that reverses the sign of numbers.

Above information may seem unnecessary when values are converted without taking in account the context in which the values appear. But when the value to convert is part of a bigger object, then above information can affect the way the converted value will be used. For example conversion of a [min … max] range is straightforward when the converter is order preserving. But if the converter is order reversing, then the minimum and maximum values need to be interchanged. For example if the converter reverses the sign of values, then the converted range is [-max … -min]. If the converter is neither order preserving or order reversing, then range conversion is not allowed at all (because it does not contain the same set of values) even if the minimum and maximum values could be converted individually.

In an architecture where a program executed on a server provides its data to multiple clients, the server’s locale conventions are not necessarily the same as those of the clients. Conventions may differ in language, but also in the way they write numeric values (even between two countries that speak the same language) as well in time zone. To produce messages that conform to the client’s conventions, SIS uses two approaches, distinguished by their level of granularity: at the level of the messages themselves, or at the level of the objects that create the messages. The approach used also determines whether it is possible to share the same instance of an object for all languages.

7.3.1. Distinct character sequences for each locale

Some classes are only designed to function according to one locale convention at a time. This is of course true for the standard implementations of java.text.Format, as they are entirely dedicated to the work of internationalization. But it is also the case for other less obvious classes like javax.imageio.ImageReader and ImageWriter. When one of these classes is implemented by SIS, we identify it by implementing the org.apache.sis.util.Localized interface. The getLocale() method of this interface can determine the locale conventions by which the instance produces its message.

Another class that provides different methods for different locales is java.lang.Throwable. The standard Java API defines two methods for getting the error message: getMessage() and getLocalizedMessage(). Usually those two methods return the same character sequences, but some exceptions thrown by Apache SIS may use different locales. The policy that SIS tries to apply on a best-effort basis is:

• getMessage() returns the message in the JVM default locale. In a client-server architecture, this is often the locale on the server side. This is the recommended language for logging messages to be read by system administrators.
• getLocalizedMessage() returns the message in a locale that depends on the context in which the exception has been thrown. This is often the locale used by a particular Format or DataStore instance, and can be presumed to be the locale on the client side. This is the recommended language to show in the user application.

The utility class org.apache.sis.util.Exceptions provides convenience methods to get messages according to the conventions of a given locale, when this information is available.

7.3.2. Single instance for all supported locales

The API conventions defined by SIS or inherited by GeoAPI favour the use of the InternationalString type when the value of a String type would likely be localized. This approach allows us to defer the internationalization process to the time when a character sequence is requested, rather than the time when the object that contains them is created. This is particularly useful for immutable classes used for creating unique instances independently of locale conventions.

When defining spatial objects independently of locale conventions, we reduce the risk of computational overload. For example, it is easier to detect that two maps use the same cartographic projection if this last is represented by the same instance of CoordinateOperation, even if the projection has a different name depending on the country. Moreover, certain types of CoordinateOperation may require coordinate transformation matrices, so sharing a single instance becomes even more preferable in order to reduce memory consumption.

7.3.3. Locale.ROOT convention

All SIS methods receiving or returning the value of a Locale type accept the Locale.ROOT value. This value is interpreted as specifying not to localize the text. The notion of a non-localized text is a little false, as it is always necessary to chose a formatting convention. This convention however, though very close to English, is usually slightly different. For example:

• Identifiers are written as they appear in UML diagrams, such as blurredImage instead of Blurred image.
• Dates are written according to the ISO 8601 format, which does not correspond to English conventions.
• Numbers are written using their toString() methods, rather than using a java.text.NumberFormat. As a result, there are differences in the number of significant digits, use of exponential notation and the absence of thousands separators.

7.3.4. Treatment of characters

In Java, sequences of characters use UTF-16 encoding. There is a direct correspondence between the values of the char type and the great majority of characters, which facilitates the use of sequences so long as these characters are sufficient. However, certain Unicode characters cannot be represented by a single char. These supplementary characters include certain ideograms, but also road and geographical symbols in the 1F680 to 1F700 range. Support for these supplementary characters requires slightly more complex interactions than the classic case, where we may assume a direct correspondence. Thus, instead of the loop on the left below, international applications must generally use the loop on the right:

for (int i=0; i<string.length(); i++) {
    char c = string.charAt(i);
    if (Character.isWhitespace(c)) {
        // A blank space was found.
    }
}

for (int i=0; i<string.length();) {
    int c = string.codePointAt(i);
    if (Character.isWhitespace(c)) {
        // A blank space was found.
    }
    i += Character.charCount(c);
}

SIS supports supplementary characters by using the loop on the right where necessary, but the loop on the left is occasionally used when it is known that the characters searched for are not supplementary characters, even if some may be present in the sequence in which we are searching.

7.3.4.1. Blank spaces interpretation

Standard Java provides two methods for determining whether a character is a blank space: Character​.isWhitespace(…) and Character​.isSpaceChar(…). These two methods differ in their interpretations of non-breaking spaces, tabs and line breaks. The first method conforms to the interpretation currently used in languages such as Java, C/C++ and XML, which considers tabs and line breaks to be blank spaces, while non-breaking spaces are read as not blank. The second method — which conforms strictly to the Unicode definition — makes the opposite interpretation.

SIS uses each of these methods in different contexts. isWhitespace(…) is used to separate the elements of a list (numbers, dates, words, etc.), while isSpaceChar(…) is used to ignore blank spaces inside a single element.

In practice, this distinction is reflected in the use of isSpaceChar(…) in the implementations of java.text.Format, or the use of isWhitespace(…) in nearly all the rest of the SIS library.

A geospatial information community is a collection of systems or individuals capable of exchanging their geospatial data through the use of common standards, allowing them to communicate with one another. As there are many ways to represent geospatial information, each community tends to structure this information in light of its areas of interest. This diversity complicates the task of Spatial Information System (SIS) users by confronting them with an apparently chaotic variety of data formats and structures. The characteristics of these structures vary according to the observed phenomenon and measurement methods, as well as the habits of the organizations producing the data. Such a variety represents an obstacle in studies that require heterogeneous combinations of data, especially when they originate in communities that are traditionally distinct. For example, a researcher studying cholera might be interested in populations of shrimp as a propagation vector of the disease. But as doctors and oceanographers may not be used to share their work, the participants of such a study may be limited by the effort required to convert the data.

We cannot impose a uniform format on all data collections, as the diversity of formats is tied to factors such as the constraints imposed by the measuring apparatus, and the statistical distribution of values. A more flexible solution is to ensure the interoperability of data across a common programming interface (API). This API is not necessarily defined in a programming language; the actual tendency is rather to define conventions that use existing web protocols, which we can translate into various programming languages. But in order for this approach to be viable, the API must be generally accepted by independent developers. In other words, the API must come as near as possible to industrial standards.

For example, one task that benefit from a successful standardization is the accessing of relational databases. The industry has established a common language — the SQL standard — that the creators of Java have embedded in standard JDBC programming interfaces. Today, these interfaces are implemented by many software programs, both free and commercial. Like databases, methods of accessing geographic information have been standardized. In this case, however, the efforts have been more recent, and their integration in software — especially in older programs — is incomplete and not always coherent. At the time of writing, no product to our knowledge has implemented all of the specifications in their entirety. However, there are many implementations that cover a fairly large spectrum. One of these is the Apache SIS™ library that is described in this document.

Apache SIS is characterized by a sustained effort to comply with standards. In general, complying with standards demands a greater effort than would be required for an isolated development, but rewards us with a double advantage: not only does it improve the interoperability of our data with that of external projects, it also points towards a robust way of elaborating the conceptual model reflected in the API. In effect, the groups of experts who conceived the standards anticipated difficulties that sometimes escape the engineer at the beginning of a project, but which risk to hit them before the end.

package org.opengis.referencing.crs;

/**
 * A 2D coordinate reference system used to approximate the shape of the earth on a planar surface.
 */
@UML(specification=ISO_19111, identifier="SC_ProjectedCRS")
public interface ProjectedCRS extends GeneralDerivedCRS {
    /**
     * Returns the coordinate system, which must be Cartesian.
     */
    @UML(obligation=MANDATORY, specification=ISO_19111, identifier="coordinateSystem")
    CartesianCS getCoordinateSystem();
}

Class<?> type = ProjectedCRS.class;
System.out.println("Standard name of type " + type.getName() + " is " + Types.getStandardName(type));

MediumName cdRom  = MediumName.CD_ROM;
MediumName usbKey = MediumName.valueOf("USB_KEY"); // There is no constraint on this value.
assert MediumName.valueOf("CD_ROM")  == cdRom  : "valueOf must return existing constants.";
assert MediumName.valueOf("USB_KEY") == usbKey : "valueOf must cache the previously requested values.";

Class<?> type   = Citation.class;
Method   method = type.getMethod("getTitle", (Class<?>[]) null);
UML      annot  = method.getAnnotation(UML.class);
String   id     = annot.identifier();
System.out.println("The standard name for the " + method.getName() + " method is " + id);

Properties isoToGeoAPI = new Properties();
try (InputStream in = UML.class.getResourceAsStream("class-index.properties")) {
    isoToGeoAPI.load(in);
}
String isoName = "CI_Citation";
String geoName = getProperty(isoName);
Class<?>  type = Class.forName(geoName);
System.out.println("The GeoAPI interface for ISO type " + isoName + " is " + type);

import org.opengis.referencing.GeodeticDatum;
import org.opengis.referencing.DatumFactory;
import java.util.ServiceLoader;

public class MyApplication {
    public void createMyDatum() {
        ServiceLoader  loader = ServiceLoader.load(DatumFactory.class);
        DatumFactory  factory = loader.iterator().next();
        GeodeticDatum myDatum = factory.createGeodeticDatum(…);
    }
}

In addition to its own tests, Apache SIS uses tests defined by GeoAPI. One advantages is that those tests provide an external source for the definition of expected results (for example the numerical values of coordinates obtained after a map projection). Such external source reduce the risk that some tests are actually anti-regression tests instead of correctness tests. Those tests can also be used by projects other than Apache SIS.

import org.opengis.metadata.Metadata;
import org.opengis.test.Validators;
import org.junit.Test;

public class MyTest {
    /*
     * Tolerate the absence of mandatory attributes in metadata and citation packages.
     * This modification applies to all tests executed by the currently running JVM.
     * If there are multiple test classes, this initialization may be performed
     * in a parent class to all test classes.
     */
    static {
        Validators.DEFAULT.metadata.requireMandatoryAttributes = false;
        Validators.DEFAULT.citation.requireMandatoryAttributes = false;
    }

    @Test
    public void testMyMetadata() {
        Metadata myObject = …; // Create an object here.
        Validators.validate(myObject);
    }
}

import org.opengis.metadata.Metadata;
import org.opengis.test.ValidatorContainer;
import org.junit.Test;

public class MyTest {
    private final ValidatorContainer validators;

    public MyTest() {
        validators = new ValidatorContainer();
        validators.metadata.requireMandatoryAttributes = false;
        validators.citation.requireMandatoryAttributes = false;
    }

    @Test
    public void testMyMetadata() {
        Metadata myObject = …; // Create an object here.
        validators.validate(myObject);
    }
}

import org.junit.*;
import org.junit.runner.RunWith;
import org.junit.runners.JUnit4;
import org.opengis.test.referencing.ParameterizedTransformTest;
import static org.junit.Assert.*;

@RunWith(JUnit4.class)
public class MyTest extends ParameterizedTransformTest {
    /**
     * Specify our own coordinate transformation factory for the GeoAPI tests.
     * GeoAPI will test the objects created by this factory.
     */
    public MyTest() {
        super(new MyMathTransformFactory());
    }

    /**
     * Changes the behaviour of a test. This example relaxes the requirements of this test a little,
     * by accepting errors of up to 10 centimetres, rather than the default value of 1 cm.
     * This change only applies to this method, and does not affect the other inherited tests.
     */
    @Test
    @Override
    public void testLambertAzimuthalEqualArea() throws FactoryException, TransformException {
        tolerance = 0.1;                    // 10 cm tolerance.
        super.testLambertAzimuthalEqualArea();
    }

    /**
     * Supplemental verification performed after each test, inherited or not.
     * In this example, we are verifying that the transformation tested
     * works correctly in two-dimensional spaces.
     */
    @After
    public void ensureAllTransformAreMath2D() {
        assertTrue(transform instanceof MathTransform2D);
    }
}

Following sections explain the rational behind some implementation choices done in Apache SIS.

8.3.3. Usages of Jacobian matrix

The Jacobian matrix (or transform derivative) was introduced in the section about coordinate operations. This section explains how derivatives are used by the Apache SIS implementation of some high-level operations.

8.3.3.1. Transform derivatives applied to envelopes

Geographic information systems often need to transform an envelope from one CRS to another. A previous section has shown why a naive approach transforming the 4 corners is not sufficient. A simple way to reduce the problem could be to sample a large number of points on each envelope border. For example we could sample 40 points on each border, transform them to the target CRS and compute an envelope containing all the transformed points. But even with 40 points, the point which is closest to the minimum in above figure can still be above that minimum. Consequently we still miss a small portion of the projected shape. It may be tempting to declare this error is negligible, but a single missing pixel can be enough for causing artifacts such as thin black lines in a mosaic or missing “piece of pie” in a projection over a pole. A workaround could be to arbitrarily increase the envelope size by some percentage, but a percentage too high will cause a larger amount of unnecessary data processing (e.g. more data to load from a file) while a percentage too low will leave some artifacts.

Map projection derivatives offer a more efficient way to resolve this problem. The figure below reproduces the same projected shape than above figure but with exaggerated deformations. The Apache SIS algorithm projects the 4 corners as in the naive approach, but opportunistically computes also their derivatives (the Jacobian matrices). Between two corners and using derivative information, there is only one possible cubic curve. We can describe that curve by the f(x) = C₀ + C₁x + C₂x² + C₃x³ equation with C coefficients computed from corner positions and derivatives. This approximation (shown by the red curve in above figure) does not match exactly the desired curve (shown by the blue shape in above figure) but is close. We do not use directly the cubic curve minimum, but instead we take the longitude λ (for above example) where this minimum is located. This is shown by the green dashed line in above figure. This position is easy to compute when the C coefficients are known. Assuming that the λ coordinate of the cubic curve minimum is close to the λ coordinates of the projected shape minimum, we can project the point on the envelope border at that location instead of projecting blindly 40 points on that border.

Legend:

• Blue: the geometric shape of the envelope after projection. This is the shape from which to get a new envelope.
• Red (with hash): The y = C₀ + C₁λ + C₂λ² + C₃λ³ approximation.
• Green (dashed line): Position λ_m of approximation minimum, obtained by resolving 0 = C₁ + 2C₂λ_m + 3C₃λ_m². The same cubic line can have two minimums.

In practice Apache SIS uses 8 points: the 4 corners plus one point at the center of each border of the envelope to project. The center points are added as a safety for map projections having symmetric deformations above and below equator. According our tests, using only those 8 points together with derivatives as described above gives more accurate results than “brute force” approach with 160 points on the 4 envelope borders. Saving 150 points seems small compared to computer performances. But above discussion used a two-dimensional envelope as an example. The same discussion applies to n-dimensional envelopes as well. Apache SIS can apply this algorithm on envelopes with any number of dimensions up to 64. The performance gain offered by this algorithm compared to “brute force” approach increases exponentially with the number of dimensions.

Apache SIS implements the algorithm described in this section by the static Envelopes​.transform(CoordinateOperation, Envelope) method. An alternative Envelopes​.transform(MathTransform, Envelope) method is also available, but should be used only when the CoordinateOperation is unknown. The reason is because MathTransform objects does not contain any information about coordinate system axes, which prevent the Envelopes​.transform(…) method to handle special cases such as envelopes containing a pole.

8.3.3.2. Transform derivatives applied to rasters

A raster can be projected by preparing an initially empty raster which will contain the resampled pixel values. Then for each pixel in the destination raster, we use the inverse of the map projection (T⁻¹) for computing the coordinates of the corresponding pixel in source raster. The transformed coordinates may not fall at the center of a source raster pixel, so source pixel values usually need to be interpolated.

Source image

Destination image

However applying the inverse transform on all pixel coordinates can be relatively slow. We can accelerate raster reprojections by pre-computing an interpolation grid. This grid contains the coordinate values of inverse transform T⁻¹ for only a few points. The coordinate values of other points are obtained by bilinear interpolations between interpolation grid points. Historically, this algorithm was implemented in the WarpGrid object of Java Advanced Imaging library. The performance gain is yet better if the same interpolation grid is reused for many rasters having their pixels at the same coordinates.

A difficulty in using interpolation grids is to determine how many points we need for keeping errors (defined as the difference between interpolated positions and actual positions computed by T⁻¹) below some threshold value (for example ¼ of pixel size). A “brute force” approach is to begin with a grid of size 3×3, then increase the number of points iteratively:

Legend:

• Blue dots: first iteration (9 points).
• Green dots: second iteration (25 points, including 16 news).
• Red dots: third iteration (81 points, including 56 news).

Continuing…

• Forth iteration: 289 points, including 208 news.
• Fifth iteration: 1089 points, including 800 news.
• Sixth iteration: 4225 points, including 3136 news.
• …

We could stop dividing the grid when, after having computed new points, we verified that the differences between coordinates computed by T⁻¹ during last iteration and coordinates computed by bilinear interpolations for the same points are lower than the threshold value. Unfortunately this approach can only tell us after computing new points… that those new points were unnecessary. This is unfortunate considering that each iteration adds an amount of points approximately equal to 3 times the sum of the number of points of all previous iterations.

Map projection derivatives help to improve the speed of interpolation grid computations. Using the derivatives, we can estimate whether another iteration is necessary before to execute it. The basic idea is to check if the derivatives at two neighbor points define two tangent lines that are almost parallel. In such case, we assume that the transform between those two points is almost linear. In order to define numerically the meaning of “almost linear”, we compute the intersection of the two tangent lines as illustrated below (the blue arrows). Then we compute the distance between that intersection and the straight line connecting the two points (the dashed line in the figure below).

In “brute force” algorithm without derivatives, iteration stops when the distance between interpolated positions (blue dashed line) and actual projected positions (red curve) is less than the threshold value. This criteria requires to compute the projected positions before the decision can be made. But with an algorithm using derivatives, the projected positions (red curve) are replaced by the tangents intersection point. If the actual projected curve does not have too much deformations (which should be the case for pairs of neighbor points close enough), then the red curve stay between the blue dashed line and the tangents intersection point. By using that intersection point, we reduce greatly the number of points to transform (3× the sum of all previous iterations).

8.3.3.3. Getting the derivative at a point

The algorithms discussed in previous section would be irrelevant if map projection derivatives were costly to compute. But analytic derivations of their formulas show that map projection and derivative formulas have many terms in common. Furthermore map projection formulas are simplified by Apache SIS implementation strategy, which takes out linear terms (delegated to affine transforms) so that NormalizedProjection subclasses can focus on the non-linear terms only. Map projection implementations in Apache SIS exploit those characteristics for computing Jacobian matrices (if required) together with projected points in a single operation, reusing many common terms between the two set of formulas. Example:

AbstractMathTransform projection = ...;         // An Apache SIS map projection.
double[] sourcePoint = {longitude, latitude};   // The geographic coordinate to project.
double[] targetPoint = new double[2];           // Where to store the projection result.
Matrix   derivative  = projection.transform(sourcePoint, 0, targetPoint, 0, true);

If only the Jacobian matrix is desired (without the projected point), then the MathTransform​.derivative(DirectPosition) method offers a more readable alternative.

Many operations on coordinate transforms can also be applied on transform derivatives: concatenation of a chain of transforms, inversion, taking a subset of the dimensions, etc. reverse projection formulas are usually more complicated than forward projections formulas, but thankfully their derivatives do not need to be computed: the Jacobian matrix of an inverse transform is the inverse of the Jacobian matrix of the forward transform. Calculation of reverse projections can be implemented as below:

@Override
public Matrix derivative(DirectPosition p) throws TransformException {
    Matrix jac = inverse().derivative(transform(p));
    return Matrices.inverse(jac);
}

MathTransform mt = ...;    // Any math transform created by Apache SIS.
if (mt instanceof AffineTransform) {
    AffineTransform at = (AffineTransform) mt;
    // Use Java2D API from here.
}