Python：如何将 geotiff 转换为 geopandas？答案

【问题标题】：Python: how to convert geotiff to geopandas?Python：如何将 geotiff 转换为 geopandas？
【发布时间】：2021-08-01 21:02:24
【问题描述】：

我有一个 geotiff 文件。

import xarray as xr
urbanData = xr.open_rasterio('myGeotiff.tif')
plt.imshow(urbanData)

这里是link 文件。

我可以将文件转换为以坐标为点的数据框

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

但是，我想获得一个数据框，其中包含像素的几何形状作为多边形而不是点。有可能吗？

【问题讨论】：

“多边形的像素几何”是什么意思？ GeoTIFF 是光栅文件，而不是矢量格式，所以我怀疑它会存储多边形数据。
@jjramsey 这就是重点。是否可以将光栅文件转换为矢量文件？
这很冒险。对于art，有一些工具可以通过执行边界跟踪并将跟踪的路径存储为矢量来从光栅文件制作矢量文件。但是，该场景中的光栅文件相对简单，通常包含曲线或纯色区域的描述。这种矢量追踪的合适候选者是扫描黑白卡通画，而不是通常在 GeoTIFF 中的那种内容。
听起来你想多边形化你的栅格。如果是这样，this 可能是一个好的开始

标签： python geopandas python-xarray geotiff rasterio

【解决方案1】：

令我惊讶的是，我还没有真正找到一个包装 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)

【讨论】：

此方法不会创建带孔的多边形。如果您想这样做，请使用geometries.append(shape(geom)) 而不是geometries.append(sg.Polygon(geom["coordinates"][0]))|