[QGIS-Developer] CRS problems with WGS84 and ETRS89

[Even Rouault]

> How come ? Or is the error-source sitting ca. 70 cm in front of the
> screen ?? (I have a 34" large flatscreen :-)

Because the data for such a transform isn't in the EPSG database.

I believe the main reason is that both WGS84 and ETRS89 are in fact datum ensembles, that is a short and convenient nickname for a number of more well-defined datum realizations. For WGS84, you have 6 different realizations from Transit to G1762. And ETRS89 has 11 realizations from ETRF89 to ETRF2014
So to get accurate data, you must research tranformations between those direct realizations. And in EPSG, accurate transformations tend not being given in terms of WGS84 realizations but from their equivalent ITRF realizations (WGS84(G1762) is identical to ITRF2014 within ~1cm),

so, while playing with PROJ' projinfo utility,

$ projinfo -s ITRF2014 -t ETRF2014

I see 2 potential coordinate transformations, and the one that makes probably more sense today is the one using epoch 2010 for the transformation:

unknown id, Conversion from ITRF2014 (geog2D) to ITRF2014 (geocentric) + ITRF2014 to ETRF2014 (2) + Conversion from ETRF2014 (geocentric) to ETRF2014 (geog2D), 0 m, Europe - ETRS89

PROJ string:
+proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=cart +ellps=GRS80 +step +proj=helmert +x=0 +y=0 +z=0 +rx=0.001785 +ry=0.011151 +rz=-0.01617 +s=0 +dx=0 +dy=0 +dz=0 +drx=8.5e-05 +dry=0.000531 +drz=-0.00077 +ds=0 +t_epoch=2010 +convention=position_vector +step +inv +proj=cart +ellps=GRS80 +step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1

And you need of course to specify the coordinate epoch of your coordinates as ITRF/WGS84 are Earth-fixed/dynamic datums.

echo "49 2 0 2019.6" | cct -d 10 +proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=cart +ellps=GRS80 +step +proj=helmert +x=0 +y=0 +z=0 +rx=0.001785 +ry=0.011151 +rz=-0.01617 +s=0 +dx=0 +dy=0 +dz=0 +drx=8.5e-05 +dry=0.000531 +drz=-0.00077 +ds=0 +t_epoch=2010 +convention=position_vector +step +inv +proj=cart +ellps=GRS80 +step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1

48.9999955186   1.9999924499  -0.0016582003     2019.6000

Actually using the first proposal using an epoch at 1989, I get exactly the same results

echo "49 2 0 2019.6" | cct -d 10 +proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=cart +ellps=GRS80 +step +proj=helmert +x=0 +y=0 +z=0 +rx=0 +ry=0 +rz=0 +s=0 +dx=0 +dy=0 +dz=0 +drx=8.5e-05 +dry=0.000531 +drz=-0.00077 +ds=0 +t_epoch=1989 +convention=position_vector +step +inv +proj=cart +ellps=GRS80 +step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1

48.9999955186   1.9999924499  -0.0016582003     2019.6000

That's a 74cm horizontal shift

-- 
Spatialys - Geospatial professional services
http://www.spatialys.com