R lidR clip_roi Error – How to Handle clip_roi Error in R lidR

lidrrsf

I am trying to clip a LAS object based on an sf object. Now both the objects have the same CRS according to ArcGIS Pro 3, but clip_roi still returns an error that the CRS don't match.

Error:

Error in st_geos_binop("intersects", x, y, sparse = sparse, prepared = prepared,  : 
  st_crs(x) == st_crs(y) is not TRUE

Code:

e1586n0472_subset = clip_roi(FL_20180602_B1_e1586n0472, clip_sf)

How can I fix this?

Now when I use the Extract LAS tool in ArcGIS Pro 3, the clip works just fine.

sf object CRS:

Coordinate Reference System:
  User input: NAD83(NSRS2007) / Florida East 
  wkt:
PROJCRS["NAD83(NSRS2007) / Florida East",
    BASEGEOGCRS["NAD83(NSRS2007)",
        DATUM["NAD83 (National Spatial Reference System 2007)",
            ELLIPSOID["GRS 1980",6378137,298.257222101,
                LENGTHUNIT["metre",1]],
            ID["EPSG",6759]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["Degree",0.0174532925199433]]],
    CONVERSION["unnamed",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",24.3333333333333,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",-81,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.999941176470588,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",200000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]]

LAS object CRS:

Coordinate Reference System:
  User input: PROJCS["NAD_1983_NSRS2007_StatePlane_Florida_East_FIPS_0901",GEOGCS["GCS_NAD_1983_NSRS2007",DATUM["D_NAD_1983_NSRS2007",SPHEROID["GRS_1980",6378137.0,298.257222101,AUTHORITY["EPSG",7019]],AUTHORITY["EPSG",6759]],PRIMEM["Greenwich",0.0,AUTHORITY["EPSG",8901]],UNIT["Degree",0.0174532925199433,AUTHORITY["EPSG",9102]],AUTHORITY["EPSG",4759]],PROJECTION["Transverse_Mercator",AUTHORITY["Esri",43006]],PARAMETER["False_Easting",200000.0,AUTHORITY["Esri",100001]],PARAMETER["False_Northing",0.0,AUTHORITY["Esri",100002]],PARAMETER["Central_Meridian",-81.0,AUTHORITY["Esri",100010]],PARAMETER["Scale_Factor",0.9999411764705882,AUTHORITY["Esri",100003]],PARAMETER["Latitude_Of_Origin",24.33333333333333,AUTHORITY["Esri",100021]],UNIT["Meter",1.0,AUTHORITY["EPSG",9001]],AUTHORITY["EPSG",3511]],VERTCS["NAVD_1988",VDATUM["North_American_Vertical_Datum_1988",AUTHORITY["EPSG",5103]],PARAMETER["Vertical_Shift",0.0,AUTHORITY["Esri",100006]],PARAMETER["Direction",1.0,AUTHORITY["Esri",100007]],UNIT["Meter",1.0,AUTHORITY["EPSG",9001]],AUTHORITY["EPSG",5703]] 
  wkt:
COMPOUNDCRS["NAD83(NSRS2007) / Florida East + NAVD88 height",
    PROJCRS["NAD83(NSRS2007) / Florida East",
        BASEGEOGCRS["NAD83(NSRS2007)",
            DATUM["NAD83 (National Spatial Reference System 2007)",
                ELLIPSOID["GRS 1980",6378137,298.257222101,
                    LENGTHUNIT["metre",1]]],
            PRIMEM["Greenwich",0,
                ANGLEUNIT["Degree",0.0174532925199433]],
            ID["EPSG",4759]],
        CONVERSION["unnamed",
            METHOD["Transverse Mercator",
                ID["EPSG",9807]],
            PARAMETER["Latitude of natural origin",24.3333333333333,
                ANGLEUNIT["degree",0.0174532925199433],
                ID["EPSG",8801]],
            PARAMETER["Longitude of natural origin",-81,
                ANGLEUNIT["degree",0.0174532925199433],
                ID["EPSG",8802]],
            PARAMETER["Scale factor at natural origin",0.999941176470588,
                SCALEUNIT["unity",1],
                ID["EPSG",8805]],
            PARAMETER["False easting",200000,
                LENGTHUNIT["metre",1],
                ID["EPSG",8806]],
            PARAMETER["False northing",0,
                LENGTHUNIT["metre",1],
                ID["EPSG",8807]]],
        CS[Cartesian,2],
            AXIS["(E)",east,
                ORDER[1],
                LENGTHUNIT["metre",1]],
            AXIS["(N)",north,
                ORDER[2],
                LENGTHUNIT["metre",1]],
        ID["EPSG",3511]],
    VERTCRS["NAVD88 height",
        VDATUM["North American Vertical Datum 1988"],
        CS[vertical,1],
            AXIS["gravity-related height (H)",up,
                LENGTHUNIT["metre",1]],
        ID["EPSG",5703]]]

Best Answer

According to what I see I can guess that your LAS file is a LAS 1.4 file with a CRS stored as WKT string instead of an EPSG code. The string stored in the LAS file is:

PROJCS["NAD_1983_NSRS2007_StatePlane_Florida_East_FIPS_0901",GEOGCS["GCS_NAD_1983_NSRS2007",DATUM["D_NAD_1983_NSRS2007",SPHEROID["GRS_1980",6378137.0,298.257222101,AUTHORITY["EPSG",7019]],AUTHORITY["EPSG",6759]],PRIMEM["Greenwich",0.0,AUTHORITY["EPSG",8901]],UNIT["Degree",0.0174532925199433,AUTHORITY["EPSG",9102]],AUTHORITY["EPSG",4759]],PROJECTION["Transverse_Mercator",AUTHORITY["Esri",43006]],PARAMETER["False_Easting",200000.0,AUTHORITY["Esri",100001]],PARAMETER["False_Northing",0.0,AUTHORITY["Esri",100002]],PARAMETER["Central_Meridian",-81.0,AUTHORITY["Esri",100010]],PARAMETER["Scale_Factor",0.9999411764705882,AUTHORITY["Esri",100003]],PARAMETER["Latitude_Of_Origin",24.33333333333333,AUTHORITY["Esri",100021]],UNIT["Meter",1.0,AUTHORITY["EPSG",9001]],AUTHORITY["EPSG",3511]],VERTCS["NAVD_1988",VDATUM["North_American_Vertical_Datum_1988",AUTHORITY["EPSG",5103]],PARAMETER["Vertical_Shift",0.0,AUTHORITY["Esri",100006]],PARAMETER["Direction",1.0,AUTHORITY["Esri",100007]],UNIT["Meter",1.0,AUTHORITY["EPSG",9001]],AUTHORITY["EPSG",5703]]

that is interpreted by sf using st_crs(). This is the content you copy/pasted above. This CRS is a compound CRS with XY being NAD83(NSRS2007) / Florida East and Z being NAVD88 height.

Your vector file on another hand only contains XY data with a simple CRS which is NAD83(NSRS2007) / Florida East.

The two CRS are indeed different. sf performs string comparisons and strings mismatch even if they are equivalent for XY. This is not something I, as the lidR developer, can improve. This is something that must come from sf.

The simplest workaround is to assign the CRS of the vector to the point cloud

st_crs(FL_20180602_B1_e1586n0472) = st_crs(clip_sf)
Related Question