Clean way to convert from US feet to meters with lidR - LiDAR point cloud

Question

I would like to convert this point cloud from US feet to meters. Metadata of the point cloud here.

I firstly though to convert it properly by re-projecting to a metric CRS but this caused me a lot of trouble. See the following posts (in their order of publication) :

But basically, what I want is just metric point clouds! So I came to the idea of simply converting the XYZ coordinates of the point cloud using las_rescale() from lidR R package:

library(lidR)
tmp_las <- readLAS('./USGS_LPC_CA_FEMA_R9_Russian_2017_LAS_2018.laz')
tmp_las_RS <- las_rescale(tmp_las,xscale = 0.3048 , yscale = 0.3048, zscale = 0.3048)
#> X coordinates were moved by 0.0762 on average
#> Y coordinates were moved by 0.0762 on average
#> Z coordinates were moved by 0.0762 on average

But the headers don’t take it into account

tmp_las_RS
#> class        : LAS (v1.2 format 3)
#> memory       : 618.7 Mb 
# extent       : 6268117, 6270217, 1998802, 2000902 (xmin, xmax, ymin, ymax)
#> coord. ref.  : +proj=lcc +lat_0=37.6666666666667 +lon_0=-122 +lat_1=39.8333333333333 +lat_2=38.3333333333333 +x_0=2000000.0001016 +y_0=500000.0001016 +ellps=GRS80 +units=us-ft +no_defs 
#> area         : 4.41 kus-ft²
#> points       : 7.05 million points
#> density      : 1.6 points/us-ft²

Is there a better way to do it through lidR?

coordinate system lidar lidr proj r

JRR · Accepted Answer

las_rescale() does not perform the task you think it does. You are not the first one to be confused. See this comprehensive answer: Rescaling and reoffsetting a point cloud with lidR?
Considering how a LAS file stores the data (header, scale factor, offset) it is actually difficult to perform what you want to do because you need to: (1) recompute the coordinates (2) quantize the coordinates (3) update the header (4) update the bounding box. All manually.
But in the version 3.1 things become easier (and better documented). See this thread. So from version 3.1 you can do
ft2m = 1200/3937
las$X = las$X * ft2m
las$Y = las$Y * ft2m
las$Z = las$Z * ft2m

But this does not update the CRS magically. So you must erase the CRS because, according to your previous posts, you don't know what CRS to input.
crs(las) <- sp::CRS()

To install lidR 3.1 which won't be released before Jan 2021.
remotes::install_github("Jean-Romain/lidR", ref = "devel")

Clean way to convert from US feet to meters with lidR - LiDAR point cloud

One Answer

Add your own answers!

Ask a Question