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

[Dmitry Mefed]

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