Python：如何将GeoTIFF转换为GeoPandas？

Question

我有一个Geotiff文件。

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

不过我希望获得一个包含像多边形一样的像素几何体而非点的数据框，这可行吗？

- emax

---

“像素几何形状”是什么意思？GeoTIFF是一种栅格文件，而不是矢量格式，因此我怀疑它是否会存储多边形数据。 - jjramsey

@jjramsey 这就是关键。能否将光栅文件转换为矢量文件？ - emax

那有点棘手。对于艺术来说，可以通过边界跟踪并将路径存储为矢量来从光栅图像生成矢量文件的工具。然而，在这种情况下的光栅文件相对简单，通常包含曲线或实色区域的描绘。适合进行此类矢量跟踪的良好候选对象是黑白卡通图画的扫描版，而不是通常出现在GeoTIFF中的内容。 - jjramsey

听起来你想将栅格转换为多边形。如果是这样，这个链接可能是一个不错的开始。 - Val

Answer 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 读取数据后，您应该先将带宽尺寸挤压成2D。

urbanData = xr.open_rasterio('myGeotiff.tif').squeeze('band', drop=True)

- Huite Bootsma

---

该方法不会创建带洞的多边形。如果您想这样做，请改用geometries.append(shape(geom))而不是geometries.append(sg.Polygon(geom["coordinates"][0]))。 - Hugo Roussaffa

Answer 2

使用 Geocube 0.4+ 版本：

import rioxarray
from geocube.vector import vectorize

data = rioxarray.open_rasterio("myGeotiff.tif").squeeze()
data.name = "myData"
gdf = vectorize(data)

- snowman2