[PROJ] PROJ 6.0.0 cs2cs: cannot instantiate source coordinate system - program abnormally terminated
Even Rouault
even.rouault at spatialys.com
Thu Mar 21 09:48:43 PDT 2019
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
More information about the PROJ
mailing list