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

[vanadovv at hetnet.nl]

Hi,
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.
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.
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.
Greetings,
Oscar van Vlijmen
----Origineel Bericht----
Van : dmitrymefed at gmail.com
Datum : 31/10/2017 07:50
Aan : proj at lists.maptools.org
Onderwerp : Re: [Proj] Grown error if re-projecting from 4269 to LCC (2285) and backward multiple times
 
Hi,
 
 
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 
 
1.       4269 -> 2285 projection, 
 
2.       then User modifies the entity (in 2285 CRS), 
 
3.       In order to store the changes, we do 2285 -> 4269 projection 
 
4.       and next time the entity is requested we do 4269 -> 2285 projection again.
 
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):
 
C:\PROJ\bin>cs2cs.exe -v +init=epsg:2285 +to +init=epsg:4269 -f %.16f
 
 
pj_open_lib(epsg): call fopen(C:\proj\share\epsg) - succeeded
 
 
 
 
 
pj_open_lib(epsg): call fopen(C:\proj\share\epsg) - succeeded
 
 
 
 
 
# ---- From Coordinate System ----
 
 
#Lambert Conformal Conic
 
 
#       Conic, Sph&Ell
 
 
#       lat_1= and lat_2= or lat_0
 
 
# +init=epsg:2285 +proj=lcc +lat_1=48.73333333333333 +lat_2=47.5 +lat_0=47
 
 
# +lon_0=-120.8333333333333 +x_0=500000.0001016001 +y_0=0 +datum=NAD83
 
 
# +units=us-ft +no_defs +ellps=GRS80 +towgs84=0,0,0
 
 
# ---- To Coordinate System ----
 
 
#Lat/long (Geodetic alias)
 
 
#
 
 
# +init=epsg:4269 +proj=longlat +datum=NAD83 +no_defs +ellps=GRS80
 
 
# +towgs84=0,0,0
 
 
2196293.3184066000 643350.39301622100
 
 
-118.5293900000000300   48.7408309999860520 0.0000000000000000
 
 
2196293.3184067495000000      
       643350.3930111383100000
 
 
-118.5293900000000300   48.7408309999721180 0.0000000000000000
 
 
2196293.3184069018000000      
       643350.3930060574800000
 
 
-118.5293900000000300   48.7408309999581850 0.0000000000000000
 
 
2196293.3184070541000000      
       643350.3930009742000000
 
 
-118.5293900000000300   48.7408309999442440 0.0000000000000000
 
 
2196293.3184072063000000      
       643350.3929958911600000
 
 
-118.5293900000000300   48.7408309999302960 0.0000000000000000
 
 
2196293.3184073586000000      
       643350.3929908055600000
 
 
-118.5293900000000300   48.7408309999163550 0.0000000000000000
 
 
2196293.3184075109000000      
       643350.3929857177400000
 
 
-118.5293900000000300   48.7408309999023930 0.0000000000000000
 
 
2196293.3184076636000000      
       643350.3929806276000000
 
 
-118.5293900000000300   48.7408309998884450 0.0000000000000000
 
 
ETC. ETC.
 
Thanks,
 
-Dmitry Mefed
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.maptools.org/pipermail/proj/attachments/20171031/150e088d/attachment.htm