[Proj] Grown error if re-projecting from 4269 to LCC (2285) and backward multiple times
Thomas Knudsen
knudsen.thomas at gmail.com
Wed Nov 1 01:58:16 EST 2017
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
2017-10-31 21:11 GMT+01:00 vanadovv at hetnet.nl <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
>
>
>
> _______________________________________________
> 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.maptools.org/pipermail/proj/attachments/20171101/fb551ac5/attachment-0001.htm
More information about the Proj
mailing list