5

I'm currently working on a better support of LAS 1.4 file format in the rlas R package. In LAS 1.4 format the CRS of the point cloud can be stored as a WKT string. For example this is what I read from a LAS 1.4 file (raw data).

"COMPD_CS[\"Projected\", PROJCS[\"UTM_10N\", GEOGCS [ \"WGS84\", DATUM [ \"WGS84\", SPHEROID [\"WGS 84\", 6378137.000, 298.257223563 ], TOWGS84 [ 0.000, 0.000, 0.000, 0.0000000000, 0.0000000000, 0.0000000000, 0.0000000000 ] ], PRIMEM [ \"Greenwich\", 0.000000 ], UNIT [ \"metres\", 1.00000000] ], PROJECTION[\"Transverse_Mercator\"], PARAMETER[\"Latitude_of_Origin\",0.0000000000], PARAMETER[\"Central_Meridian\",-123.0000000000], PARAMETER[\"Scale_Factor\",0.9996000000], PARAMETER[\"False_Easting\",500000.000], PARAMETER[\"False_Northing\",0.000], UNIT [ \"metres\", 1.00000000]] ], VERT_CS[\"NAVD88 (Geoid03) ContUS\", VERT_DATUM[\"./Resources/CoordSysData/navd88_geo03_contus.bin\", 1 ], UNIT [ \"metres\", 1.00000000] ] ]"

How can I interpret this string in R? In the R spatial ecosystem all packages rely on PROJ4. How I can convert this WKT string into PROJ4 string?

geos::readWKT does not understand this string.

  • There's some python solutions here that don't require lookups of any data base and work by converting the string format: gis.stackexchange.com/questions/8547/… - the osgeo lib works on your example. If you can call python from R then that might be a solution? – Spacedman Jan 5 at 16:38
5

Use the rgdal package and showP4:

Your string:

> ps
[1] "COMPD_CS[\"Projected\", PROJCS[\"UTM_10N\", GEOGCS [ \"WGS84\", DATUM [ \"WGS84\", SPHEROID [\"WGS 84\", 6378137.000, 298.257223563 ], TOWGS84 [ 0.000, 0.000, 0.000, 0.0000000000, 0.0000000000, 0.0000000000, 0.0000000000 ] ], PRIMEM [ \"Greenwich\", 0.000000 ], UNIT [ \"metres\", 1.00000000] ], PROJECTION[\"Transverse_Mercator\"], PARAMETER[\"Latitude_of_Origin\",0.0000000000], PARAMETER[\"Central_Meridian\",-123.0000000000], PARAMETER[\"Scale_Factor\",0.9996000000], PARAMETER[\"False_Easting\",500000.000], PARAMETER[\"False_Northing\",0.000], UNIT [ \"metres\", 1.00000000]] ], VERT_CS[\"NAVD88 (Geoid03) ContUS\", VERT_DATUM[\"./Resources/CoordSysData/navd88_geo03_contus.bin\", 1 ], UNIT [ \"metres\", 1.00000000] ] ]"

Your string as a Proj4 string:

> showP4(ps)
[1] "+proj=tmerc +lat_0=0 +lon_0=-7047.380880109124 +k=0.9996 +x_0=500000 +y_0=0 +ellps=WGS84 +towgs84=0.000,0.000,0.000,0.0000000000,0.0000000000,0.0000000000,0.0000000000 +units=m +no_defs "

Or using the sf package with st_crs and the wkt= parameter:

> st_crs(wkt=ps)
Coordinate Reference System:
  No EPSG code
  proj4string: "+proj=tmerc +lat_0=0 +lon_0=-7047.380880109124 +k=0.9996 +x_0=500000 +y_0=0 +ellps=WGS84 +towgs84=0.000,0.000,0.000,0.0000000000,0.0000000000,0.0000000000,0.0000000000 +units=m +no_defs"

Note I'm not convinced every parameter is always handled or converted properly... Caveat emptor.

  • I searched in rgdal but showP4 is not a trivial name. Thanks. That should make the job. – JRR Jan 5 at 17:34
  • Its not obvious - the sf package has an improved interface and names. ??wkt finds showWKT as "Show Well-Known Text spatial reference system metadata" and the help for that lists showP4 as well. Think thats the roundabout route that led me to it. – Spacedman Jan 5 at 17:41
  • There is a problem in the WKT definition, at GEOGCS [ ... , UNIT [ "metres" , 1.00000000 ] ], .... The units of the GEOGCS must be ["degree", 0.01745329251994328]. Thats why the central meridian is converted to +lon_0=-7047.380880109124. – Gabriel De Luca Feb 23 at 0:10
0

Since the WKT tells you about WGS84 UTM 10N, try http://epsg.io/32610 as a start-off.

For the vertical datum, see https://github.com/OSGeo/proj.4/wiki/VerticalDatums

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.