<div dir="ltr"><div><div><div><div><div>Thank you all for the interesting discussion.<br><br></div>I&#39;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="gmail_extra"><br><div class="gmail_quote">On Wed, Nov 1, 2017 at 4:33 PM, Thomas Knudsen <span dir="ltr">&lt;<a href="mailto:knudsen.thomas@gmail.com" target="_blank">knudsen.thomas@gmail.com</a>&gt;</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">&lt;<a href="mailto:vanadovv@hetnet.nl" target="_blank">vanadovv@hetnet.nl</a>&gt;</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_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 &quot;gie&quot; 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 &#39;..\..\..\proj\lcc-lcca.gie&#39;
  </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">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>