[Proj] Grown error if re-projecting from 4269 to LCC (2285) and backward multiple times
Dmitry Mefed
dmitrymefed at gmail.com
Mon Nov 6 09:11:52 PST 2017
Hi Thomas,
All you said about number of iterations sounds interesting to me, but I do
not have projections algorithms background, so I will not be able to help
here in any way.
Thank you all for your help.
-Dmitry Mefed
On Fri, Nov 3, 2017 at 12:55 AM, Thomas Knudsen <knudsen.thomas at gmail.com>
wrote:
> 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
>>
>
>
> _______________________________________________
> 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/20171106/e8094fd5/attachment.html>
More information about the Proj
mailing list