Am kind of new to the whole Geo Info systems, I am trying to clip USGS lidar data which is (.laz) format using polygons in a (.shp) file. I am using Python 3.10: with packages laspy, geopandas. I use laspy for reading the pointcloud, then geopandas for clipping and other processing. But I am unable to clip them for some reason. I suspect it is because of the CRS values, which for some reason I am not able to find in the USGS points. But the polygons have a CRS. This is the code I am using for this.

with laspy.open('USGS_LPC_MA_CentralEastern_2021_B21_19TCH324723.laz') as fh:
    print('Points from Header:', fh.header.point_count)
    las_file = fh.read()
polygons = gpd.read_file("/home/tarun/Downloads/structures_poly/structures_poly_9.shp")
# print(polygons.crs)
points = gpd.GeoDataFrame(geometry=gpd.points_from_xy(las_file.x, las_file.y))
clipped_points = points.clip(polygons)


This code shows an error :

UserWarning: CRS mismatch between the CRS of left geometries and the CRS of right geometries. Use to_crs() to reproject one of the input geometries to match the CRS of the other.

Left CRS: None
Right CRS: EPSG:26986

How can I get the CRS from the LAZ file ?

Update from a previous stack: I tried using the header header.vlrs string, but that had multiple entries, was not sure which one to use as the CRS. Attaching the header string below header.vlrs

  • Both datasets need to be in (or set) to the same projection 'NAD83 / Massachusetts Mainland' epsg.io/26986
    – Mapperz
    Commented Jun 29, 2023 at 20:58
  • try gis.stackexchange.com/a/434219/276
    – Mapperz
    Commented Jun 29, 2023 at 21:01
  • Yeah I tried to follow that stack as it is, but the header file has multiple EPSG values. I was confused on which one to use. I have updated the question with the string output. Thanks for the help !! Commented Jun 30, 2023 at 16:44


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.