[gdal-dev] GDAL 3.3.0beta1 available

Even Rouault even.rouault at spatialys.com
Sat Apr 24 08:38:19 PDT 2021


oh.... so I was wrong and you've spotted the reason for what Greg 
notices. So this is (was) a GDAL 3.3 specific bug (3.2 is fine).

GDAL should have passed HUGE_VAL when no explicit time is specified, 
which would be interpreted by PROJ as ignoring the time-dependent parts 
(or equivalently using the reference epoch of the transformation, which 
nullifies the effect of those time-dependent parts), but it incorrectly 
passes 0 currently

Transformation at epoch 0:

$ echo 100 -50 0 0 | gdaltransform -s_srs EPSG:6319 -t_srs EPSG:7912
100.000520678775 -49.9999650963593 1.20454746671021

==> ~ 100 m shift

Transformation at epoch 2010;

$ echo 100 -50 0 2010  | gdaltransform -s_srs EPSG:6319 -t_srs EPSG:7912
100.00001055598 -49.9999754213041 0.900845746509731

Current incorrect behavior is the same as at epoch 0

$ echo 100 -50 | gdaltransform -s_srs EPSG:6319 -t_srs EPSG:7912
100.000520678775 -49.9999650963593 1.20454746671021

With fix queued in https://github.com/OSGeo/gdal/pull/3738

we now get the same as epoch 2010

$ echo 100 -50 | gdaltransform -s_srs EPSG:6319 -t_srs EPSG:7912
100.00001055598 -49.9999754213041 0.900845746509731

Even

Le 24/04/2021 à 17:13, Javier Jimenez Shaw a écrit :
> Just one idea: which time is setting GDAL (and/or PROJ understanding) 
> on the transformation if nothing is specified? Is it the same as in cs2cs?
> That transformation is time dependent.
> I read somewhere that the "default" value is not the same (but my 
> memory is not very reliable).
>
> Cheers.
> .___ ._ ..._ .. . ._. .___ .. __ . _. . __..  ... .... ._ .__
> Entre dos pensamientos racionales
> hay infinitos pensamientos irracionales.
>
>
>
> On Sat, 24 Apr 2021 at 16:58, Even Rouault <even.rouault at spatialys.com 
> <mailto:even.rouault at spatialys.com>> wrote:
>
>     Greg,
>
>     There's a high chance for this to be completely PROJ behavior. GDAL
>     hardly knows about CRS those days :-)
>
>     Could you share your input and output files ?
>
>     Also set CPL_DEBUG=PROJ and PROJ_DEBUG=2 as environment variables
>     to see
>     which transformations PROJ pickup : look for lines starting with
>     "PROJ:
>     Using coordinate operation ..."
>
>     With latest PROJ master (and also 6.3.2), I see that
>
>     projinfo -s EPSG:6319 -t EPSG:7912 -o PROJ
>
>     Candidate operations found: 1
>     -------------------------------------
>     Operation No. 1:
>
>     unknown id, Conversion from NAD83(2011) (geog3D) to NAD83(2011)
>     (geocentric) + Inverse of ITRF2014 to NAD83(2011) (1) + Conversion
>     from
>     ITRF2014 (geocentric) to ITRF2014 (geog3D), 0 m, Puerto Rico -
>     onshore
>     and offshore. United States (USA) onshore and offshore - Alabama;
>     Alaska; Arizona; Arkansas; California; Colorado; Connecticut;
>     Delaware;
>     Florida; Georgia; Idaho; Illinois; Indiana; Iowa; Kansas; Kentucky;
>     Louisiana; Maine; Maryland; Massachusetts; Michigan; Minnesota;
>     Mississippi; Missouri; Montana; Nebraska; Nevada; New Hampshire; New
>     Jersey; New Mexico; New York; North Carolina; North Dakota; Ohio;
>     Oklahoma; Oregon; Pennsylvania; Rhode Island; South Carolina; South
>     Dakota; Tennessee; Texas; Utah; Vermont; Virginia; Washington; West
>     Virginia; Wisconsin; Wyoming. US Virgin Islands - onshore and
>     offshore.
>
>     PROJ string:
>     +proj=pipeline
>        +step +proj=axisswap +order=2,1
>        +step +proj=unitconvert +xy_in=deg +z_in=m +xy_out=rad +z_out=m
>        +step +proj=cart +ellps=GRS80
>        +step +inv +proj=helmert +x=1.0053 +y=-1.9092 +z=-0.5416
>     +rx=0.0267814
>              +ry=-0.0004203 +rz=0.0109321 +s=0.00037 +dx=0.0008
>     +dy=-0.0006
>              +dz=-0.0014 +drx=6.67e-05 +dry=-0.0007574 +drz=-5.13e-05
>     +ds=-7e-05
>              +t_epoch=2010 +convention=coordinate_frame
>        +step +inv +proj=cart +ellps=GRS80
>        +step +proj=unitconvert +xy_in=rad +z_in=m +xy_out=deg +z_out=m
>        +step +proj=axisswap +order=2,1
>
>     So it applies a time-dependent Helmert transformation, which has also
>     non-time dependent terms.
>
>     However trying a bit, the transformation seems to introduce a
>     shift of
>     about ~1 m or so, not 35 m
>
>     $ echo 50 -100 0 | cs2cs -d 8 EPSG:6319 EPSG:7912
>     50.00000754    -100.00001400 -0.68652183
>
>     $ echo 50 -100 50.00000754 -100.00001400 | geod -I +ellps=GRS80
>     -50d7'11.377"    129d52'48.584"    1.308
>
>     Even
>
>     Le 24/04/2021 à 16:39, Greg Troxel a écrit :
>     > Even Rouault <even.rouault at spatialys.com
>     <mailto:even.rouault at spatialys.com>> writes:
>     >
>     >> - http://download.osgeo.org/gdal/3.3.0/gdal-3.3.0beta1.tar.gz
>     <http://download.osgeo.org/gdal/3.3.0/gdal-3.3.0beta1.tar.gz>
>     > The tl;dr is:
>     >
>     >    - I'm seeing bad transforms NAD83/ITRF with 3.3.0beta1 but I
>     do not
>     >      have any basis to blame 3.3.0 because a lot of other things
>     are in
>     >      flux.
>     >
>     >    - If you are using 3.3.0beta1 and transforms I suggest
>     looking at your
>     >      output data.
>     >
>     >
>     >
>     > I have noticed a problem but I have little idea where it is.  I
>     will be
>     > looking into it myself but I wanted to mention it so that others
>     could
>     > be on the lookout.
>     >
>     > Besides updating to gdal, I'm doing a general software
>     freshening and
>     > rebuilding lots of things, but aside from the gdal beta my geo
>     packages
>     > haven't changed versions.  I have proj 6.3.2 with installed
>     grids, which
>     > was just rebuilt.
>     >
>     > The background is that I'm mostly working in EPSG:6319 which is
>     > NAD83(2011).
>     >
>     > I transform to EPSG:7912 which is ITRF2014 to display the data on
>     > leaflet webmaps (as a proxy for WGS84 to avoid null transforms).
>     >
>     > I noticed that my points on the webmap were very far off,
>     displaced very
>     > roughly 35m ESE.  I then put an NAD83 version in the webspace,
>     and that
>     > made things look ok (with the ~1 meter error from treating NAD83 and
>     > ITRF2014 as equal, surely).
>     >
>     > I'm extracting things from a geopackage in 6319 like
>     >
>     >      ogr2ogr -f GeoJSON -t_srs "EPSG:7912" -s_srs "EPSG:6319"
>     FOO/issues.geojson foo.gpkg issues
>     >      ogr2ogr -f GeoJSON -t_srs "EPSG:6319" -s_srs "EPSG:6319"
>     FOO/issues-nad83.geojson foo.gpkg issues
>     >
>     > which has previously been fine.  I'm definitely getting very large
>     > offsets from two files extracted seconds apart from the very same
>     > geopackage.
>     >
>     > So this is not a request for help, but if you are using ogr2ogr to
>     > transform datums, and you are using gdal 3.3.0beta1 especially
>     with proj
>     > 6, I'd suggest a quick look at your output data to see if it
>     seems ok.
>     >
>     > I will be retesting after all the software rebuilds are done, and
>     > rolling back to released gdal.  (Yes, I know proj 6.3 is old and
>     that I
>     > should update.)
>     >
>     > Thanks,
>     > Greg
>
>     -- 
>     http://www.spatialys.com <http://www.spatialys.com>
>     My software is free, but my time generally not.
>
>     _______________________________________________
>     gdal-dev mailing list
>     gdal-dev at lists.osgeo.org <mailto:gdal-dev at lists.osgeo.org>
>     https://lists.osgeo.org/mailman/listinfo/gdal-dev
>     <https://lists.osgeo.org/mailman/listinfo/gdal-dev>
>
-- 
http://www.spatialys.com
My software is free, but my time generally not.

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20210424/f2076a21/attachment-0001.html>


More information about the gdal-dev mailing list