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.

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.

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

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