<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">https://github.com/OSGeo/proj.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>-------------------------------------------------------------------------------</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>-------------------------------------------------------------------------------</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>-------------------------------------------------------------------------------</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>-------------------------------------------------------------------------------</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>-------------------------------------------------------------------------------</div><div>Reading file '..\..\..\proj\lcc-lcca.gie'</div><div>-------------------------------------------------------------------------------</div><div>+proj=lcc +lat_1=48.73333333333333 +lat_2=47.5 +lat_0=47 +lon_0=-120....</div><div>-------------------------------------------------------------------------------</div><div> FAILURE in lcc-lcca.gie(11):</div><div> roundtrip deviation: 1.550 mm, expected: 0.001 mm</div><div>-------------------------------------------------------------------------------</div><div>total: 1 tests succeeded, 1 tests FAILED!</div><div>-------------------------------------------------------------------------------</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><div class="gmail_extra"><br><div class="gmail_quote">2017-10-31 21:11 GMT+01:00 <a href="mailto:vanadovv@hetnet.nl">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"><div>Hi,</div><div><br></div><div>With my independent software, using the EPSG Guidance Note 7-2 formulae, I arrive for this problem at a difference of 9.3e-10 m in the easting and 1.7e-9 m in the northing for the first round-trip calculation.</div><div><br></div><div>Thirty years ago I would expect that with an accuracy of merely around 1e-7 someone could have calculated things with single instead of double precision.<br>Nowadays I could expect fewer iteration steps or a looser iteration stopping criterion, or accuracy problems with transcendental functions (sin, cos, atan, log etc.), which are not part of PROJ.</div><div><br></div><div>Greetings,</div><div>Oscar van Vlijmen</div><div><br></div><br><blockquote style="margin-right:0px;margin-left:15px">----Origineel Bericht----<br>Van : <a href="mailto:dmitrymefed@gmail.com" target="_blank">dmitrymefed@gmail.com</a><br>Datum : 31/10/2017 07:50<br>Aan : <a href="mailto:proj@lists.maptools.org" target="_blank">proj@lists.maptools.org</a><br>Onderwerp : Re: [Proj] Grown error if re-projecting from 4269 to LCC (2285) and backward multiple times<br><br><div dir="ltr"><div><div class="h5">
<p class="MsoNormal"><span>Hi,</span></p>
<p class="MsoNormal"><span><br></span></p>
<p class="MsoNormal"><span>Our software stores the data with EPSG 4269 and present to the User in a local projected CRS, EPSG 2285 in examples described below. In case if an entity is modified on the UI we need to push the changes back. In that scenario we have to do </span></p>
<p class="m_5500118343644927306gmail-m_-5491966359676502627gmail-MsoListParagraphCxSpFirst"><span><span>1.<span style="font:7.0pt "Times New Roman""> </span></span></span><span>4269 -> 2285 projection, </span></p>
<p class="m_5500118343644927306gmail-m_-5491966359676502627gmail-MsoListParagraphCxSpMiddle"><span><span>2.<span style="font:7.0pt "Times New Roman""> </span></span></span><span>then User modifies the entity (in 2285 CRS), </span></p>
<p class="m_5500118343644927306gmail-m_-5491966359676502627gmail-MsoListParagraphCxSpMiddle"><span><span>3.<span style="font:7.0pt "Times New Roman""> </span></span></span><span>In order to store the changes, we do 2285 -> 4269 projection </span></p>
<p class="m_5500118343644927306gmail-m_-5491966359676502627gmail-MsoListParagraphCxSpLast"><span><span>4.<span style="font:7.0pt "Times New Roman""> </span></span></span><span>and next time the entity is requested we do 4269 -> 2285 projection again.</span></p>
<p class="MsoNormal"><span>The problem is that proj4 does not produce same results in step 2 and 4. It would be expected behavior as 2285 provides more precision due to used units US-ft, in comparison to 4269 which units are degrees. But the real issue for us is that subsequent transformations between 2285 -> 4269 -> 2285… produces growing error. Please consider the following output of cs2cs program, (in parallel I did inverse projection to obtain coordinates in EPSG 4269):</span></p>
<p style="line-height:normal" class="MsoNormal"><font size="1"><span style="font-family:monospace,monospace">C:\PROJ\bin>cs2cs.exe -v +init=epsg:2285 +to +init=epsg:4269 -f %.16f</span></font></p>
<font size="1"><span style="font-family:monospace,monospace"></span></font>
<p style="line-height:normal" class="MsoNormal"><font size="1"><span style="font-family:monospace,monospace">pj_open_lib(epsg): call fopen(C:\proj\share\epsg) - succeeded</span></font></p>
<font size="1"><span style="font-family:monospace,monospace"></span></font>
<p style="line-height:normal" class="MsoNormal"><font size="1"><span style="font-family:monospace,monospace"> </span></font></p>
<font size="1"><span style="font-family:monospace,monospace"></span></font>
<p style="line-height:normal" class="MsoNormal"><font size="1"><span style="font-family:monospace,monospace">pj_open_lib(epsg): call fopen(C:\proj\share\epsg) - succeeded</span></font></p>
<font size="1"><span style="font-family:monospace,monospace"></span></font>
<p style="line-height:normal" class="MsoNormal"><font size="1"><span style="font-family:monospace,monospace"> </span></font></p>
<font size="1"><span style="font-family:monospace,monospace"></span></font>
<p style="line-height:normal" class="MsoNormal"><font size="1"><span style="font-family:monospace,monospace"># ---- From Coordinate System ----</span></font></p>
<font size="1"><span style="font-family:monospace,monospace"></span></font>
<p style="line-height:normal" class="MsoNormal"><font size="1"><span style="font-family:monospace,monospace">#Lambert Conformal Conic</span></font></p>
<font size="1"><span style="font-family:monospace,monospace"></span></font>
<p style="line-height:normal" class="MsoNormal"><font size="1"><span style="font-family:monospace,monospace">#<span> </span>Conic, Sph&Ell</span></font></p>
<font size="1"><span style="font-family:monospace,monospace"></span></font>
<p style="line-height:normal" class="MsoNormal"><font size="1"><span style="font-family:monospace,monospace">#<span> </span>lat_1= and lat_2= or lat_0</span></font></p>
<font size="1"><span style="font-family:monospace,monospace"></span></font>
<p style="line-height:normal" class="MsoNormal"><font size="1"><span style="font-family:monospace,monospace"># +init=epsg:2285 +proj=lcc +lat_1=48.73333333333333 +lat_2=47.5 +lat_0=47</span></font></p>
<font size="1"><span style="font-family:monospace,monospace"></span></font>
<p style="line-height:normal" class="MsoNormal"><font size="1"><span style="font-family:monospace,monospace"># +lon_0=-120.8333333333333 +x_0=500000.0001016001 +y_0=0 +datum=NAD83</span></font></p>
<font size="1"><span style="font-family:monospace,monospace"></span></font>
<p style="line-height:normal" class="MsoNormal"><font size="1"><span style="font-family:monospace,monospace"># +units=us-ft +no_defs +ellps=GRS80 +towgs84=0,0,0</span></font></p>
<font size="1"><span style="font-family:monospace,monospace"></span></font>
<p style="line-height:normal" class="MsoNormal"><font size="1"><span style="font-family:monospace,monospace"># ---- To Coordinate System ----</span></font></p>
<font size="1"><span style="font-family:monospace,monospace"></span></font>
<p style="line-height:normal" class="MsoNormal"><font size="1"><span style="font-family:monospace,monospace">#Lat/long (Geodetic alias)</span></font></p>
<font size="1"><span style="font-family:monospace,monospace"></span></font>
<p style="line-height:normal" class="MsoNormal"><font size="1"><span style="font-family:monospace,monospace">#</span></font></p>
<font size="1"><span style="font-family:monospace,monospace"></span></font>
<p style="line-height:normal" class="MsoNormal"><font size="1"><span style="font-family:monospace,monospace"># +init=epsg:4269 +proj=longlat +datum=NAD83 +no_defs +ellps=GRS80</span></font></p>
<font size="1"><span style="font-family:monospace,monospace"></span></font>
<p style="line-height:normal" class="MsoNormal"><font size="1"><span style="font-family:monospace,monospace"># +towgs84=0,0,0</span></font></p>
<font size="1"><span style="font-family:monospace,monospace"></span></font>
<p style="line-height:normal" class="MsoNormal"><font size="1"><span style="font-family:monospace,monospace">2196293.3184066000 643350.39301622100</span></font></p>
<font size="1"><span style="font-family:monospace,monospace"></span></font>
<p style="line-height:normal" class="MsoNormal"><font size="1"><span style="font-family:monospace,monospace">-118.5293900000000300<span> </span>48.7408309999860520 0.0000000000000000</span></font></p>
<font size="1"><span style="font-family:monospace,monospace"></span></font>
<p style="line-height:normal" class="MsoNormal"><font size="1"><span style="font-family:monospace,monospace">2196293.3184067495000000<span>
<u></u><u></u> </span>643350.3930111383100000</span></font></p>
<font size="1"><span style="font-family:monospace,monospace"></span></font>
<p style="line-height:normal" class="MsoNormal"><font size="1"><span style="font-family:monospace,monospace">-118.5293900000000300<span> </span>48.7408309999721180 0.0000000000000000</span></font></p>
<font size="1"><span style="font-family:monospace,monospace"></span></font>
<p style="line-height:normal" class="MsoNormal"><font size="1"><span style="font-family:monospace,monospace">2196293.3184069018000000<span>
<u></u><u></u> </span>643350.3930060574800000</span></font></p>
<font size="1"><span style="font-family:monospace,monospace"></span></font>
<p style="line-height:normal" class="MsoNormal"><font size="1"><span style="font-family:monospace,monospace">-118.5293900000000300<span> </span>48.7408309999581850 0.0000000000000000</span></font></p>
<font size="1"><span style="font-family:monospace,monospace"></span></font>
<p style="line-height:normal" class="MsoNormal"><font size="1"><span style="font-family:monospace,monospace">2196293.3184070541000000<span>
<u></u><u></u> </span>643350.3930009742000000</span></font></p>
<font size="1"><span style="font-family:monospace,monospace"></span></font>
<p style="line-height:normal" class="MsoNormal"><font size="1"><span style="font-family:monospace,monospace">-118.5293900000000300<span> </span>48.7408309999442440 0.0000000000000000</span></font></p>
<font size="1"><span style="font-family:monospace,monospace"></span></font>
<p style="line-height:normal" class="MsoNormal"><font size="1"><span style="font-family:monospace,monospace">2196293.3184072063000000<span>
<u></u><u></u> </span>643350.3929958911600000</span></font></p>
<font size="1"><span style="font-family:monospace,monospace"></span></font>
<p style="line-height:normal" class="MsoNormal"><font size="1"><span style="font-family:monospace,monospace">-118.5293900000000300<span> </span>48.7408309999302960 0.0000000000000000</span></font></p>
<font size="1"><span style="font-family:monospace,monospace"></span></font>
<p style="line-height:normal" class="MsoNormal"><font size="1"><span style="font-family:monospace,monospace">2196293.3184073586000000<span>
<u></u><u></u> </span>643350.3929908055600000</span></font></p>
<font size="1"><span style="font-family:monospace,monospace"></span></font>
<p style="line-height:normal" class="MsoNormal"><font size="1"><span style="font-family:monospace,monospace">-118.5293900000000300<span> </span>48.7408309999163550 0.0000000000000000</span></font></p>
<font size="1"><span style="font-family:monospace,monospace"></span></font>
<p style="line-height:normal" class="MsoNormal"><font size="1"><span style="font-family:monospace,monospace">2196293.3184075109000000<span>
<u></u><u></u> </span>643350.3929857177400000</span></font></p>
<font size="1"><span style="font-family:monospace,monospace"></span></font>
<p style="line-height:normal" class="MsoNormal"><font size="1"><span style="font-family:monospace,monospace">-118.5293900000000300<span> </span>48.7408309999023930 0.0000000000000000</span></font></p>
<font size="1"><span style="font-family:monospace,monospace"></span></font>
<p style="line-height:normal" class="MsoNormal"><font size="1"><span style="font-family:monospace,monospace">2196293.3184076636000000<span>
<u></u><u></u> </span>643350.3929806276000000</span></font></p>
<font size="1"><span style="font-family:monospace,monospace"></span></font>
<p style="line-height:normal" class="MsoNormal"><font size="1"><span style="font-family:monospace,monospace">-118.5293900000000300<span> </span>48.7408309998884450 0.0000000000000000</span></font></p>
<p style="line-height:normal" class="MsoNormal"><span style="font-family:Consolas"> </span></p></div></div><p style="line-height:normal" class="MsoNormal"><span style="font-family:Consolas">ETC. ETC.</span></p><p style="line-height:normal" class="MsoNormal"><span style="font-family:Consolas"><br></span></p><p style="line-height:normal" class="MsoNormal"><br></p>
<p class="MsoNormal"><span>Thanks,</span></p>
<p class="MsoNormal"><span>-Dmitry Mefed</span></p>
</div><br></blockquote><br><p></p><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>