<div dir="ltr"><div><div><div>Hi Thomas,<br><br></div>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.<br></div><br>Thank you all for your help.<br><br></div>-Dmitry Mefed<br></div><div class="gmail_extra"><br><div class="gmail_quote">On Fri, Nov 3, 2017 at 12:55 AM, Thomas Knudsen <span dir="ltr"><<a href="mailto:knudsen.thomas@gmail.com" target="_blank">knudsen.thomas@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div>Hi Dmitry,</div><div><br></div><div>I googled around for the “alternative” Lambert, without luck, but then this afternoon, I accidentally misspelled Lambert as “Lamber”, and ran into this thread: <a href="http://lists.maptools.org/pipermail/proj/2003-March/000644.html" target="_blank">http://lists.maptools.org/<wbr>pipermail/proj/2003-March/<wbr>000644.html</a> 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.</div><div><br></div><div>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.</div><div><br></div><div>So enough about LCCA - it should not be used for anything but to maintain interoperability with legacy reference systems.</div><div><br></div><div>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 </div><div><br></div><div>+proj=lcc +lat_1=33 +lat_2=45 +lat_0=33 +lon_0=-90 +x_0=0 +y_0=0 +ellps=GRS80 +no_defs,</div><div><br></div><div>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:</div><div><br></div><div>TOL=1e-10: avg iterations=3.64, dt=1716 ns, maxdiff=2000 nm</div><div>TOL=1e-12: avg iterations=4.48, dt=1923 ns, maxdiff=28 nm</div><div>TOL=1e-13: avg iterations=4.75, dt=1977 ns, maxdiff=5 nm</div><div><br></div><div>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.</div><div><br></div><div>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).</div><div><br></div><div>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.</div><div><br></div><div>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).</div><div><br></div><div>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.</div><div><br></div><div>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.</div><div><br></div><div>/Thomas</div><div><br></div></div><div class="gmail_extra"><br><div class="gmail_quote">2017-11-02 18:25 GMT+01:00 Dmitry Mefed <span dir="ltr"><<a href="mailto:dmitrymefed@gmail.com" target="_blank">dmitrymefed@gmail.com</a>></span>:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div><div><div><div><div><div><div>Thank you Thomas,<br><br></div>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.<br></div>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:<br><br><font size="1"><span style="font-family:monospace,monospace">LCC EPSG:2285<br>cs2cs.exe -v +init=epsg:4269 +to =+proj=lcc +lat_1=48.73333333333333 +lat_2=47.5 +lat_0=47 +lon_0=-120.8333<br>333333333 +x_0=500000.0001016001 +y_0=0 +ellps=GRS80 +datum=NAD83 +to_meter=0.3048006096012192 +no_defs -f %.16f<br>-118.5293900000000300 48.7408309999023930<br>2196293.3184076622000000 <wbr> 643350.3929806272500000 0.0000000000000000<br><br>LCCA same as above but lcca<br>cs2cs.exe -v +init=epsg:4269 +to =+proj=lcca +lat_1=48.73333333333333 +lat_2=47.5 +lat_0=47 +lon_0=-120.8333<br>333333333 +x_0=500000.0001016001 +y_0=0 +ellps=GRS80 +datum=NAD83 +to_meter=0.3048006096012192 +no_defs -f %.16f<br>-118.5293900000000300 48.7408309999860520<br>2196554.2242959067000000 <wbr> 643311.0110160053000000 0.0000000000000000</span></font><br><br></div>Thank you all for your help.<br><br></div>PS<br></div>I did not observed any sensible performance issues with updated criterion value. However I did not do precise measurements.<br><br></div>Thanks,<br></div>-Dmitry Mefed <br></div><div class="m_1973532511232938381HOEnZb"><div class="m_1973532511232938381h5"><div class="gmail_extra"><br><div class="gmail_quote">On Wed, Nov 1, 2017 at 11:43 PM, Thomas Knudsen <span dir="ltr"><<a href="mailto:knudsen.thomas@gmail.com" target="_blank">knudsen.thomas@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr">Dmitry,<div><br></div><div>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.</div><div><br></div><div>/Thomas</div></div><div class="gmail_extra"><br><div class="gmail_quote">2017-11-01 19:47 GMT+01:00 Dmitry Mefed <span dir="ltr"><<a href="mailto:dmitrymefed@gmail.com" target="_blank">dmitrymefed@gmail.com</a>></span>:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div><div><div><div><div>Thank you all for the interesting discussion.<br><br></div>I'm going to try both suggestions:<br></div>1. try to use LCCA instead of LCC<br></div>2. strengthen convergence criterion, I believe it is EPS10 define in case of LCC<br><br></div>Thanks,<br></div>-Dmitry<br></div><div class="m_1973532511232938381m_-8205083118294684756m_930862448995323276HOEnZb"><div class="m_1973532511232938381m_-8205083118294684756m_930862448995323276h5"><div class="gmail_extra"><br><div class="gmail_quote">On Wed, Nov 1, 2017 at 4:33 PM, Thomas Knudsen <span dir="ltr"><<a href="mailto:knudsen.thomas@gmail.com" target="_blank">knudsen.thomas@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr">Hi Oscar,<div><br></div><div>Regarding LCCA:</div><div><br></div><div>I can see that the functions for computing dr and its derivative can be</div><div>easily inlined, by defining them as macros, rather than functions.</div><div><br></div><div>The call to fabs() is probably, on most platforms, converted to a builtin,</div><div>so I would not touch that.</div><div><br></div><div>The meridional arc length function, pj_mlfn(), is fairly simple and COULD</div><div>be inlined, but at the cost of future agony if the function is one day updated,</div><div>so I would prefer to keep it as is until actual profiling shows it to be a</div><div>serious problem.</div><div><br></div><div>Are there any other obvious speed problems in PJ_lcca (ignoring everything</div><div>under the PROJECTION() call, which is setup functionality, executed only</div><div>once, and hence not so speed critical).</div><div><br></div><div><br></div><div>Regarding LCC:</div><div><br></div><div>The iteration you mention goes on in the pj_phi2() function in pj_phi2.c.</div><div><br></div><div>By strengthening the convergence criterion to 1e-12, the 1000 iteration</div><div>roundtrip deviation for LCC was reduced to 0.005 mm, which certainly</div><div>is more acceptable.</div><div><br></div><div>By strengthening to 1e-13, the 1000 roundtrip deviation was below 0.0001 mm.</div><div><br></div><div>I have, however, not yet had time to measure what the time cost of this improvement</div><div>amounts to. Once I have hard numbers, I will post them to the list, and</div><div>probably suggest that we switch - either from LCC to LCCA, or from stop</div><div>criterion of 1e-10 to 1e-12.</div><div><br></div><div>/Thomas </div><div><br></div></div><div class="gmail_extra"><br><div class="gmail_quote">2017-11-01 12:29 GMT+01:00 <a href="mailto:vanadovv@hetnet.nl" target="_blank">vanadovv@hetnet.nl</a> <span dir="ltr"><<a href="mailto:vanadovv@hetnet.nl" target="_blank">vanadovv@hetnet.nl</a>></span>:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">Hi Thomas,<div><br></div><div><br></div><div>LCCA: good suggestion!</div><div>As far as I can tell, there are two significant differences between LCC and LCCA.</div><div>LCC has an iteration loop criterion (in the inverse) of 1e-10, whereas LCCA is somewhat more accurate with 1e-12.</div><div>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.</div><div><br></div><div><br></div><div>Oscar van Vlijmen</div><div><br></div><div><br></div><div><br></div><div><blockquote style="margin-right:0px;margin-left:15px">----Origineel Bericht----<br>Van : <a href="mailto:knudsen.thomas@gmail.com" target="_blank">knudsen.thomas@gmail.com</a><br>Datum : 01/11/2017 07:58<br>Aan : <a href="mailto:vanadovv@hetnet.nl" target="_blank">vanadovv@hetnet.nl</a>, <a href="mailto:proj@lists.maptools.org" target="_blank">proj@lists.maptools.org</a><span><br>Onderwerp : Re: [Proj] Grown error if re-projecting from 4269 to LCC (2285) and backward multiple times<br><br></span><div><div class="m_1973532511232938381m_-8205083118294684756m_930862448995323276m_-1897303595459080174m_7246666681453626205h5"><div dir="ltr">
Oscar,
<div>
<br>
</div>
<div>
I certainly agree that the deviation is quite high - my comment was more related
</div>
<div>
to the expection of exact roundtrips, which I find unrealistic.
</div>
<div>
<br>
</div>
<div>
Nevertheless, looking into the PROJ.4 code, I see there is an alternative
</div>
<div>
implementation of LCC, called LCCA, which I had forgotten about, which
</div>
<div>
actually seems to roundtrip exactly.
</div>
<div>
<br>
</div>
<div>
In the master branch of PROJ.4, over at
<a href="https://github.com/OSGeo/proj.4" target="_blank">https://github.com/OSGeo/proj.<wbr>4</a>,
</div>
<div>
and comming in the next release, there is a tool "gie" for doing (a.o.)
</div>
<div>
roundtrip tests. It is still lacking in docs, but I think you can follow this example:
</div>
<div>
<br>
</div>
<div>
With this input file:
</div>
<div>
<br>
</div>
<div>
$ cat lcc-lcca.gie
</div>
<div>
<br>
</div>
<div>
<div>
BEGIN
</div>
<div>
<br>
</div>
<div>
------------------------------<wbr>------------------------------<wbr>-------------------
</div>
<div>
operation +proj=lcc +lat_1=48.73333333333333 +lat_2=47.5 \
</div>
<div>
+lat_0=47 +lon_0=-120.8333333333333 \
</div>
<div>
+x_0=500000.0001016001 +y_0=0 \
</div>
<div>
+units=us-ft +ellps=GRS80 +towgs84=0,0,0 +no_defs
</div>
<div>
------------------------------<wbr>------------------------------<wbr>-------------------
</div>
<div>
tolerance 0.0010 mm
</div>
<div>
accept -118.5293900000000300 48.7408309999860520
</div>
<div>
roundtrip 1000
</div>
<div>
<br>
</div>
<div>
<br>
</div>
<div>
------------------------------<wbr>------------------------------<wbr>-------------------
</div>
<div>
operation +proj=lcca +lat_1=48.73333333333333 +lat_2=47.5 \
</div>
<div>
+lat_0=47 +lon_0=-120.8333333333333 \
</div>
<div>
+x_0=500000.0001016001 +y_0=0 \
</div>
<div>
+units=us-ft +ellps=GRS80 +towgs84=0,0,0 +no_defs
</div>
<div>
------------------------------<wbr>------------------------------<wbr>-------------------
</div>
<div>
tolerance 0.0 mm
</div>
<div>
accept -118.5293900000000300 48.7408309999860520
</div>
<div>
roundtrip 1000
</div>
<div>
END
</div>
</div>
<div>
<br>
</div>
<div>
I get this output:
</div>
<div>
<br>
</div>
<div>
$ gie lcc-lcca.gie
</div>
<div>
<div>
<br>
</div>
<div>
------------------------------<wbr>------------------------------<wbr>-------------------
</div>
<div>
Reading file '..\..\..\proj\lcc-lcca.gie'
</div>
<div>
------------------------------<wbr>------------------------------<wbr>-------------------
</div>
<div>
+proj=lcc +lat_1=48.73333333333333 +lat_2=47.5 +lat_0=47 +lon_0=-120....
</div>
<div>
------------------------------<wbr>------------------------------<wbr>-------------------
</div>
<div>
FAILURE in lcc-lcca.gie(11):
</div>
<div>
roundtrip deviation: 1.550 mm, expected: 0.001 mm
</div>
<div>
------------------------------<wbr>------------------------------<wbr>-------------------
</div>
<div>
total: 1 tests succeeded, 1 tests FAILED!
</div>
<div>
------------------------------<wbr>------------------------------<wbr>-------------------
</div>
</div>
<div>
<br>
</div>
<div>
i.e. lcca roundtrips exactly, while lcc diverges heavily.
</div>
<div>
<br>
</div>
<div>
So Dmitry: This seems to be your workaround - define your projection using lcca,
</div>
<div>
rather than letting the epsg-list select lcc for you.
</div>
<div>
<br>
</div>
<div>
/Thomas
</div>
<div>
<br>
</div><div><br></div><div><br></div><div><br></div></div><br></div></div></blockquote><br><p></p></div><br>______________________________<wbr>_________________<br>
Proj mailing list<br>
<a href="mailto:Proj@lists.maptools.org" target="_blank">Proj@lists.maptools.org</a><br>
<a href="http://lists.maptools.org/mailman/listinfo/proj" rel="noreferrer" target="_blank">http://lists.maptools.org/mail<wbr>man/listinfo/proj</a><br></blockquote></div><br></div>
<br>______________________________<wbr>_________________<br>
Proj mailing list<br>
<a href="mailto:Proj@lists.maptools.org" target="_blank">Proj@lists.maptools.org</a><br>
<a href="http://lists.maptools.org/mailman/listinfo/proj" rel="noreferrer" target="_blank">http://lists.maptools.org/mail<wbr>man/listinfo/proj</a><br></blockquote></div><br></div>
</div></div><br>______________________________<wbr>_________________<br>
Proj mailing list<br>
<a href="mailto:Proj@lists.maptools.org" target="_blank">Proj@lists.maptools.org</a><br>
<a href="http://lists.maptools.org/mailman/listinfo/proj" rel="noreferrer" target="_blank">http://lists.maptools.org/mail<wbr>man/listinfo/proj</a><br></blockquote></div><br></div>
<br>______________________________<wbr>_________________<br>
Proj mailing list<br>
<a href="mailto:Proj@lists.maptools.org" target="_blank">Proj@lists.maptools.org</a><br>
<a href="http://lists.maptools.org/mailman/listinfo/proj" rel="noreferrer" target="_blank">http://lists.maptools.org/mail<wbr>man/listinfo/proj</a><br></blockquote></div><br></div>
</div></div><br>______________________________<wbr>_________________<br>
Proj mailing list<br>
<a href="mailto:Proj@lists.maptools.org" target="_blank">Proj@lists.maptools.org</a><br>
<a href="http://lists.maptools.org/mailman/listinfo/proj" rel="noreferrer" target="_blank">http://lists.maptools.org/mail<wbr>man/listinfo/proj</a><br></blockquote></div><br></div>
<br>______________________________<wbr>_________________<br>
Proj mailing list<br>
<a href="mailto:Proj@lists.maptools.org">Proj@lists.maptools.org</a><br>
<a href="http://lists.maptools.org/mailman/listinfo/proj" rel="noreferrer" target="_blank">http://lists.maptools.org/<wbr>mailman/listinfo/proj</a><br></blockquote></div><br></div>