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.

Improve this question
1
  • 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, 2019 at 16:38

2 Answers 2

Reset to default
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.

Improve this answer
3
  • I searched in rgdal but showP4 is not a trivial name. Thanks. That should make the job.
    – JRR
    Jan 5, 2019 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, 2019 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, 2019 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

Improve this answer

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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