How to efficiently reproject a LAS point cloud to EPSG:2056 (not known by LAStools)?

Question

I have a large LAS pointcloud in LV03 (EPSG:21781) that I want to reproject to LV95 (EPSG:2056). My approaches so far bear no fruit.

What is a good way to re-project LAS data to a different CRS (optionally efficient)?

---

I've tried so far:

• The excellent LAStools can't do this in my case, since it doesn't know EPSG:2056, which is a known issue. Normally las2las -epsg WXYZ -target_epsg ABCD would do the trick.

• The swisstopo Geosuite-Software (link, suggested here) produces a lot of errors for me. Results are only partially reprojected.

• I wanted to do in lidR, as suggested here with this but that results in the following error:

  LASfile <- system.file("extdata", "Megaplot.laz", package="lidR")
  las = readLAS(LASfile)
  las_transformed = st_transform(las, st_crs(2056))
  #> Error in UseMethod("st_transform"): not applicable method 'st_transform' applied to object of class "c('LAS', 'Spatial')"
  

with lidR 3.2.1

Answer 1

The only difference in those coordinate systems seems to be a change in the false origin value, which means all you really need to do is add the offset change to your coordinates. There's also another change to the datum but it appears to be some sort of administrative name change with no difference in the values. But we can test this empirically by creating some points in the first CRS, transforming to the second using st_transform, and seeing if all we have done is add a constant to the X and Y coordinates

Here's a test...

Create a grid of points over the projected bounds of 21781 as shown on epsg.io:

xy = expand.grid(x=seq(485072, 837118, len=50), y=seq(74261,299941,len=50))
pts21781 = st_as_sf(data.frame(xy), coords=1:2, crs=21781)

Check I've got this right by making sure these lat-long coords are as per the WGS84 bounds listed on epsg.io:

st_transform(pts21781, 4326)

OK good. Now transform to 2056 and compare value of original and projected coordinates:

pts2056 = st_transform(pts21781, 2056)
cc = st_coordinates(pts2056) - st_coordinates(pts21781)
range(cc[,1])
range(cc[,2])

This should show both the X and Y offset. How constant is this offset over the space? Let's look at the difference between the largest and smallest values in X and Y:

> min(cc[,1]) - max(cc[,1])
[1] -9.313226e-10
> min(cc[,2]) - max(cc[,2])
[1] -1.210719e-08

These values are in metres. Doing the full st_transform will differ from adding or subtracting a constant by 10^-8 metres - unless you care about precision down to a few hundred atoms, this looks good enough!