1

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]]]
2
  • 1
    You are showing us the CRS as interpreted by arcgis and they indeed look similar. But the problem is from R. Please show us the CRS as interpreted in R. For the las file you can show the output of readLASheader(). For both you can print the output of st_crs()
    – JRR
    Sep 1 at 22:44
  • @JRR, done, please check Sep 2 at 13:23

1 Answer 1

1

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)
0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.