7

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.

Andre Silva
  • 10,259
  • 12
  • 54
  • 106
JRR
  • 9,389
  • 1
  • 13
  • 28
  • There's some python solutions here that don't require lookups of any data base and work by converting the string format: https://gis.stackexchange.com/questions/8547/what-is-the-best-way-to-programmatically-convert-between-wkt-and-proj4-string - the osgeo lib works on your example. If you can call python from R then that might be a solution? – Spacedman Jan 05 '19 at 16:38

2 Answers2

8

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.

Spacedman
  • 63,755
  • 5
  • 81
  • 115
  • I searched in rgdal but showP4 is not a trivial name. Thanks. That should make the job. – JRR Jan 05 '19 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 05 '19 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 '19 at 00: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

AndreJ
  • 76,698
  • 5
  • 86
  • 162