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() print(las_file.header.vlrs[-1].string) 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) print(clipped_points.shape)
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 ?