在Python中通过光栅范围剪裁地理数据帧点

Question

背景

感谢善良的灵魂在这里回答了我之前的问题：从Dask或Pandas数据帧创建地理数据帧的更快方法，我成功地将我的大型块模型导入为地理数据帧。这是用以下代码完成的：

import dask_geopandas 
import dask
from dask import dataframe as dd
import geopandas as gpd  
BM = dd.read_csv(BM_path, skiprows=2, names=["X","Y","Z","Lith"])
BM["geometry"] = dask_geopandas.points_from_xy(BM,"X","Y")
BM = dask_geopandas.from_dask_dataframe(BM,geometry="geometry")
BM = BM.compute()

我甚至能够在这里找到代码，将我的光栅转换为地理数据框架！：

import rasterio as rio
with rio.Env():
with rio.open(RAS_path) as src:
crs = src.crs
xmin, ymax = np.around(src.xy(0.00,0.00),9)
xmax, ymin = np.around(src.xy(src.height-1, src.width-1),9)
x = np.linspace(xmin, xmax, src.width)
y = np.linspace(ymax, ymin, src.height)

xs, ys = np.meshgrid(x,y)
zs = src.read(1)

mask = src.read_masks(1) > 0
xs, ys, zs = xs[mask], ys[mask], zs[mask]

data = {'X': pd.Series(xs.ravel()),
'Y': pd.Series(ys.ravel()),
'Z': pd.Series(zs.ravel())}

df = pd.DataFrame(data=data)
geometry = gpd.points_from_xy(df.X,df.Y)
gdf = gpd.GeoDataFrame(df, crs=crs, geometry=geometry)

虽然我不知道如何使用这个地理数据框架来剪裁我的块模型，所以它可能没有用。。。

我知道如何在ArcGIS中做到这一点，但它又慢又笨重。我想在python中简化整个过程，甚至不打开Esri产品。

问题

我一直在努力寻找一种可靠的方法，根据我从ArcGIS导出的光栅的X-Y范围来剪裁我的块模型(BM((该光栅具有浮动类型值，不确定是否相关(。

使用clip()我运气不好，因为我认为光栅不能用于在此工具中剪裁点(必须转换为多边形？(。我想知道是否有人知道将我的地理数据帧点(作为常规数据帧或地理数据帧(剪裁到光栅范围的可靠方法？

起初，我认为我可以从光栅范围创建一个平坦的单个多边形，并将其与clip()工具一起使用，但遗憾的是，我也不知道如何打开光栅并更改值，而不将其作为列表读取，并使用rasterio.open()和rasterio.read()

Answer 1

您希望如何使用网格？是否要使用网格的整个延伸？还是网格内部有遮罩？

如果你想进行空间剪辑，你可能需要创建一个表示网格扩展边界的向量(其他方法需要同时使用两个数据集的X和Y(。如果要使用整个光栅的延伸，请使用rasterio(例如使用src.transform(获取光栅的角点。使用这些值创建一个形状优美的多边形。https://shapely.readthedocs.io/en/stable/manual.html#Polygon

则BM.clip([polygon])

如果要在地理标志内使用不规则遮罩进行剪裁，请使用rasterio.features.shapes从光栅创建多边形。https://rasterio.readthedocs.io/en/latest/api/rasterio.features.html#rasterio.features.shapes你可以将这些转换为地质标准。易于使用的GeoDataFrame。geopandas.GeoDataFrame.from_features(geoms)

您还需要决定如何处理网格内任何可能的孔。如果你只想要外壳，你可能可以尝试只得到最大的多边形。