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

Dmitry Mefed dmitrymefed at gmail.com
Thu Nov 2 12:25:41 EST 2017


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
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.maptools.org/pipermail/proj/attachments/20171102/f7cfedef/attachment.htm 


More information about the Proj mailing list