aboutsummaryrefslogtreecommitdiff
path: root/test/unit
AgeCommit message (Collapse)Author
2020-11-04Code formattingEven Rouault
2020-11-01Merge pull request #2371 from rouault/epsg10_part2Even Rouault
EPSG v10 update part 2: ingest DatumEnsemble from the database
2020-11-01proj.h: add PJ_CATEGORY_DATUM_ENSEMBLE for proj_create_from_database()Even Rouault
2020-11-01projinfo / createObjectsFromName(): support returning a datum ensembleEven Rouault
2020-11-01When reading from database, possibly return VerticalCRS with a DatumEnsembleEven Rouault
Only occurence for now is EPSG:9451 'BI height' using the 'British Isles height ensemble'
2020-11-01When reading from database, possibly return Geographic/GeodeticCRS with a ↵Even Rouault
DatumEnsemble, typically for WGS 84 and ETRS89 ('breaking change')
2020-11-01Merge pull request #2397 from cffk/merc-updateCharles Karney
Update Mercator projection, more accurate, faster
2020-10-31Fill remarks for PROJ-based operation mixing horizontal and vertical ↵Even Rouault
transformations
2020-10-27Use nm units in builtins.gie. Remove backward looking comments in code.Charles Karney
2020-10-26Address comments by @schwehrCharles Karney
2020-10-26Update Mercator projectionCharles Karney
Introduction ------------ The existing formulation for the Mercator projection is "satisfactory"; it is reasonably accurate. However for a core projection like Mercator, I think we should strive for full double precision accuracy. This commit uses cleaner, more accurate, and faster methods for computing the forward and inverse projections. These use the formulation in terms of hyperbolic functions that are manifestly odd in latitude psi = asinh(tan(phi)) - e * atanh(e * sin(phi)) (phi = latitude; psi = isometric latitude = Mercator y coordinate). Contrast this with the existing formulation psi = log(tan(pi/4 - phi/2)) - e/2 * log((1 + e * sin(phi)) / (1 - e * sin(phi))) where psi(-phi) isn't exactly equal to -psi(phi) and psi(0) isn't guaranteed to be 0. Implementation -------------- There's no particular issue implementing the forward projection, just apply the formulas above. The inverse projection is tricky because there's no closed form solution for the inverse. The existing code for the inverse uses an iterative method from Snyder. This is the usual hokey function iteration, and, as usual, the convergence rate is linear (error reduced by a constant factor on each iteration). This is OK (just) for low accuracy work. But nowadays, something with quadratic convergence (e.g., Newton's method, number of correct digits doubles on each iteration) is preferred (and used here). More on this later. The solution for phi(psi) I use is described in my TM paper and I lifted the specific formulation from GeographicLib's Math::tauf, which uses the same underlying machinery for all conformal projections. It solves for tan(phi) in terms of sinh(psi) which as a near identity mapping is ideal for Newton's method. For comparison I also look at the approach adopted by Poder + Engsager in their TM paper and implemented in etmerc. This uses trigonometric series (accurate to n^6) to convert phi <-> chi. psi is then given by psi = asinh(tan(chi)) Accuracy -------- I tested just the routines for transforming phi <-> psi from merc.cpp and measured the errors (converted to true nm = nanometers) for the forward and inverse mapping. I also included in my analysis the method used by etmerc. This uses a trigonometric series to convert phi <-> chi = atan(sinh(psi)), the conformal latitude. forward inverse max rms max rms old merc 3.60 0.85 2189.47 264.81 etmerc 1.82 0.38 1.42 0.37 new merc 1.83 0.30 2.12 0.31 1 nm is pretty much the absolute limit for accuracy in double precision (1 nm = 10e6 m / 2^53, approximately), and 5 nm is probably the limit on what you should routinely expect. So the old merc inverse is considerably less accurate that it could be. The old merc forward is OK on accuracy -- except that if does not preserve the parity of the projection. The accuracy of etmerc is fine (the truncation error of the 6th order series is small compared with the round-off error). However, situation reverses as the flattening is increased. E.g., at f = 1/150, the max error for the inverse projection is 8 nm. etmerc is OK for terrestrial applications, but couldn't be used for Mars. Timing ------ Here's what I get with g++ -O3 on various Linux machines with recent versions of g++. As always, you should take these with a grain of salt. You might expect the relative timings to vary by 20% or so when switching between compilers/machines. Times per call in ns = nanoseconds. forward inverse old merc 121 360 etmerc 4e-6 1.4 new merc 20 346 The new merc method is 6 times faster at the forward projection and modestly faster at the inverse projection (despite being more accurate). The latter result is because it only take 2 iterations of Newton's method to get full accuracy compared with an average of 5 iterations for the old method to get only um accuracy. A shocking aspect of these timings is how fast etmerc is. Another is that forward etmerc is streaks faster that inverse etmerc (it made be doubt my timing code). Evidently, asinh(tan(chi)) is a lot faster to compute than atan(sinh(psi)). The hesitation about adopting etmerc then comes down to: * the likelihood that Mercator may be used for non-terrestrial bodies; * the question of whether the timing benefits for the etmerc method would be noticeable in a realistic application; * need to duplicate the machinery for evaluating the coefficients for the series and for Clenshaw summation in the current code layout. Ripple effects ============== The Mercator routines used the the Snyder method, pj_tsfn and pj_phi2, are used in other projections. These relate phi to t = exp(-psi) (a rather bizarre choice in my book). I've retrofitted these to use the more accurate methods. These do the "right thing" for phi in [-pi/2, pi/2] , t in [0, inf], and e in [0, 1). NANs are properly handled. Of course, phi = pi/2 in double precision is actually less than pi/2, so cos(pi/2) > 0. So no special handling is needed for pi/2. Even if angles were handled in such a way that 90deg were exactly represented, these routines would still "work", with, e.g., tan(pi/2) -> inf. (A caution: with long doubles = a 64-bit fraction, we have cos(pi/2) < 0; and now we would need to be careful.) As a consequence, there no need for error handling in pj_tsfn; the HUGE_VAL return has gone and, of course, HUGE_VAL is a perfectly legal input to tsfn's inverse, phi2, which would return -pi/2. This "error handling" was only needed for e = 1, a case which is filtered out upstream. I will note that bad argument handling is much more natural using NAN instead of HUGE_VAL. See issue #2376 I've renamed the error condition for non-convergence of the inverse projection from "non-convergent inverse phi2" to "non-convergent sinh(psi) to tan(phi)". Now that pj_tsfn and pj_phi2 now return "better" results, there were some malfunctions in the projections that called them, specifically gstmerc, lcc, and tobmerc. * gstmerc invoked pj_tsfn(phi, sinphi, e) with a value of sinphi that wasn't equal to sin(phi). Disaster followed. I fixed this. I also replaced numerous occurrences of "-1.0 * x" by "-x". (Defining a function with arguments phi and sinphi is asking for trouble.) * lcc incorrectly thinks that the projection isn't defined for standard latitude = +/- 90d. This happens to be false (it reduces to polar stereographic in this limit). The check was whether tsfn(phi) = 0 (which only tested for the north pole not the south pole). However since tsfn(pi/2) now (correctly) returns a nonzero result, this test fails. I now just test for |phi| = pi/2. This is clearer and catches both poles (I'm assuming that the current implementation will probably fail in these cases). * tobmerc similarly thinks that phi close to +/- pi/2 can't be transformed even though psi(pi/2) is only 38. I'm disincline to fight this. However I did tighten up the failure condition (strict equality of |phi| == pi/2). OTHER STUFF =========== Testing ------- builtins.gei: I tightened up the tests for merc (and while I was about it etmerc and tmerc) to reflect full double precision accuracy. My test values are generated with MPFR enabled code and so should be accurate to all digits given. For the record, for GRS80 I use f = 1/298.2572221008827112431628366 in these calculations. pj_phi2_test: many of the tests were bogus testing irrelevant input parameters, like negative values of exp(-psi), and freezing in the arbitrary behavior of phi2. I've reworked most for the tests to be semi-useful. @schwehr can you review. Documentation ------------- I've updated merc.rst to outline the calculation of the inverse projection. phi2.cpp includes detailed notes about applying Newton's method to find tan(phi) in terms of sinh(psi). Future work ----------- lcc needs some tender loving care. It can easily (and should) be modified to allow stdlat = +/- 90 (reduces to polar stereographic), stdlat = 0 and stdlat_1 + stdlat_2 = 0 (reduces to Mercator). A little more elbow grease will allow the treatment of stdlat_1 close to stdlat_2 using divided differences. (See my implementation of the LambertConformalConic class in GeographicLib.) All the places where pj_tsfn and pj_phi2 are called need to be reworked to cut out the use of Snyder's t = exp(-psi() variable and instead use sinh(psi). Maybe include the machinery for series conversions between all auxiliary latitudes as "support functions". Then etmerc could use this (as could mlfn for computing meridional distance). merc could offer the etmerc style projection via chi as an option when the flattening is sufficiently small.
2020-10-23Database: add interpolation_crs_auth_name and interpolation_crs_code columns ↵Even Rouault
to other_transformation table
2020-10-23WKT1_ESRI export: try to export Geographic3D and Projected3D CRS when we can ↵Even Rouault
find a corresponding ellipsoidal vertical datum
2020-10-23WKT1_ESRI export: generate VERTCS[...,DATUM[...,SPHEROID[]] syntax when ↵Even Rouault
ellipsoidal height is found
2020-10-22WKT1_ESRI export: export CompoundCRS as PROJCS[...],VERTCS[...] or ↵Even Rouault
GEOGCS[...],VERTCS[...]
2020-10-22WKT parser: accept ESRI VERTCS[...,DATUM[...,SPHEROID[]] syntax to express ↵Even Rouault
ellipsoidal heights
2020-10-22WKT parser: accept implicit compoundCRS from ESRI WKT, like ↵Even Rouault
"PROJCS[...],VERTCS[...]"
2020-10-20VerticalCRS: morph CRS and datum name using ESRI aliases on import from / ↵Even Rouault
export to WKT1:ESRI
2020-10-20Improve identification of compound CRS from ESRI WKT1, and for compound CRS ↵Even Rouault
whose result is not in the DB but whose horiz and vertical parts are known
2020-10-20Orthographic projection: do not add f=0 to PROJ string if the ellipsoid is a ↵Even Rouault
sphere (fixes GDAL PDS4 tests)
2020-10-20test_c_api.cpp: reformatEven Rouault
2020-10-19C API: add proj_context_clone() (#2383)Alan D. Snow
Fixes #2382
2020-10-16Add multi-line PROJ string export capability, and use it by default in ↵Even Rouault
projinfo (unless --single-line is specified) (fixes #1543)
2020-10-16Merge pull request #2370 from rouault/epsg10Even Rouault
Update to EPSG 10.003 and make code base robust to dealing with WKT CRS with DatumEnsemble
2020-10-11C API: add proj_dynamic_datum_get_frame_reference_epoch()Even Rouault
2020-10-11Database query: add ↵Even Rouault
AuthorityFactory::ObjectType::DYNAMIC_GEODETIC_REFERENCE_FRAME and DYNAMIC_VERTICAL_REFERENCE_FRAME, and make corresponding C API work
2020-10-11Database: add a frame_reference_epoch column in vertical_datum to be able to ↵Even Rouault
handle dynamic vertical datums, and instanciate them properly from database
2020-10-10WKT2:2019 import/export: handle DATUM (at top level object) with PRIMEMEven Rouault
This is a peculiarity of the WKT grammar. Despite ISO 19111 saying that the prime meridian is a component of the datum, in WKT, they are placed at the same level, for backward compatibility with earlier WKT versions. So handle exporting and importing that. The fix is only for situation where DATUM is the top level object (was working fine otherwise), which is a uncommon use case. And to limit the amount of issue, on export emit the prime meridian only if it is not Greenwich.
2020-10-08Add C API to work with datum ensembleEven Rouault
Add: - proj_crs_get_datum_ensemble() - proj_crs_get_datum_forced() - proj_datum_ensemble_get_member_count() - proj_datum_ensemble_get_accuracy() - proj_datum_ensemble_get_member() Make proj_create_geographic_crs_from_datum() and proj_create_geocentric_crs_from_datum() accept a datum ensemble.
2020-10-08Make createOperations() work with DatumEnsembleEven Rouault
2020-10-08Make CRS identification work with CRS with DatumEnsembleEven Rouault
2020-10-08Add DatumEnsemble::asDatum() and use it in exportToWkt()Even Rouault
to allow exporting DatumEnsemble to WKT < 2019.
2020-10-08Add a AuthorityFactory::createDatumEnsemble() method, and make it inherit ↵Even Rouault
from ObjectUsage as mandated by ISO 19111:2019
2020-10-08Database: import datum ensemble accuracy and members (but do not use them)Even Rouault
2020-10-08getCRSInfoList(): make it use area descriptionEven Rouault
2020-10-08Database: use extended description for extent/area of use, as done by ↵Even Rouault
epsg.org WKT export
2020-10-06Database: instanciate DynamicGeodeticReferenceFrame (things like ITRFxxx, ↵Even Rouault
WGS 84 (Gxxxx), etc.) when possible
2020-10-06Database: add a reference_frame_epoch column to the geodetic_datum for ↵Even Rouault
dynamic datums, but not yet used
2020-10-06Database: "minimal" update to EPSG v10.003Even Rouault
Content mostly unchanged since v9.9 This update is "minimal" in that it mostly reflects the removal of the 'area' table, replaced now by 'extent', 'scope' and 'usage' Other new aspects of EPSG v10 are left aside.
2020-10-05proj_crs_create_bound_crs_to_WGS84(): make it work on ↵Even Rouault
verticalCRS/compoundCRS such as EPSG:4326+5773 and EPSG:4326+3855
2020-10-02createOperations(): avoid elimination of ballpark transformation that can ↵Even Rouault
help for NAD83->WGS84->NAD83(2011) hops
2020-09-30Merge pull request #2361 from rouault/ortho_ellipsoidalEven Rouault
Implement ellipsoidal formulation of +proj=ortho (fixes #397)
2020-09-30Merge pull request #2344 from rouault/tinshiftEven Rouault
Add a +proj=tinshift for triangulation-based transformations
2020-09-30Add a +proj=tinshift for triangulation-based transformationsEven Rouault
Implements RFC-6
2020-09-26Implement ellipsoidal formulation of +proj=ortho (fixes #397)Even Rouault
- Map ESRI 'Local' to +proj=ortho when Scale_Factor = 1 and Azimuth = 0 - Map ESRI 'Orthographic' to a PROJ WKT2 'Orthographic (Spherical)' which maps to +proj=ortho +f=0 to froce spherical evaluation
2020-09-17createOperations(): tune sorting of transformations so that the ones with ↵Even Rouault
greater 'version numbers' are prefered over other ones, when all other comparison criteria are equal. Helps with Amersfoort RD New to EPSG:4326
2020-09-17Adjust createBoundCRSToWGS84IfPossible() and operation filtering (for POSGAR ↵Even Rouault
2007 to WGS84 issues) (fixes #2356) - We make createBoundCRSToWGS84IfPossible() more restrictive. If there are more than one Helmert transformation from the CRS to WGS 84 covering the area of use of the CRS, we do not create a BoundCRS / +towgs84 - In createOperations() filtering, we are less aggressive in discarding operations that have the same area of use but worse accuracy. We do it only if they involve more transformation steps. We now get: ``` $ projinfo EPSG:5340 -o PROJ PROJ.4 string: +proj=longlat +ellps=GRS80 +no_defs +type=crs $ projinfo -s EPSG:5340 -t EPSG:4326 --spatial-test intersects --summary Candidate operations found: 2 EPSG:9264, POSGAR 2007 to WGS 84 (2), 0.5 m, Argentina EPSG:5351, POSGAR 2007 to WGS 84 (1), 1.0 m, Argentina ```
2020-09-09Database: update to EPSG 9.9Even Rouault
2020-08-26proj_create_vertical_crs_ex(): add a ACCURACY option to provide an explicit ↵Even Rouault
accuracy, or derive it from the grid name if it is known
2020-08-21createObjectsFromName(): add back case insensitivity when comparing namesEven Rouault
Fixes #2332 Didn't affect released versions
generated by cgit v1.2.3 (git 2.25.1) at 2025-12-15 19:23:15 +0000