[PROJ] PROJ 6.0.0 cs2cs: cannot instantiate source coordinate system - program abnormally terminated
Lesparre, Jochem
Jochem.Lesparre at kadaster.nl
Thu Mar 28 08:22:55 PDT 2019
Hi Even,
Soon we will publish a new NTv2 and Vdatum grids and new 7 parameters for EPSG:: 7415 to EPSG:: 4937. So I wouldlike to know how to avoid that it is unclear to Proj. what this is realated to. But maybe I should first attend the Proj. tutorial at the EUREF symposium in May, because I do not yet understand the Proj.6 syntax that you list.
By the way, ETRS89 (e.g. EPSG:4258) is bound to ITRS en WGS84 is linked to ITRS too.
Regards, Jochem
-----Original Message-----
From: Even Rouault <even.rouault at spatialys.com>
Sent: Thursday, March 21, 2019 5:49 PM
To: Lesparre, Jochem <Jochem.Lesparre at kadaster.nl>
Cc: proj at lists.osgeo.org; Sebastiaan Couwenberg <sebastic at xs4all.nl>
Subject: Re: [PROJ] PROJ 6.0.0 cs2cs: cannot instantiate source coordinate system - program abnormally terminated
Jochem,
> What is unclear here? I still don't understand the early and late
> binding terminology so I can't understand your explanation.
Yeah, sorry. I'm not sure there's even a universally admitted definition for those. So in my mind,
- early-binding implies using WGS84 as a pivot when going from CRS A to CRS B.
This was the compulsory approach when using PROJ < 6. Actually if using geoidgrids/nadgrids in CRS A and towgs84=0,0,0 in CRS B, the pivot CRS would be something else than WGS84, but was a bit of cheating
- late-binding implies that if you want to go from CRS A to CRS B, the software will first research if there is a direct transformation between both, before researching tranformations through a potential pivot, such as WGS84.
This is the default behaviour of PROJ 6.
When transforming from EPSG:4258 to "+proj=sterea +lat_0=52.15616055555555
+lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000
++ellps=bessel nadgrids=rdtrans2008.gsb +geoidgrids=naptrans2008.gtx
++units=m +no_defs type=crs", we are in a middle situation.
EPSG:4258 is a CRS definition, that is not bound to a particular pivot
Whereas "+proj=sterea ... +nadgrids=.... +geoidgrids=...." implies normally a
WGS84 pivot, but as you mentionned here, in fact the nadgrids here was mean for going to the Netherlands' datum to ETRS89 directly, and PROJ has no way of figuring that out from the +proj= expression.
With PROJ master or latest state of 6.0 branch (not 6.0.0), if you expand
EPSG:4258 definition to "+proj=longlat +ellps=GRS80 +towgs84=0,0,0 +type=crs", you'll get the expected result:
$ src/projinfo -s "+proj=longlat +ellps=GRS80 +towgs84=0,0,0 +type=crs" -t "+proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079
+x_0=155000 +y_0=463000 +ellps=bessel +nadgrids=@rdtrans2008.gsb
+geoidgrids=@naptrans2008.gtx +units=m +no_defs +type=crs" -o PROJ
Candidate operations found: 1
-------------------------------------
Operation n°1:
unknown id, Transformation from unknown to WGS84 + Null geographic offset from WGS 84 to WGS 84 + unknown + unknown to WGS84 + unknown to WGS84 ellipsoidal height, unknown accuracy, World, at least one grid missing
PROJ string:
+proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +inv
+proj=vgridshift +grids=@naptrans2008.gtx +multiplier=1 +step +inv
+proj=hgridshift +grids=@rdtrans2008.gsb +step +proj=sterea
+lat_0=52.1561605555556 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000
+y_0=463000 +ellps=bessel
Grid @naptrans2008.gtx needed but not found on the system.
Grid @rdtrans2008.gsb needed but not found on the system.
~~~~~
Another approach here would be to use also a WKT definition from the +proj= part or a EPSG code, since the rdtrans2008.gsb grid is referenced in EPSG:
For example, if using EPSG:4289 (Amersfoort geographic CRS), you get:
$ src/projinfo -s EPSG:4258 -t EPSG:4289 -o PROJ Candidate operations found: 5
-------------------------------------
Operation n°1:
INVERSE(EPSG):15739, Inverse of Amersfoort to ETRS89 (3), 0.5 m, Netherlands - onshore
PROJ string:
+proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert
+xy_in=deg +xy_out=rad +step +proj=push +v_3 +step +proj=cart
++ellps=GRS80 step +inv +proj=helmert +x=565.2369 +y=50.0087 +z=465.658
+rx=0.406857330322398 +ry=-0.350732676542563 +rz=1.8703473836068
++s=4.0812 convention=coordinate_frame +step +inv +proj=cart
++ellps=bessel +step proj=pop +v_3 +step +proj=unitconvert +xy_in=rad
++xy_out=deg +step proj=axisswap +order=2,1
-------------------------------------
Operation n°2:
INVERSE(EPSG):15740, Inverse of Amersfoort to ETRS89 (4), 0.5 m, Netherlands - onshore
PROJ string:
+proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert
+xy_in=deg +xy_out=rad +step +proj=push +v_3 +step +proj=cart
++ellps=GRS80 step +inv +proj=molobadekas +x=593.0297 +y=26.0038
++z=478.7534
+rx=0.406857330322398 +ry=-0.350732676542563 +rz=1.8703473836068
++s=4.0812
+px=3903453.1482 +py=368135.3134 +pz=5012970.3051
++convention=coordinate_frame step +inv +proj=cart +ellps=bessel +step
++proj=pop +v_3 +step proj=unitconvert +xy_in=rad +xy_out=deg +step
++proj=axisswap +order=2,1
-------------------------------------
Operation n°3:
INVERSE(EPSG):4830, Inverse of Amersfoort to ETRS89 (5), 0.5 m, Netherlands - onshore
PROJ string:
+proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert
+xy_in=deg +xy_out=rad +step +proj=push +v_3 +step +proj=cart
++ellps=GRS80 step +inv +proj=helmert +x=565.4171 +y=50.3319 +z=465.5524
+rx=0.398957388243134 +ry=-0.343987817378283 +rz=1.87740163998045
++s=4.0725 convention=coordinate_frame +step +inv +proj=cart
++ellps=bessel +step proj=pop +v_3 +step +proj=unitconvert +xy_in=rad
++xy_out=deg +step proj=axisswap +order=2,1
-------------------------------------
Operation n°4:
INVERSE(EPSG):4831, Inverse of Amersfoort to ETRS89 (6), 0.5 m, Netherlands - onshore
PROJ string:
+proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert
+xy_in=deg +xy_out=rad +step +proj=push +v_3 +step +proj=cart
++ellps=GRS80 step +inv +proj=molobadekas +x=593.0248 +y=25.9984
++z=478.7459
+rx=0.398957388243134 +ry=-0.343987817378283 +rz=1.87740163998045
++s=4.0725
+px=3903453.1482 +py=368135.3134 +pz=5012970.3051
++convention=coordinate_frame step +inv +proj=cart +ellps=bessel +step
++proj=pop +v_3 +step proj=unitconvert +xy_in=rad +xy_out=deg +step
++proj=axisswap +order=2,1
-------------------------------------
Operation n°5:
INVERSE(EPSG):7000, Inverse of Amersfoort to ETRS89 (7), 0.001 m, Netherlands
- onshore, at least one grid missing
PROJ string:
+proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert
+xy_in=deg +xy_out=rad +step +inv +proj=hgridshift
++grids=rdtrans2008.gsb step +proj=unitconvert +xy_in=rad +xy_out=deg
++step +proj=axisswap +order=2,1
Grid rdtrans2008.gsb needed but not found on the system.
The vertical part is also there:
$ src/projinfo -s EPSG:5709 -t EPSG:4937 -o PROJ Candidate operations found: 1
-------------------------------------
Operation n°1:
INVERSE(EPSG):7001, Inverse of ETRS89 to NAP height (1), 0.01 m, Netherlands - onshore, at least one grid missing
PROJ string:
+proj=pipeline +step +inv +proj=vgridshift +grids=naptrans2008.gtx
+multiplier=1
Grid naptrans2008.gtx needed but not found on the system.
~~~~~~
So it is a matter of creating a CompoundCRS WKT with the ProjectedCRS using the Amersfoort geographic CRS and the vertical CRS using EPSG:5709, and asking to transform that to EPSG:4937
COMPOUNDCRS["unknown",
PROJCRS["unknown",
BASEGEOGCRS["Amersfoort",
DATUM["Amersfoort",
ELLIPSOID["Bessel 1841",6377397.155,299.1528128,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]]],
CONVERSION["unknown",
METHOD["Oblique Stereographic",
ID["EPSG",9809]],
PARAMETER["Latitude of natural origin",52.1561605555556,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",5.38763888888889,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",0.9999079,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",155000,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",463000,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["(E)",east,
ORDER[1],
LENGTHUNIT["metre",1,
ID["EPSG",9001]]],
AXIS["(N)",north,
ORDER[2],
LENGTHUNIT["metre",1,
ID["EPSG",9001]]]],
VERTCRS["NAP height",
VDATUM["Normaal Amsterdams Peil"],
CS[vertical,1],
AXIS["gravity-related height (H)",up,
LENGTHUNIT["metre",1]],
USAGE[
SCOPE["unknown"],
AREA["Netherlands - onshore"],
BBOX[50.75,3.2,53.7,7.22]],
ID["EPSG",5709]]]
The vertical part is correctly handled, but for the horizontal part, there's currently an issue to detect the Amersfoort datum and the associated grid transformation to ETRS89. I've just filed this as
https://github.com/OSGeo/proj.4/issues/1343
Even
--
Spatialys - Geospatial professional services http://www.spatialys.com
Disclaimer:
De inhoud van dit bericht is uitsluitend bestemd voor geadresseerde.
Gebruik van de inhoud van dit bericht door anderen zonder toestemming van het Kadaster
is onrechtmatig. Mocht dit bericht ten onrechte bij u terecht komen, dan verzoeken wij u
dit direct te melden aan de verzender en het bericht te vernietigen.
Aan de inhoud van dit bericht kunnen geen rechten worden ontleend.
Disclaimer:
The content of this message is meant to be received by the addressee only.
Use of the content of this message by anyone other than the addressee without the consent
of the Kadaster is unlawful. If you have received this message, but are not the addressee,
please contact the sender immediately and destroy the message.
No rights can be derived from the content of this message
More information about the PROJ
mailing list