[Proj] Grown error if re-projecting from 4269 to LCC (2285) and backward multiple times

Wed Nov 1 13:43:51 PDT 2017

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
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/proj/attachments/20171101/cb7f9c58/attachment.html>