[Proj] Grown error if re-projecting from 4269 to LCC (2285) and backward multiple times
Thomas Knudsen
knudsen.thomas at gmail.com
Thu Nov 2 14:55:12 PDT 2017
Hi Dmitry,
I googled around for the “alternative” Lambert, without luck, but then this
afternoon, I accidentally misspelled Lambert as “Lamber”, and ran into this
thread: http://lists.maptools.org/pipermail/proj/2003-March/000644.html
where PROJ.4 founder Gerald Evenden, and long time PROJ.4 maintainer Frank
Warmerdam discusses the LCCA, and Gerald mentions that it is a truncated
series, only introduced to be compatible with some legacy French/North
African systems, that were defined using the truncated series.
Gerald explicitly states that it should not be used for new work, so my
recommendation of looking at it was terribly bad advice. Also, I had not
realized that LCCA does not support the secant (“2 parallel”) case, and
hence lat_0 and not lat_1 & lat_2 are used for parameterizing the
projection.
So enough about LCCA - it should not be used for anything but to maintain
interoperability with legacy reference systems.
But I have been doing some timing experiments with a stripped down version
of the “real” LCC code. Varying the “TOL” parameter in pj_phi2.c. Using the
setup string
+proj=lcc +lat_1=33 +lat_2=45 +lat_0=33 +lon_0=-90 +x_0=0 +y_0=0
+ellps=GRS80 +no_defs,
I ran roundtrip tests for each 0.1 deg covering a latitude range of 10 N to
68 N, and a longitude range of 135 W to 45 W (522000 points), and got the
following connection between TOL, average number of iterations, average
duration for one projection call, and the maximum roundtrip difference
encountered:
TOL=1e-10: avg iterations=3.64, dt=1716 ns, maxdiff=2000 nm
TOL=1e-12: avg iterations=4.48, dt=1923 ns, maxdiff=28 nm
TOL=1e-13: avg iterations=4.75, dt=1977 ns, maxdiff=5 nm
Setting TOL lower than 1e-13 resulted in longer running time, but not
better maximum roundtrip difference, but the extra 50 ns for a fivefold
improvement of the worst case still seems worthwhile, compared to the 200
ns giving a 70-fold improvement.
So I find it well argued, that the choice of tolerance parameter in
pj_phi2.c must be considered a bug, and corrected: The enormous improvement
in roundtrip precision comes at a time cost of 15% on the stripped down
code, and a bit less for the code running inside PROJ.4, due to plumbing
overhead (which will, however approximately balance out the fact that I’m
measuring roundtrips above, not just inverse calls. But forward calls are
approximately 5 times faster than inverse - so give and take a bit, 15% is
not a bad estimate).
Also, we do not change anything in the forward projection, only improve the
consistency between forward and inverse. So this improvement should, all
other things being equal, be an obvious candidate for implementation.
But all other things are not necessarily entirely equal: The pj_phi2
function is also called from the Mercator variants calcofi, gstmerc, omerc,
and merc, and I have not tested the timing impact on them (their precision
is not altered - the test suite completes nicely).
I also took a look at the histograms of the number of iterations for the
1e-10 and 1e-13 cases. Interestingly they were largely unpopulated: The
1e-10 case was alive at 3 iterations (40% of calls) and 4 (60% of all
calls), while the 1e-13 case had approx 50/50 split between 4 and 5
iterations.
No cases of (1,2), resp. (1,2,3) could indicate that the initial guess
could be improved, but I have not found any obvious way of doing it (except
for things that really amounts to moving one iteration out of the loop,
which is hardly interesting). If you have any ideas, I would be very
interested.
/Thomas
2017-11-02 18:25 GMT+01:00 Dmitry Mefed <dmitrymefed at gmail.com>:
> Thank you Thomas,
>
> I've tried both solutions. The one with strengthen the criterion does work
> for me. Now I can push changes to the DB and query it back with no changes
> to the coordinates with maximum ACAD precision.
> As for LCCA it does not work for us, as in our example it produces a bit
> different results, and it is critical for us. Please see the below output:
>
> LCC EPSG:2285
> cs2cs.exe -v +init=epsg:4269 +to =+proj=lcc +lat_1=48.73333333333333
> +lat_2=47.5 +lat_0=47 +lon_0=-120.8333
> 333333333 +x_0=500000.0001016001 +y_0=0 +ellps=GRS80 +datum=NAD83
> +to_meter=0.3048006096012192 +no_defs -f %.16f
> -118.5293900000000300 48.7408309999023930
> 2196293.3184076622000000 643350.3929806272500000 0.0000000000000000
>
> LCCA same as above but lcca
> cs2cs.exe -v +init=epsg:4269 +to =+proj=lcca +lat_1=48.73333333333333
> +lat_2=47.5 +lat_0=47 +lon_0=-120.8333
> 333333333 +x_0=500000.0001016001 +y_0=0 +ellps=GRS80 +datum=NAD83
> +to_meter=0.3048006096012192 +no_defs -f %.16f
> -118.5293900000000300 48.7408309999860520
> 2196554.2242959067000000 643311.0110160053000000 0.0000000000000000
>
> Thank you all for your help.
>
> PS
> I did not observed any sensible performance issues with updated criterion
> value. However I did not do precise measurements.
>
> Thanks,
> -Dmitry Mefed
>
> On Wed, Nov 1, 2017 at 11:43 PM, Thomas Knudsen <knudsen.thomas at gmail.com>
> wrote:
>
>> Dmitry,
>>
>> No, it's not the EPS10 define in PJ_lcc.c - it's the #define TOL 1e-10 in
>> pj_phi2.c you need to modify.
>>
>> /Thomas
>>
>> 2017-11-01 19:47 GMT+01:00 Dmitry Mefed <dmitrymefed at gmail.com>:
>>
>>> Thank you all for the interesting discussion.
>>>
>>> I'm going to try both suggestions:
>>> 1. try to use LCCA instead of LCC
>>> 2. strengthen convergence criterion, I believe it is EPS10 define in
>>> case of LCC
>>>
>>> Thanks,
>>> -Dmitry
>>>
>>> On Wed, Nov 1, 2017 at 4:33 PM, Thomas Knudsen <knudsen.thomas at gmail.com
>>> > wrote:
>>>
>>>> Hi Oscar,
>>>>
>>>> Regarding LCCA:
>>>>
>>>> I can see that the functions for computing dr and its derivative can be
>>>> easily inlined, by defining them as macros, rather than functions.
>>>>
>>>> The call to fabs() is probably, on most platforms, converted to a
>>>> builtin,
>>>> so I would not touch that.
>>>>
>>>> The meridional arc length function, pj_mlfn(), is fairly simple and
>>>> COULD
>>>> be inlined, but at the cost of future agony if the function is one day
>>>> updated,
>>>> so I would prefer to keep it as is until actual profiling shows it to
>>>> be a
>>>> serious problem.
>>>>
>>>> Are there any other obvious speed problems in PJ_lcca (ignoring
>>>> everything
>>>> under the PROJECTION() call, which is setup functionality, executed only
>>>> once, and hence not so speed critical).
>>>>
>>>>
>>>> Regarding LCC:
>>>>
>>>> The iteration you mention goes on in the pj_phi2() function in
>>>> pj_phi2.c.
>>>>
>>>> By strengthening the convergence criterion to 1e-12, the 1000 iteration
>>>> roundtrip deviation for LCC was reduced to 0.005 mm, which certainly
>>>> is more acceptable.
>>>>
>>>> By strengthening to 1e-13, the 1000 roundtrip deviation was below
>>>> 0.0001 mm.
>>>>
>>>> I have, however, not yet had time to measure what the time cost of this
>>>> improvement
>>>> amounts to. Once I have hard numbers, I will post them to the list, and
>>>> probably suggest that we switch - either from LCC to LCCA, or from stop
>>>> criterion of 1e-10 to 1e-12.
>>>>
>>>> /Thomas
>>>>
>>>>
>>>> 2017-11-01 12:29 GMT+01:00 vanadovv at hetnet.nl <vanadovv at hetnet.nl>:
>>>>
>>>>> Hi Thomas,
>>>>>
>>>>>
>>>>> LCCA: good suggestion!
>>>>> As far as I can tell, there are two significant differences between
>>>>> LCC and LCCA.
>>>>> LCC has an iteration loop criterion (in the inverse) of 1e-10, whereas
>>>>> LCCA is somewhat more accurate with 1e-12.
>>>>> Furthermore LCCA works with a Newton iteration scheme. This could be
>>>>> faster than iteration by successive approximation, but there are a couple
>>>>> of minor inefficiencies in the code, like function calls instead of inline,
>>>>> but YMMV.
>>>>>
>>>>>
>>>>> Oscar van Vlijmen
>>>>>
>>>>>
>>>>>
>>>>> ----Origineel Bericht----
>>>>> Van : knudsen.thomas at gmail.com
>>>>> Datum : 01/11/2017 07:58
>>>>> Aan : vanadovv at hetnet.nl, proj at lists.maptools.org
>>>>> Onderwerp : Re: [Proj] Grown error if re-projecting from 4269 to LCC
>>>>> (2285) and backward multiple times
>>>>>
>>>>> Oscar,
>>>>>
>>>>> I certainly agree that the deviation is quite high - my comment was
>>>>> more related
>>>>> to the expection of exact roundtrips, which I find unrealistic.
>>>>>
>>>>> Nevertheless, looking into the PROJ.4 code, I see there is an
>>>>> alternative
>>>>> implementation of LCC, called LCCA, which I had forgotten about, which
>>>>> actually seems to roundtrip exactly.
>>>>>
>>>>> In the master branch of PROJ.4, over at
>>>>> https://github.com/OSGeo/proj.4,
>>>>> and comming in the next release, there is a tool "gie" for doing
>>>>> (a.o.)
>>>>> roundtrip tests. It is still lacking in docs, but I think you can
>>>>> follow this example:
>>>>>
>>>>> With this input file:
>>>>>
>>>>> $ cat lcc-lcca.gie
>>>>>
>>>>> BEGIN
>>>>>
>>>>> -------------------------------------------------------------------------------
>>>>>
>>>>> operation +proj=lcc +lat_1=48.73333333333333 +lat_2=47.5 \
>>>>> +lat_0=47 +lon_0=-120.8333333333333 \
>>>>> +x_0=500000.0001016001 +y_0=0 \
>>>>> +units=us-ft +ellps=GRS80 +towgs84=0,0,0 +no_defs
>>>>> -------------------------------------------------------------------------------
>>>>>
>>>>> tolerance 0.0010 mm
>>>>> accept -118.5293900000000300 48.7408309999860520
>>>>> roundtrip 1000
>>>>>
>>>>>
>>>>> -------------------------------------------------------------------------------
>>>>>
>>>>> operation +proj=lcca +lat_1=48.73333333333333 +lat_2=47.5 \
>>>>> +lat_0=47 +lon_0=-120.8333333333333 \
>>>>> +x_0=500000.0001016001 +y_0=0 \
>>>>> +units=us-ft +ellps=GRS80 +towgs84=0,0,0
>>>>> +no_defs
>>>>> -------------------------------------------------------------------------------
>>>>>
>>>>> tolerance 0.0 mm
>>>>> accept -118.5293900000000300 48.7408309999860520
>>>>> roundtrip 1000
>>>>> END
>>>>>
>>>>> I get this output:
>>>>>
>>>>> $ gie lcc-lcca.gie
>>>>>
>>>>> -------------------------------------------------------------------------------
>>>>>
>>>>> Reading file '..\..\..\proj\lcc-lcca.gie'
>>>>> -------------------------------------------------------------------------------
>>>>>
>>>>> +proj=lcc +lat_1=48.73333333333333 +lat_2=47.5 +lat_0=47
>>>>> +lon_0=-120....
>>>>> -------------------------------------------------------------------------------
>>>>>
>>>>> FAILURE in lcc-lcca.gie(11):
>>>>> roundtrip deviation: 1.550 mm, expected: 0.001 mm
>>>>> -------------------------------------------------------------------------------
>>>>>
>>>>> total: 1 tests succeeded, 1 tests FAILED!
>>>>> -------------------------------------------------------------------------------
>>>>>
>>>>>
>>>>> i.e. lcca roundtrips exactly, while lcc diverges heavily.
>>>>>
>>>>> So Dmitry: This seems to be your workaround - define your projection
>>>>> using lcca,
>>>>> rather than letting the epsg-list select lcc for you.
>>>>>
>>>>> /Thomas
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>> _______________________________________________
>>>>> Proj mailing list
>>>>> Proj at lists.maptools.org
>>>>> http://lists.maptools.org/mailman/listinfo/proj
>>>>>
>>>>
>>>>
>>>> _______________________________________________
>>>> Proj mailing list
>>>> Proj at lists.maptools.org
>>>> http://lists.maptools.org/mailman/listinfo/proj
>>>>
>>>
>>>
>>> _______________________________________________
>>> Proj mailing list
>>> Proj at lists.maptools.org
>>> http://lists.maptools.org/mailman/listinfo/proj
>>>
>>
>>
>> _______________________________________________
>> Proj mailing list
>> Proj at lists.maptools.org
>> http://lists.maptools.org/mailman/listinfo/proj
>>
>
>
> _______________________________________________
> Proj mailing list
> Proj at lists.maptools.org
> http://lists.maptools.org/mailman/listinfo/proj
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/proj/attachments/20171102/3ff3638c/attachment.html>
More information about the Proj
mailing list