Proj.4 versus Apache SIS: a performance comparison
- 28/08/2017
In August 17th 2017, we presented an introduction to Apache Spatial Information System (SIS) in the Free and Open Source Software for Geospatial (FOSS4G) conference. Discussions about performance and accuracy were planed, but skipped because of lack of time. A first part of the results that we intended to show are below.
Disclaimer: the author of this blog is an Apache SIS contributor.
Summary
The performance of two map projection libraries were compared:
- Proj.4, a library in the C language.
- Apache Spatial Information System (SIS), a library in the Java language.
The benchmark code was written in Java. The Proj.4 functions were invoked through the Java Native Interfaces (JNI). The use of JNI introduces a bias in benchmark measurements for Proj.4, but this bias is estimated lower than one standard deviation in the time measurements. In our results, Proj.4 is faster except for the Mercator inverse projection. The Proj.4 performance advantage is explained by the cost of Java 8 trigonometric functions like Math.sin(φ) and Math.asin(y)compared to their C counterparts. Java 9 is known to be faster than Java 8 but has not been tested yet. The Mercator inverse projection exception is explained by the mathematical work done in Apache SIS, where some formulas have been rewritten in more efficient ways using mathematical equivalences; in the Mercator case those gains outweigh the handicap of slower trigonometric functions. In a few extreme cases, Apache SIS is 400 times faster than Proj.4. Those extreme cases are explained by Apache SIS capability to detect when a chain of coordinate operations can be simplified as an affine transform. In all tested map projections, Proj.4 and Apache SIS results differ by a few micrometres or less. In datum shift tests, Proj.4 and Apache SIS are sometime in agreement and sometime apart by one or two meters. Those differences are explained by the way the two libraries use the EPSG geodetic dataset. More details are given in the discussion after “Material and method” section.
Material and method
The benchmark compared Proj.4 release 4.9.3 (August 2016) with Apache SIS 0.8-jdk8-SNAPSHOT (August 2017). The environment is Java 1.8.0_144-b01 on MacOS 10.12.6.
The tests use the GeoAPI interfaces. GeoAPI 3.0 is an OGC standard which allows to write code without knowledge of the underlying implementation. It is similar in this respect to Java DataBase Connectivity (JDBC) interfaces. This approach allows us to write benchmark or test codes only once, then execute it on arbitrary GeoAPI implementations. Apache SIS is one such implementations. The sis-gdal module provides another implementation as wrappers around the Proj.4 library (another variant is available in geoapi-proj4 module provided by the GeoAPI project).
Coordinate Reference System (CRS) instantiations
All Coordinate Reference System (CRS) objects are created from an EPSG code through the GeoAPI CRSAuthorityFactory interface. The Apache SIS CRS class provides convenience static methods which delegate to GeoAPI implementations for performing the real work. Coordinate Reference Systems backed by Proj.4 are created as below. Note that despite the “epsg” part, this is considered a Proj4 definition (indicated by the “Proj4::” prefix) rather than an EPSG definition. Those definitions differ in axis order and sometime in units of measurement.
CoordinateReferenceSystem crs = CRS.forCode(“Proj4::+init=epsg:4326”);Coordinate Reference Systems backed by Apache SIS are created as below. The first line creates a CRS as defined by EPSG. The second line modifies the CRS for the same axis order than Proj.4 (note that the CRS intentionally lost its “EPSG::4326” identifier in this process).
CoordinateReferenceSystem crs = CRS.forCode(“EPSG::4326”);
crs = AbstractCRS.castOrCopy(crs).forConvention(AxesConvention.RIGHT_HANDED);Note: alternatively, we could have created the Proj.4 CRS with “Proj4::+init=epsg:4326 +axis=neu”definition string. But this approach forces us to maintain a list of axis orientations for all supported EPSG codes. This is the approach implemented by GeoAPI wrappers for Proj.4, as a way to get an EPSG factory closer to authoritative definitions. Apache SIS takes a different approach where EPSG codes are handled by Apache SIS itself and the codes provided by Proj.4 are considered to be in a different, non-EPSG, name space.
Coordinate Operation instantiations
Many coordinate operations may exist for the same pair of source and target CRS. For example the EPSG geodetic dataset contains about 85 operations from NAD27 to WGS84, each operation using a different set of parameters for different state or geographic area. For making a choice, the Apache SIS CRS class provides another convenience method:
CoordinateOperation op = CRS.findOperation(sourceCRS, targetCRS, areaOfInterest);The areaOfInterest option is set to null in this test, which default to the widest domain of validity. Note that in the case of NAD27 to WGS84 transformations, the “widest domain of validity” criterion causes Apache SIS to select a transformation for Canada, not for USA. In order to avoid ambiguity, “NAD27 to WGS84” transformations are not compared in this benchmark. The issue of Proj.4 and Apache SIS selecting different transformations will be discussed in more details in another blog post, to be named “Proj.4 versus Apache SIS: an accuracy comparison”. For the performance comparison discussed in this post, suffice to say that we verified that the benchmark compares the same operation methods.
All coordinate conversions or transformations are executed through the GeoAPI MathTransforminterface. In the Proj.4 wrapper case, MathTransform delegates to Proj.4 pj_transformfunction. All those functions can work on arbitrary number of coordinates. Our benchmark uses this capability for applying coordinate operations on groups of 65536 coordinates. In Proj.4 case, this means one single JNI call and one single pj_transform call per group of 65536 coordinates. The cost of those calls has been estimated by performing an identity transform (where source and destination CRS are the same) on 65536 random coordinates with Proj.4 and assuming that all the execution time were due to JNI overhead. We measured 0.2 milliseconds, which is about 0.6% of the average execution time of other coordinate operations tested in this benchmark.
Conversions or transformations have been tested between the following pairs of CRS. Those pairs have been selected in order to use different operation methods: Cylindrical Equal Area (spherical and ellipsoidal cases), Pseudo-Mercator, Mercator, Lambert Conic Conformal, Polar Stereographic, Albers Equal Area, Geocentric translations and Position Vector transformation. This list does not include any Transverse Mercator (TM) test because there is at least 3 different TM approximations which are considered by EPSG as 3 different projection methods. Since Proj.4 and Apache SIS do not use the same approximation (Apache SIS uses the method recommended in EPSG guidance note 7-2), tests involving the TM projection would not be comparing the same thing.
| SOURCE CRS | TARGET CRS |
|---|---|
|
EPSG:4053 – International 1924 Authalic Sphere |
EPSG:3410 – NSIDC EASE-Grid Global |
|
EPSG:4326 – WGS 84 |
EPSG:6933 – WGS 84 / NSIDC EASE-Grid 2.0 Global |
|
EPSG:4326 – WGS 84 |
EPSG:3857 – WGS 84 / Pseudo-Mercator |
|
EPSG:4326 – WGS 84 |
EPSG:3395 – WGS 84 / World Mercator |
|
EPSG:4269 – NAD83 |
EPSG:3978 – NAD83 / Canada Atlas Lambert |
|
EPSG:4326 – WGS 84 |
EPSG:3031 – WGS 84 / Antarctic Polar Stereographic |
|
EPSG:4269 – NAD83 |
EPSG:5070 – NAD83 / Conus Albers |
|
EPSG:3994 – WGS 84 / Mercator 41 |
EPSG:3395 – WGS 84 / World Mercator |
|
EPSG:4301 – Tokyo |
EPSG:4612 – JGD2000 |
|
EPSG:4277 – OSGB 1936 |
EPSG:4230 – ED50 |
|
EPSG:4625 – Martinique 1938 |
EPSG:5489 – RGAF09 |
For each pair of CRS, a geographic bounding box is determined as below:
- If the EPSG geodetic dataset has an explicit coordinate operation for a given pair of CRS, then the domain of validity of that operation is used.
- Otherwise the intersection of the domain of validity of source and target CRS is used.
The geographic bounding box is divided in a grid of 256×256 equally-spaced cells. A coordinate is created at a random location inside each cell. Those coordinates are stored in an array of 65536 coordinates, then the coordinates are shuffled in random order. The result is named “original coordinates” and will not change anymore, in order to give us a stable reference for step 3 below. Then:
- Convert all coordinates from source CRS to target CRS (forward operation) with one single
MathTransformmethod call. Execution time is measured. - Convert all step 1 results from target CRS to source CRS (inverse operation) with one single
MathTransformmethod call. Execution time is measured. - The two above steps should cancel each other. This is verified by computing the average distance (in metres) between original coordinates and step 2 results, together with standard deviation. This distance is named “drift” in following discussion.
- Steps 1 to 3 are repeated 100 times (as recommended by Geospatial Integrity of Geoscience Software tests), using step 2 outputs as step 1 inputs. Execution times are ignored for the 10 first executions in order to give a chance to the Java Virtual Machine to “warm up” (i.e. compile the Java code). After 10 iterations, statistics on execution time are accumulated.
Note that the drift after 100 iterations is not necessarily equal to 100 times the drift after one iteration. For some map projections, the drift increases linearly as we would intuitively expect. But for some other projections, the drift increases asymptomatically. An example will be given in the “Proj.4 versus Apache SIS: an accuracy comparison” blog post.
Performances can be impacted by the precision of map projection algorithms. For example many map projections use iterative algorithms, where iterations are executed until the error falls below some arbitrary threshold value. Implementations with smaller threshold will iterate more often and consequently will appear slower (if everything else is equal). To be strict, we should compare implementations of equal precision. But for this benchmark we only verify that we are comparing implementations having at least a centimetric precision, as required by EPSG. We do not check the absolute precision (that would be the topic of other tests) but only the internal consistency as given by the drift measurement described at step 3 above.
Results and discussion
For all tested projections except Proj.4 Cylindrical Equal Area (ellipsoidal case), the drift after 100 iterations (i.e. after executing a “forward – inverse” cycle one hundred times) is less than one millimetre. The sole exception is the Proj.4 "cea" projection, which drifted by (6 ± 5) centimetres (for all measurements in this page, the value on right side of ± is one standard deviation. Assuming a gaussian distribution, (6 ± 5) cm means that about 68% of points drifted by a value between 1 and 11 cm). This drift occurs because the Cylindrical Equal Area inverse projection is approximated by a series expansion (φ ≈ c1⋅sin(2β) + c2⋅sin(4β) + c3⋅sin(8β)) which would need one more term. This is not Proj.4 “fault” since it implements the series expansion as published in both Snyder book (equation 3-18) and EPSG guidance note. The drift does not happen with Apache SIS because its high eccentricity support (for planetary CRS) fixed the imprecision as a side effect. That later topic will be discussed in the “Proj.4 versus Apache SIS: an accuracy comparison” blog post.
The performance comparison results is show below. Dark blue and red bars are Apache SIS average execution times; light blue and orange bars are Proj.4 average execution times. Smaller is better. Discussion follow.
The ⟝⟞ symbols at the tip of bars represent one standard deviation above and below the cumulated (forward + inverse) execution time. For Proj.4, execution time measurements have standard deviations ranging from 0.3 to 2 milliseconds. Since we estimated JNI cost to be less than 0.2 milliseconds, that cost appear to be smaller than our time measurement uncertainties.
Proj.4 appears faster than Apache SIS in most cases. The cause for this performance difference can be isolated more easily by looking at the Cylindrical Equal Area (Spherical) projection. This particular operation method contains a single sin(φ) in the forward projection case and a single asin(y) in the inverse projection case, ignoring multiplications and additions. A micro-benchmark invoking sin(φ) and asin(y) on 65536 random values in a loop, without any other arithmetic operations, reproduced approximatively the same execution time. This suggests that the Apache SIS execution time is mostly due to slower trigonometric functions in Java 8 compared to C. The same micro-benchmark on Java 9 release candidate 1 shows a sin(φ) function 3 to 4 time faster than in Java 8. The Java 9 performance improvement brings the forward projection at the same speed level than Proj.4. However the asin(y) performance stay the same in Java 9 than in Java 8. We could not yet test what will be the net performance gain on Java 9 since Apache SIS needs some work for Java 9 compatibility.
Note: the slower trigonometric functions in Java 8 come from stricter requirements. Java does not require the exact same result on all platforms, but requires the computed results to be within 1 ULP of exact results for most trigonometric functions. By contrast, the C/C++ language allows more flexibility in favour of performance. The precision requirement could be a reason why the Apache SIS drift is smaller than Proj.4 drift for map projections based on exact formulas, but this hypothesis has not been verified.
While Proj.4 is generally faster, there is some exceptions. The Mercator inverse projection is faster with Apache SIS than with Proj.4. The reason is that Proj.4 uses an iterative algorithm while Apache SIS uses a series expansion. That series expansion is not given explicitly in Snyder’s book, but can be inferred from the understanding that Mercator projection is the limit case of Lambert Conic Conformal projection when the standard parallels are at the equator. The difference between Apache SIS and Proj.4 results is less than a nanometre.
Note: series expansions are generally faster but their precision depends on the number of terms, which is hard-coded. This number of terms is usually selected for planets having an eccentricity not greater than Earth’s eccentricity. On the other hand, iterative algorithm are slower but can adapt to planets having greater eccentricity. Proj.4 implementation follows Snyder’s book, which gives an iterative algorithm for the Mercator projection and a series expansion for Lambert Conic Conformal. Apache SIS uses a mix of both: it starts with a series expansion for performance reasons, then continue with iterative algorithm if needed. This approach will be discussed in the future “Proj.4 versus Apache SIS: an accuracy comparison” blog post.
An extreme case is the conversion from WGS 84 / Mercator 41 (EPSG:3994) to WGS 84 / World Mercator (EPSG:3395). Those two projections differ only by the scale factor or latitude of standard parallel. In the Mercator projection case, a change of standard parallel can be mapped to a change of scale factor (such equivalence does not exist for every map projections). Consequently the “Mercator 41” to “World Mercator” conversion can be represented by the following chain of operations:
In this chain, the green boxes represent linear operations. We show only the scale factors for simplicity, but those green boxes can also contain longitude rotations, unit conversions, everything that can be represented by an affine transform. Having extracted linear operations out, the remaining non-linear parts in yellow boxes depend only on eccentricity (for the Mercator projection; other map projections have more non-linear terms). Since in this particular case nothing stands between the “inverse Mercator” and “Mercator” operations (i.e. there is no datum shift) and those two operations use the same eccentricity, they cancel each other:
By performing this decomposition, Apache SIS has been able to determine that the “WGS 84 / Mercator 41” to “WGS 84 / World Mercator” conversion simplifies to an affine transform. Such affine transform are much faster to execute than any map projection except Plate Carré; at least 400 times faster for the particular conversion tested here.
Apache SIS performs such decomposition not only for the Mercator projection, but for allcoordinate operations implemented in SIS: all map projections, all datum shift methods, any operation with any number of dimensions. In practice, there is not so many cases where a whole chain of operations can be simplified to an affine transform, so performance gain as spectacular as 400 times are not so common. However there is a lot of cases where some intermediate steps can be simplified, for example conversions from radians to degrees (e.g. the output of inverse projection) followed by a conversion from degrees to radians (e.g. the input of forward projection). Those simplifications do not always bring noticeable performance gains since saving a few multiplications has negligible effect compared to trigonometric function costs. But they contribute to reducing rounding errors, open the door to more optimizations in the future and make much easier the calculation of Jacobian matrices. The Jacobian matrices help to speedup higher-level operations discussed in previous blog posts (e.g. envelope transformations and raster resampling).
Conclusion
When operating on 65536 coordinates by method call, Proj.4 is generally faster than Apache SIS on Java 8. This difference in performance is due in large part to trigonometric functions. Faster trigonometric functions in Java 9 should help to reduce the performance gap in the general case, but not eliminate it since not all trigonometric functions have been improved. However there is cases where Apache SIS appears faster than Proj.4 despite the trigonometric function costs. Those cases can be explained by the mathematical analytic done in Apache SIS and the aggressive decomposition of some coordinate operation steps into (affine transform + non linear function) tupples. This Apache SIS architecture also gives other mathematical tools for more optimizations, like Jacobian matrices used for faster and more accurate envelope transformations.
Proj.4 and Apache SIS numerical results are in agreement for coordinate conversions (including map projections), but can differ by one or two metres for coordinate transformations (including datum shifts). Those differences are caused by different ways to use the EPSG geodetic dataset (early-binding versus late-binding approaches) and will be discussed in more details in another blog post, “Proj.4 versus Apache SIS: an accuracy comparison”.
