| Age | Commit message (Collapse) | Author |
|---|
|
utm: error out when value of +zone= is not an integer (fixes #2671)
|
|
NGO48 to ETRS89 TIN transfrm (#2554)
* Correction grid NKG:ETRF14 to EPSG:7922
* Added NKG:ITRF_TO_NO GIE test
* Correction grid no_kv_NKGETRF14_EPSG7922_2000 added to grid_alternatives.sql
* proj_method 'velocity_grid' added in check_grid_alternatives_proj_method. NKG velocity grid added to grid_alternatives.sql
* also include follow-up fix to master #2548: e10a849cd01c5b15e3123da1f6df845006322843
* Add NGO48 (EPSG:4273) to ETRS89 (EPSG:4258) triangulation-based transformation (#2554)
Co-authored-by: Even Rouault <even.rouault@spatialys.com>
|
|
This adds the NKG 2008 and 2020 transformations to proj.db. The NKG
transformations offers transformations between global reference frames
and the national realisations of ETRS89 in Denmark, Estonia, Finland,
Latvia, Lithuania, Norway and Sweden.
The 2008 transformations are already implemented in the NKG 2008 file
but will now be more accessible with the modern API.
The 2020 transformations are new to PROJ and offers and updated version
of the 2008 transformations using a new and improved deformation model
(eu_nkg_nkgrf17vel.tif). A 2020 version of the NKG transformations are
currently not available for Norway but will in all likelyhood be
included at a later point in time.
|
|
Fixes #2482
And also add proj_context_errno_string()
Revise gie 'expect failure errno XXXX' strings
|
|
Add +proj=topocentric geocentric->topocentric conversion (fixes #500)
|
|
The restriction was a copy&paste from the Evenden/Snyder approximate
ellipsoidal implementation, but the spherical one is exact, so this
restriction isn't needed.
Also tune a bit the handling of lat=0, |lon| > 90
|
|
(fixes #2468)
Corrected formula given by @evanmiller
|
|
|
|
Due to the improved accuracy of pj_tsfn(), it no longer returns
0 when phi=90° due to the conversion in radians.
Some GDAL tests are very sensitive to the pole transforming to (0,0)
exactly, so add a special case for that.
master only
|
|
Update Mercator projection, more accurate, faster
|
|
|
|
|
|
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.
|
|
by a number of projected CRS in Colombia (fixes #589)
|
|
|
|
Implement ellipsoidal formulation of +proj=ortho (fixes #397)
|
|
Implements RFC-6
|
|
to (0,0)
|
|
improve initial guessing
|
|
equatorial
|
|
|
|
- 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
|
|
Fixes #2333
Was due to 3dc92ad7014e1cf20a3bc95c8c21a34a371fcb78
Doesn't affect released versions
|
|
|
|
- Move the generic method initiated from adams_ws2 to a
pj_generic_inverse_2d() method
- Use it in adams_ws2
- Use it in wink2
Fixes https://github.com/qgis/QGIS/issues/35512
|
|
Limit peirce_q to northern hemisphere, and fix images for adams_hemi, guyou and peirce_q
|
|
- Central lon for zone 2 should be -d10, not d10
- Extra lobe was missing for zone 11
- New figure generated
- New test suite values generated
|
|
- The current implementation of the Interrupted Goode Homolosine
projection emphasizes land area. This is a compliment projection
that emphasizes ocean area.
- A value of lon0=-160 produces a reasonable real-world map.
|
|
Add a +proj=defmodel transformation for multi-component time-based deformation models
|
|
deformation models
Fixes #1001
Co-authored-by: Chris Crook <ccrook@linz.govt.nz>
|
|
Otherwise it returns junk (negative latitudes are mapped to
the same location as positive latitudes)
I'm a bit confused by PROJ peirce_q implementation. Looking at
the projection of the world, it looks like this matches the diamond
formulation of the right map shown at
https://desktop.arcgis.com/en/arcmap/latest/map/projections/peirce-quincuncial.htm,
but limited to the inner square of this diamond (which corresponds
to the northern hemisphere). We lack the 4 triangles on left, top,
right and bottom for the southern hemisphere. Furthermore, this
formulation of peirce_q does not seem to have the quincuncial
property of the square formulation (left images of the above ESRI
doc), or the one at https://en.wikipedia.org/wiki/Peirce_quincuncial_projection ...
|
|
|
|
In that mode:
* All non-comment/decoration lines must start with a valid tag
* Commands split on several lines should be terminated with " \"
|
|
tmerc/utm: add a +algo=auto/evenden_snyder/poder_engsager parameter
|
|
|
|
Add square conformal projections from libproject
|
|
This commit adds five new projections to PROJ:
adams_hemi: Adams Hemisphere in a Square
adams_wsI: Adams World in a Square I
adams_wsII: Adams World in a Square II
guyou: Guyou
peirce_q: Pierce Quincuncial
The code originates from Gerry Evendens libproject and has been adapted
to work with modern PROJ. To ensure that the modified code works as
intended extensive test data has been created using libproject and
sproj so that no errors occured when porting from libproject to PROJ.
The test data is wrapped in a gie files. All test cases reproduce
results from libproject at the mm level.
|
|
The default remains +alg=poder_engsager.
This default value can be changed in proj.ini
+algo=auto will use Evenden Synder implementation if the error in
doing so remains below 0.1 mm on Earth-sized ellipsoid
|
|
|
|
Saves 2 sincos() and 1 atan2() calls.
With the following bench
```
#include "proj.h"
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
int main(int argc, char* argv[])
{
if( argc != 2 ) {
fprintf(stderr, "Usage: bench_inv_cart fwd/inv\n");
exit(1);
}
PJ* p = proj_create(0, "+proj=cart");
const int dir = strcmp(argv[1], "inv") == 0 ? PJ_INV : PJ_FWD;
PJ_COORD coord;
if( dir == PJ_FWD )
{
coord.xyz.x = 3.14159/2;
coord.xyz.y = 3.14159/2;
coord.xyz.z = 100;
}
else
{
coord.xyz.x = 3e6;
coord.xyz.y = 3e6;
coord.xyz.z = 3e6;
}
for(int i = 0; i < 10* 1024 * 1024; i++ )
{
proj_trans(p, dir, coord);
}
proj_destroy(p);
return 0;
}
```
On Intel(R) Core(TM) i7-6700HQ CPU @ 2.60GHz
Time before: 2.37s
Time after: 1.57s
|
|
ENOMEM was wrongly set after setting PJD_ERR_ELLIPSOID_USE_REQUIRED
Note: it is a bit strange to forbid the pure spherical case whereas the
maths would allow it. I presume this is due to the typical usage of those
methods.
|
|
Fixes #1984
- Copy BETA2007.gsb, MD, alaska, conus, ntf_r93.gsb, ntv1_can.dat grids
from proj-datumgrid to data/tests.
- Replace a couple uses of nzgd2kgrid0005.gsb in tests by ntf_r93.gsb
- Add downsampled/subsetted versions of egm96_15.gtx as tests/egm96_15_downsampled.gtx
and ntv2_0.gsb as tests/ntv2_0_downsampled.gsb
This results in a few changes in expected results
- Simpify travis/install.sh due to less configurations to test
This results in a hopefully acceptable increase of the proj-X.Y.Z.tar.gz
from 2.9 to 5.3 MB
|
|
(fixes #1938)
And add testing of both little-endian and big-endian NTv2 files
|
|
Fix test issues on i386
|
|
Fix a few issues of #1906 found when running the test suite
on Ubuntu 16.04 with gcc 5.5 -m32. When applied on top of
the fix of #1912, make check succeeds
|
|
This should avoid issues with numerical stability as uncovered in
https://github.com/OSGeo/PROJ/issues/1906.
Practically speaking this change isn't going to affect real life
scenarios since the position of the center of the Earth is rarely
expressed in geodetic coordinates.
|
|
Fixes #1846
|
|
|
|
This commit is the result of the squashing of rfc4_dev branch in a single
commit. It implements mostly RFC 4 related work.
* Grid handling:
- remove obsolete and presumably unfinished implementation of grid catalog functionality
- all grid functionality is in grids.cpp/.hpp
- vertical and horizontal grid shift: rework to no longer load whole grid into memory
- remove hgrids and vgrids member from PJ structure, and store them in hgridshift/vgridshift/deformation structures
- build systems: add optional libtiff dependency. Must be explicitly disabled if not desired
- add support for horizontal and vertical grids in GeoTIFF, if libtiff is available
- add GenericShiftGridSet and GenericShiftGrid classes, relying on TIFF grids, that can be used for generic purpose grid-based adjustment
- add a +proj=xyzgridshift method to perform geocentric translation by grid. Used for French NTF to RGF93 transformation using gr3df97a.tif grid
- deformation: add support for +grids= for GeoTIFF grids
- horizontal grid shift: fix failures on points slightly outside a subgrid (fixes #209)
* File management:
- add a filemanager.cpp/.hpp to deal with file related work
- test for legacy proj_api.h fileapi
- proj.h: add proj_context_set_fileapi() and proj_context_set_sqlite3_vfs_name() (fixes #866)
- add capability to read resource files from the user writable directory
* Network access:
- build systems: add optional curl dependency
- add a curl-based default implementation for network related functionality
- proj.h: add C API to control network functionality, and optionaly provide network callbacks
- add data/proj.ini with default settings
- add a SQLite3 local cache of downloaded chunks
- add proj_is_download_needed() and proj_download_file()
* Use Win32 Unicode APIs and expect all strings to be UTF-8 (fixes #1765)
For backward compatibility, if PROJ_LIB content is found to be not UTF-8 or
pointing to a non existing directory, then an attempt at interpretating it
in the ANSI page encoding is done.
proj_context_set_search_paths() now assumes strings to be in UTF-8, and
functions returning paths will also return values in UTF-8.
|
|
Should fix the issue reported in https://lists.osgeo.org/pipermail/proj/2020-January/009188.html
Some extra north-american grids present in data/ can affect the results of
some tests, so create a data/for_tests/ subdirectory in which we copy only
select grids.
|
