search

r-lidar / lidR

Airborne LiDAR data manipulation and visualisation for forestry application
https://CRAN.R-project.org/package=lidR
GNU General Public License v3.0
582 stars 130 forks source link

Error when trying to merge pointclouds - Different CRS #725

Closed wiesehahn closed 10 months ago

wiesehahn commented 10 months ago

I try to merge two pointclouds of the same extent but run into problems related to the projection.

LAS1 has no CRS stored internally and is of class LAS (v1.2 format 1).

> st_crs(las1)
Coordinate Reference System: NA

LAS2 has a CRS stored and is of class LAS (v1.4 format 8).

> st_crs(las2)
Coordinate Reference System:
  User input: COMPD_CS["ETRS89 / ETRS89 / UTM 32N + DHHN2016",PROJCS["ETRS89 / ETRS89 / UTM 32N",GEOGCS["ETRS89",DATUM["European_Terrestrial_Reference_System_1989",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","6258"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4258"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",9],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","25832"]],VERT_CS["DHHN2016",VERT_DATUM["Deutsches Haupthoehennetz 2016",2005,AUTHORITY["EPSG","7837"]],UNIT["metre",1.0,AUTHORITY["EPSG","9001"]],AXIS["Gravity-related height",UP],AUTHORITY["EPSG","7837"]]] 
  wkt:
COMPOUNDCRS["ETRS89 / ETRS89 / UTM 32N + DHHN2016",
    PROJCRS["ETRS89 / ETRS89 / UTM 32N",
        BASEGEOGCRS["ETRS89",
            DATUM["European Terrestrial Reference System 1989",
                ELLIPSOID["GRS 1980",6378137,298.257222101,
                    LENGTHUNIT["metre",1]]],
            PRIMEM["Greenwich",0,
                ANGLEUNIT["degree",0.0174532925199433]],
            ID["EPSG",4258]],
        CONVERSION["UTM zone 32N",
            METHOD["Transverse Mercator",
                ID["EPSG",9807]],
            PARAMETER["Latitude of natural origin",0,
                ANGLEUNIT["degree",0.0174532925199433],
                ID["EPSG",8801]],
            PARAMETER["Longitude of natural origin",9,
                ANGLEUNIT["degree",0.0174532925199433],
                ID["EPSG",8802]],
            PARAMETER["Scale factor at natural origin",0.9996,
                SCALEUNIT["unity",1],
                ID["EPSG",8805]],
            PARAMETER["False easting",500000,
                LENGTHUNIT["metre",1],
                ID["EPSG",8806]],
            PARAMETER["False northing",0,
                LENGTHUNIT["metre",1],
                ID["EPSG",8807]]],
        CS[Cartesian,2],
            AXIS["easting",east,
                ORDER[1],
                LENGTHUNIT["metre",1]],
            AXIS["northing",north,
                ORDER[2],
                LENGTHUNIT["metre",1]],
        ID["EPSG",25832]],
    VERTCRS["DHHN2016",
        VDATUM["Deutsches Haupthoehennetz 2016"],
        CS[vertical,1],
            AXIS["gravity-related height",up,
                LENGTHUNIT["metre",1]],
        ID["EPSG",7837]]]

Theoretically I think they should be in the same CRS. So I try to assign the CRS of LAS2 to LAS1 which gives a warning that the header was not updated but otherwise seems to work at the first glance.

> st_crs(las1) <- st_crs(las2)
Warning message:
EPSG code not found: header not updated. Try to provide the ESPG code instead of a crs 

> st_crs(las1)
Coordinate Reference System:
  User input: COMPD_CS["ETRS89 / ETRS89 / UTM 32N + DHHN2016",PROJCS["ETRS89 / ETRS89 / UTM 32N",GEOGCS["ETRS89",DATUM["European_Terrestrial_Reference_System_1989",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","6258"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4258"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",9],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","25832"]],VERT_CS["DHHN2016",VERT_DATUM["Deutsches Haupthoehennetz 2016",2005,AUTHORITY["EPSG","7837"]],UNIT["metre",1.0,AUTHORITY["EPSG","9001"]],AXIS["Gravity-related height",UP],AUTHORITY["EPSG","7837"]]] 
  wkt:
COMPOUNDCRS["ETRS89 / ETRS89 / UTM 32N + DHHN2016",
    PROJCRS["ETRS89 / ETRS89 / UTM 32N",
        BASEGEOGCRS["ETRS89",
            DATUM["European Terrestrial Reference System 1989",
                ELLIPSOID["GRS 1980",6378137,298.257222101,
                    LENGTHUNIT["metre",1]]],
            PRIMEM["Greenwich",0,
                ANGLEUNIT["degree",0.0174532925199433]],
            ID["EPSG",4258]],
        CONVERSION["UTM zone 32N",
            METHOD["Transverse Mercator",
                ID["EPSG",9807]],
            PARAMETER["Latitude of natural origin",0,
                ANGLEUNIT["degree",0.0174532925199433],
                ID["EPSG",8801]],
            PARAMETER["Longitude of natural origin",9,
                ANGLEUNIT["degree",0.0174532925199433],
                ID["EPSG",8802]],
            PARAMETER["Scale factor at natural origin",0.9996,
                SCALEUNIT["unity",1],
                ID["EPSG",8805]],
            PARAMETER["False easting",500000,
                LENGTHUNIT["metre",1],
                ID["EPSG",8806]],
            PARAMETER["False northing",0,
                LENGTHUNIT["metre",1],
                ID["EPSG",8807]]],
        CS[Cartesian,2],
            AXIS["easting",east,
                ORDER[1],
                LENGTHUNIT["metre",1]],
            AXIS["northing",north,
                ORDER[2],
                LENGTHUNIT["metre",1]],
        ID["EPSG",25832]],
    VERTCRS["DHHN2016",
        VDATUM["Deutsches Haupthoehennetz 2016"],
        CS[vertical,1],
            AXIS["gravity-related height",up,
                LENGTHUNIT["metre",1]],
        ID["EPSG",7837]]]
> st_crs(las1) == st_crs(las2)
[1] TRUE

When I try to merge the pointclouds now I receive an error that the CRS are different

> las3 = rbind(las1, las2, fill=T)
Error: Different CRS.

Another option I tried is to set the crs manually for both files which has the same outcome

> st_crs(las1) <- 25832
> st_crs(las2) <- 25832
> st_crs(las1) == st_crs(las2)
[1] TRUE
> las3 = rbind(las1, las2, fill=T)
Error: Different CRS.

Yet another option I tried is to set the CRS of LAS1 manually and then transform to the CRS of LAS2 with the same outcome

> st_crs(las1) <- 25832
> las1 <- sf::st_transform(las1, crs=st_crs(las2))
Warning message:
EPSG code not found: header not updated. Try to provide the ESPG code instead of a crs
> st_crs(las1) == st_crs(las2)
[1] TRUE
> las3 = rbind(las1, las2, fill=T)
Error: Different CRS.

And vice versa which gives no warning but otherwise the same outcome

> st_crs(las1) <- 25832
> las2 <- sf::st_transform(las2, crs=st_crs(las1))
> st_crs(las1) == st_crs(las2)
[1] TRUE
> las3 = rbind(las1, las2, fill=T)
Error: Different CRS.

Do you have a clue what is the problem? I guess this might be related to the fact that the files are in different LAS standards, since the documentation states that

There are two ways to store the CRS of a point cloud in a LAS file:

  • Store an EPSG code (for LAS 1.0 to 1.3)
  • Store a WTK string (for LAS 1.4)
Jean-Romain commented 10 months ago

Can you produce a minimal reproducible example with data. The following MRE works for me.

LASfile <- system.file("extdata", "Megaplot.laz", package="lidR")
las1 = readLAS(LASfile)
las2 = readLAS(LASfile) 
st_crs(las1) <- 25832
st_crs(las2) <- 25832
st_crs(las1) == st_crs(las2)
rbind(las1,las2)
wiesehahn commented 10 months ago

Yes sure, I created leightweight thinned versions of the files. 

las1 <- readLAS(here::here("las1.laz"))
las2 <- readLAS(here::here("las2.laz"))
st_crs(las1) <- 25832
st_crs(las2) <- 25832
st_crs(las1) == st_crs(las2)
las3 = rbind(las1, las2, fill=T)

lasfiles.zip

Jean-Romain commented 10 months ago

Ho, I got it! fill is not a LAS object. You are trying to use the data.table version of rbind but you are trying to bind LAS objects, not data.table. By the way, your point clouds are not bindable, They have a different number of attributes. The first one is format 1, the second one is format 8 (with Overlap flag, RGB, NIR and 2 extrabytes attributes)

wiesehahn commented 10 months ago

ah Ok thanks, I tried it with the same attributes before but it didnt work for some reason.

data.table's warning confused me as it mentioned to use fill=TRUE and after implementing that the warning was gone and the above mentioned error occured, why I thought it had worked on LAS files to fill the missing columns. 

> las3 = rbind(las1, las2)
Error in data.table::rbindlist(lapply(dots, function(x) { : 
  Item 2 has 25 columns, inconsistent with item 1 which has 16 columns. To fill missing columns use fill=TRUE.
In addition: Warning message:
Different LAS objects have different scales and/or offsets. The first object was used as reference to quantize the others. 
> las3 = rbind(las1, las2, fill=T)
Error: Different CRS.

So filtering the data in advance does the trick, e.g.:

las1 <- readLAS(here::here("las1.laz"), select = "xyzi")
las2 <- readLAS(here::here("las2.laz"), select = "xyzi")
st_crs(las1) <- 25832
st_crs(las2) <- 25832
las3 = rbind(las1, las2)
Jean-Romain commented 10 months ago

I see. I could add an error upstream to prevent data.table to throw a confusing message.