How to assign/project las-file using R?

Question

I have got a LAS file from a MLS, which has following specifications:

test <- lidR::readLAS("<file.las>")
print(test)
#> class        : LAS (v1.0 format 1)
#> memory       : 11.4 Gb 
#> extent       : -135.0345, 139.0671, -120.6662, 140.9799 (xmin, xmax, ymin, ymax)
#> coord. ref.  : NA 
#> area         : 55785.25 units²
#> points       : 152.99 million points
#> density      : 2742.47 points/units²

As you see, the LAS file does not come with any CRS. However, I have the coordinates for a specific point in the point cloud in both, the internal x/y/z dimensions as well as in the CRS I want to assign/project my dataset into.

Furthermore I know the units of the internal x/y/z dimensions as well as that the point cloud has got a northing.

How to assing/project the LAS file into a specific CRS, when the LAS file does not have a specific CRS yet? In my case EPSG:32632 see URL:https://epsg.io/32632 )

see also:

• How to assign (or project) CRS to LiDAR file in lidR?
• Assigning CRS using las2las from libLAS?

Answer 1

Your question is ill-defined or at least not clear enough. To answer I'll assume that you know the long/lat of the (0,0) in your point-cloud. Without this assumption the problem is not solvable. Also I'll assume that your point cloud is already projected with meters units.

First I'm generating a reproducible dataset

library(lidR)

X = runif(10000, -140, 180)
Y = runif(10000, -140, 180)
Z = runif(10000, 0, 20)
data = data.frame(X,Y,Z)
las = LAS(data)

Then I'm defining where on earth is (0,0) of the point cloud. Here I'm starting from long/lat but maybe you already have coordinates in your target CRS.

library(sf)

target = st_crs(32617)
zero = c(-79.6, 48.4)
zero = st_point(zero)
zero = st_sfc(zero, crs = 4326)
zero = st_transform(zero, target)

Now I'm extracting the X0 and Y0 that will be used to translate your coordinates

X0 <- st_coordinates(zero)[,1]
Y0 <- st_coordinates(zero)[,2]

And I'm applying the translation + assigning a CRS

las@header@PHB[["X offset"]] <- floor(X0)
las@header@PHB[["Y offset"]] <- floor(Y0)
las$X = las$X + X0
las$Y = las$Y + Y0
projection(las) <- target
las
#> class        : LAS (v1.2 format 0)
#> memory       : 248 Kb 
#> extent       : 603482.7, 603802.6, 5361567, 5361887 (xmin, xmax, ymin, ymax)
#> coord. ref.  : WGS 84 / UTM zone 17N 
#> area         : 102076.2 metre²
#> points       : 10 thousand points
#> density      : 0.1 points/metre²

---

For future reader: it is likely that future versions of lidR include a function translate_points() or similar name. Something like.

las <- translate_points(las, X0, Y0)