tmerc/utm: add a +algo=auto/evenden_snyder/poder_engsager parameter

[rouault]

RFC for now as I guess this will raise questions/concerns. Would need documentation adjustments in particular if we decided for it.

Currently, tmerc and utm use the 'exact' algorithm (not so exact by
the way from what I could see from @cffk work), that is the Poder/Engsager
one. The price to pay for that is measurable: it is twice slower than
the approx Evenden/Snyder one, when used through proj_trans().

With the following bench.c

 #include "proj.h"
 #include <stdlib.h>
 #include <math.h>

 int main(int argc, char* argv[])
 {
    PJ* p = proj_create(0, argv[1]);
    PJ_COORD coord;
    coord.xy.x = atof(argv[2]) / 180. * M_PI;
    coord.xy.y = 0;
    for(int i = 0; i < 20 * 1000* 1000; i++)
        proj_trans(p, PJ_FWD, coord);
 }

./bench "+proj=utm +zone=31" 4 runs in 6.1s, whereas
with +approx it runs in 3.0s

With GRS80, at +/- 3 degrees far from the central meridian,
the difference between both algorithms is lesser than 0.1 mm:

echo "6 0" | src/proj -d 5 +proj=utm +zone=31 +approx
833978.55690	0.00000

$ echo "6 0" | src/proj -d 5 +proj=utm +zone=31
833978.55692	0.00000

I bet that most people will be happy with that, and operating
within +/- 3 deg of the central meridian is likely to cover 75% of
usages.

(If we allowed for an error of the order of 1mm, we could go up to +/- 6 deg,
but that's probably too much.)

Thus this PR implements an auto-selection mode of the algorithm,
based on the coordinate values.

In the forward path, Evenden/Snyder is selected if the distance to
the central meridian is below 3 degrees, otherwise it goes to Poder/Engsager.
In the reverse path, I chose 0.05 metre, which is a bit below that
the the projected distance of 3 degree on the unit sphere.

The cost of this extra test is within the measurement noise (likely < 1%)

Actually, I've added a +alg=evenden_snyder/poder_engsager/auto new
parameter to be able to precisely control the algorithm in cases where
it needs to be.

[cffk]

I believe the adjective describing the Poder/Engsager method is "extended" not "exact". Also in the code, it should be "ellipsoidal" not "elliptical".

[rouault]

method is "extended" not "exact"

I mentionned "exact", because this is what is mentioned currently in the code in function names etc to describe Poder/Engsager

Also in the code, it should be "ellipsoidal" not "elliptical".

In the doc you meant I guess. From grepping, I see this is used quite extensively in many pages of the doc with the "spherical and elliptical projection" expression. I guess a better phrasing would be just "spherical and ellipsoidal" (as this is after a "Available forms" header)

[busstoptaktik]

As @cffk mentions, the original meaning of etmerc was extended transverse Mercator, incidentally also celebrating one of the original authors of the code, by overloading the e to also refer to Karsten Engsager.

The word "exact" used in the comments in tmerc.cpp is unfortunate. It would make sense if we one day jump through the hoops of elliptical integrals, and also include the exact version by @cffk.

I think the exact/approx blunder may result from using "+approx" to indicate that the Snyder version should be used. Unfortunate, because both versions are approximate (truncated series expansions), but Snyder's version diverges faster than Poder/Engsager's extension of prior work by Krüger (1912) and König/Weise (1951).

We probably cannot change the name of +approx, but changing the comments to reflect this is less of a problem :-)

[cffk]

The criterion for using the approximate algorithm on the reverse
projection is too lax. E.g.,

echo 8e5 9e6 | ./proj -d 12 +proj=utm +zone=31 -I +alg=auto
  => 19.801327810051 80.669072129444
echo 8e5 9e6 | ./proj -d 12 +proj=utm +zone=31 -I +alg=poder_engsager
  => 19.801446664213 80.669070863629
echo 80.669072129444 19.801327810051 80.669070863629 19.801446664213 |
  ./geod -I +ellps=WGS84
  => 93d45'26.742" -86d14'32.835" 2.157

So the error is 2m. The problem is, of course, that the error in
approximate algorithm is a function of longitude from the central
meridian and not the distance.

My own preference would be for the default algorithms to yield close to
full double precision accuracy wherever possible. I know that 0.1 mm is
a tiny distance. But when proj is embedded in a larger system, little
errors can accumulate and cause failures which are difficult to track
down.

Even though speeding up UTM projections by a factor of 2 seems like a
big deal, 0.3 us is still very fast. So I'm left wondering whether this
is the bottleneck in many real life applications.

[rouault]

So the error is 2m. The problem is, of course, that the error in
approximate algorithm is a function of longitude from the central
meridian and not the distance.

Good point. Perhaps the criterion on the inverse path could be 0.05 * cos(lat) ~= 0.05 * (1 - 0.405 * lat *lat) to get something that evaluates fast and remains positive.
Admittedly, when combined with #2036 improvements, the speed difference between inverse Snyder / Extended methods shrinks dramatically (extended being then something like 16% slower than Snyder), so the above becomes less interesting.

However there's still a big gap between the forward paths that could make it still relevant: even after PR #2036, extended forward is 2.4x slower (with UTM31 and x=450e3 and y=5e6) (note: my initial measurement using y == 0 lead to a big bias, I believe due to a number of transcendent functions being evaluated very fast by glibc for very small arguments)

But when proj is embedded in a larger system, little errors can accumulate

Yes, but I'd tend to think they will become significant for real life only if you do something like hundeds or thousands of operations on the same point. Most workflow will probably do less than 10 steps on a given point.

So I'm left wondering whether this is the bottleneck in many real life applications.

There are applications where time spent in PROJ tmerc has been spotted as a bottleneck. For example raster reprojection with GDAL involving RPC georeferencing with DEM orthorectification, where the usual strategy of using piecewise linear coordinate transformation to reduce the number of coordinate transformation cannot be used, because of the irregularities of the DEM. In that case you need to perform one coordinate transformation for each pixel of the raster, and you can get then easily in the territory of several tens millions of such transformation per image. For folks processing thousands or more images per day, that starts showing in bills. I guess this is the same for foks using massive point clouds.
I'm not saying my current proposed implementation is the best regarding the threshold values, but something adaptative could probably be of interest.

[kbevers]

@rouault what's the status on this PR? Since you opened it you have made a bunch of other speedups to the transverse mercator projections but I reckon this is still relevant although I'm guessing the gains are a bit smaller now?

[rouault]

Approx fwd tmerc is still 2.06 times faster than extended. And inv is 1.61 times faster. But not sure we can find some consensus on the benefit (speed up) vs drawbacks (apparent "non determinism" and accuracy guarantee) of this approach. (I'm right now playing with SIMD computation of etmerc to compute several coordinates at the same time...)

$ time LD_LIBRARY_PATH=src/.libs ./bench "+proj=tmerc +approx" fwd

real	0m5.814s
user	0m5.812s
sys	0m0.000s

$ time LD_LIBRARY_PATH=src/.libs ./bench "+proj=tmerc" fwd

real	0m12.013s
user	0m12.000s
sys	0m0.000s

$ time LD_LIBRARY_PATH=src/.libs ./bench "+proj=tmerc +approx" inv

real	0m8.103s
user	0m8.092s
sys	0m0.000s

$ time LD_LIBRARY_PATH=src/.libs ./bench "+proj=tmerc" inv

real	0m13.067s
user	0m13.052s
sys	0m0.000s

[kbevers]

But not sure we can find some consensus on the benefit (speed up) vs drawbacks (apparent "non determinism" and accuracy guarantee) of this approach

There must be a way forward with this that can work for the majority of our users. I think this "adaptive algorithm" approach is a great idea. We could employ some settings (envvars or in proj.ini) that allow users with specific needs to configure the behavior of projection to their needs. A few ideas for that could be options for

1. always using a specific algorithm (force use of either tmerc or etmerc)
2. user definable threshold values

I think it should be possible to come up with a set up sensible default threshold values that will suit most needs but it is always nice for users to be able to control things completely. Personally I would set mine to "etmerc always on" but I have specific reasons for doing that.

I'm right now playing with SIMD computation of etmerc to compute several coordinates at the same time...

Very cool. Looking forward to seeing you results!

[cffk]

My vote (in case it wasn't clear already!) is that the default should be etmerc. This is reinforced by considering what the situation is likely to be in 5 or 10 years when error requirements are likely to be more stringent and computers are certainly going to be much faster.

[busstoptaktik]

@cffk said:

My vote (in case it wasn't clear already!) is that the default should be etmerc.

I fully agree: etmerc is extremely important for Greenland systems, where the entire country is mapped in one utm zone. So etmerc is important today, and as Charles say, probably even more so in 5 years.

[kbevers]

I fully agree: etmerc is extremely important for Greenland systems, where the entire country is mapped in one utm zone. So etmerc is important today, and as Charles say, probably even more so in 5 years.

I don't see how that is an argument against auto selection of the algorithm? For the Greenland case, coordinates near the UTM zone center would be projected using Snyder tmerc and coordinates further away using Poder/Engsager etmerc. Best of both worlds: fast and accurate in the center, slower but still accurate further away from the center. In this case etmerc will be used where it really shines and the faster tmerc is used where etmerc's extra series expansion terms are not needed.

[busstoptaktik]

I don't see how that is an argument against auto selection of the algorithm?

It is not. That was not what I was talking about. The thing is that, since the first derivatives are probably not continuous across the breakline, the auto-selection should be a (probably very attractive) opt-in.

[kbevers]

since the first derivatives are probably not continuous across the breakline

Right, I see. Let's check if that's the case then. Perhaps there's a fix for that if their indeed is a discontinuity.

[busstoptaktik]

indeed is a discontinuity

There surely is, and it will be terribly small. But you/we know from elevation model experience that discontinuities in the processing algorithm at way below millimeter level suddenly become visible in the large picture, if not smoothed over (i.e. forcing first derivatives to be contiuous).

So I'd prefer a Poder/Engsager default, with switches for +auto and +snyder, and would expect +auto to be widely used

[rouault]

I've force pushed an update for this pull request. Changes are:

• the default remains alg=poder_engsager
• this default can be changed through proj.ini
• the criterion for the inverse path has been improved to be dependent on the input y value

[busstoptaktik]

Looks good - but if you'll allow me to nitpick, I think "algo" will be more recognizable as a mnemonic for algorithm, than "alg" (cf. https://stackoverflow.com/questions/5372129/short-name-for-algorithm)

[rouault]

I think "algo" will be more recognizable

alg renamed to algo, and other similar references as well

[rouault]

@kbevers @cffk Does the latest changes look good to you ?

[kbevers]

@kbevers @cffk Does the latest changes look good to you ?

Yes, I am happy with your solution. Thanks for working through this!

[cffk]

@kbevers @cffk Does the latest changes look good to you ?

Yes, this is fine!