问题描述
import xarray as xr
urbanData = xr.open_Rasterio('myGeotiff.tif')
plt.imshow(urbanData)
我可以将文件转换为数据框,坐标为点
ur = xr.DataArray(urbanData,name='myData')
ur = ur.to_dataframe().reset_index()
gdfur = gpd.GeoDataFrame(ur,geometry=gpd.points_from_xy(ur.x,ur.y))
但是我想获得一个数据框,其中包含像素的几何形状作为多边形而不是点。可能吗?
解决方法
有点出乎我的意料,我还没有真正找到包装 rasterio.features
以获取 DataArrays 并生成 GeoDataFrames 的包。
这些可能非常有用:
https://corteva.github.io/geocube/stable/
https://corteva.github.io/rioxarray/stable/
我通常使用这样的东西:
import affine
import geopandas as gpd
import rasterio.features
import xarray as xr
import shapely.geometry as sg
def polygonize(da: xr.DataArray) -> gpd.GeoDataFrame:
"""
Polygonize a 2D-DataArray into a GeoDataFrame of polygons.
Parameters
----------
da : xr.DataArray
Returns
-------
polygonized : geopandas.GeoDataFrame
"""
if da.dims != ("y","x"):
raise ValueError('Dimensions must be ("y","x")')
values = da.values
transform = da.attrs.get("transform",None)
if transform is None:
raise ValueError("transform is required in da.attrs")
transform = affine.Affine(*transform)
shapes = rasterio.features.shapes(values,transform=transform)
geometries = []
colvalues = []
for (geom,colval) in shapes:
geometries.append(sg.Polygon(geom["coordinates"][0]))
colvalues.append(colval)
gdf = gpd.GeoDataFrame({"value": colvalues,"geometry": geometries})
gdf.crs = da.attrs.get("crs")
return gdf
请注意,在使用 xr.open_rasterio
读取后,您应该先从 xarray 中挤出波段维度以使其成为 2D:
urbanData = xr.open_rasterio('myGeotiff.tif').squeeze('band',drop=True)