使用并行 dask GeoPandas 检查两个数据帧之间的重叠区域

Question

我有两个不同的 GeoDataFrames：其中一个包含一个大网格中的多边形正方形。另一个包含更大和更少的多边形。 我希望计算每个网格正方形与其他较大正方形的重叠面积。

为此，我做了一个简单的循环方法

for _, patch in tqdm(layer.iterrows(), total=layer.shape[0], desc=name):
    # Index of intersecting squares
    idx = joined.intersects(patch.geometry)
    intersection_polygon = joined[idx].intersection(patch.geometry)
    area_of_intersection = intersection_polygon.area
    joined.loc[idx, "value"] += area_of_intersection

为了加快这种方法的速度，我将包含较大补丁的layer DataFrame 转换为 Dask-DataFrame。

我通过以下方式实现它：

def multi_area(patch, joined=None):
    # Index of intersecting squares
    idx = joined.intersects(patch.geometry)
    intersection_polygon = joined[idx].intersection(patch.geometry)
    area_of_intersection = intersection_polygon.area
    joined.loc[idx, "value"] += area_of_intersection
    return joined["value"]

layer_dask = dask_geopandas.from_geopandas(layer, npartitions=8)

with ProgressBar():
    joined["value"] = layer_dask.apply(multi_area, meta=joined, joined=joined, axis=1).compute(scheduler='multiprocessing')

但是，这会返回错误 AttributeError: 'GeoDataFrame' object has no attribute 'name'，此时我不确定这是否是最佳方法，以及我做错了什么。

我将要做的工作将有 4 亿个方格，所以我计划稍后在较小的区域上批量计算，因为我想不出更聪明的方法来做这件事......

Answer 1

正如迈克尔在 cmets 中所建议的那样，我设法使用空间连接和覆盖来加速这个过程。 此外，我实现了 Dask Dataframes，因此最终代码变为：

import dask_geopandas as dg
import geopandas as gpd

def dissolve_shuffle(ddf, by=None, **kwargs):
    """Shuffle and map partition"""
    meta = ddf._meta.dissolve(by=by, as_index=False, **kwargs)

    shuffled = ddf.shuffle(
        by, npartitions=ddf.npartitions, shuffle="tasks", ignore_index=True
    )

    return shuffled.map_partitions(
        gpd.GeoDataFrame.dissolve, by=by, as_index=False, meta=meta, **kwargs
    )


def calculate_area_overlap_dask(
    df_grid,
    layer,
    nthreads=8,
) -> gpd.GeoDataFrame:
    """This function calculates the area of overlap in each grid cell for a given map-layer
    """

    layer = layer[["geometry"]]
    df_grid = df_grid[["geometry"]]

    # Split up the layer using the grid
    _overlay = gpd.overlay(layer, df_grid, how="intersection")
    
    # Convert the overlay to a dask geopandas dataframe and calculate the area of each new polygon
    _overlay = dg.from_geopandas(_overlay, npartitions=nthreads)
    _overlay["area"] = _overlay.area
    _overlay = _overlay.compute()
    
    # Convert the grid to a dask geopandas dataframe and spatial join all split layer polygons to corresponding grid cells
    df_grid = dg.from_geopandas(df_grid, npartitions=nthreads)
    joined = dg.sjoin(df_grid, _overlay, how="inner").reset_index()

    # Faster dissolve of area within each grid cell
    scored_grid = dissolve_shuffle(
        joined,
        "index",
    )
    scored_grid = scored_grid.compute()
    return scored_grid

def polygon_to_grid(name: str, gdf) -> gpd.GeoDataFrame:
    """This function converts a geodataframe to a grid of polygons
    """
  
    gdf["value"] = range(len(gdf.index))

    # Rasteriser polygonet
    out_grid: xr.Dataset = make_geocube(
        vector_data=gdf,
        measurements=["value"],
        resolution=(-100, 100),
        fill=np.nan,
    )

    vals: xr.DataArray = out_grid.value.values
    vals[~np.isnan(vals)] = np.arange(len(vals[~np.isnan(vals)]), dtype=np.int32)
    vals[np.isnan(vals)] = -9999
    out_grid.value.values = vals
    out_grid.rio.to_raster( f"{name}_raster.tif")
   
    # Read saved raster
    src: xr.Dataset = rasterio.open(f"{name}_raster.tif")
    r = src.read(1).astype(np.int32)

    # Convert polygons
    shapes = features.shapes(r, mask=r != -9999, transform= src.transform)
    polygons: list[Polygon] = list(shapes)
    geom: list[Polygon] = [shapely.geometry.shape(i[0]) for i in polygons]

    # Convert to geodataframe
    grid = gpd.GeoDataFrame(
        geometry=gpd.GeoSeries(
            geom,
        ),
    )
    return grid

if __name__=="__main__":
    area = gpd.read_file("some_area.shp")
    layer = gpd.read_file("some_map_layer.shp")
    area_grid = polygon_to_grid("area", area)
    grid_evaluated = calculate_area_overlap_dask(area_grid, layer)

这种混乱最终奏效了，但它很容易出现大型数据集的内存问题。所以我选择了一个不太精确但速度更快的解决方案。