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 <[email protected] <mailto:[email protected]>> 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 <[email protected]
    <mailto:[email protected]>> 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
    [email protected] <mailto:[email protected]>
    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.

_______________________________________________
gdal-dev mailing list
[email protected]
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Reply via email to