pypsa-eur/scripts/build_sequestration_potentials.py

import pandas as pd
import geopandas as gpd

def area(gdf):
    """Returns area of GeoDataFrame geometries in square kilometers."""
    return gdf.to_crs(epsg=3035).area.div(1e6)


def allocate_sequestration_potential(gdf, regions, attr='conservative estimate Mt', threshold=3):
    gdf = gdf.loc[gdf[attr] > threshold, [attr, "geometry"]]
    gdf["area_sqkm"] = area(gdf)
    overlay = gpd.overlay(regions, gdf, keep_geom_type=True)
    overlay["share"] = area(overlay) / overlay["area_sqkm"]
    adjust_cols = overlay.columns.difference({"name", "area_sqkm", "geometry", "share"})
    overlay[adjust_cols] = overlay[adjust_cols].multiply(overlay["share"], axis=0)
    gdf_regions = overlay.groupby("name").sum()
    gdf_regions.drop(["area_sqkm", "share"], axis=1, inplace=True)
    return gdf_regions.squeeze()


if __name__ == "__main__":
    if 'snakemake' not in globals():
        from helper import mock_snakemake
        snakemake = mock_snakemake(
            'build_sequestration_potentials',
            simpl='',
            clusters="181"
        )

    cf = snakemake.config["sector"]["regional_co2_sequestration_potential"]

    gdf = gpd.read_file(snakemake.input.sequestration_potential[0])

    regions = gpd.read_file(snakemake.input.regions_offshore)
    if cf["include_onshore"]:
        onregions = gpd.read_file(snakemake.input.regions_onshore)
        regions = pd.concat([regions, onregions]).dissolve(by='name').reset_index()

    s = allocate_sequestration_potential(gdf, regions, attr=cf["attribute"], threshold=cf["min_size"])

    s = s.where(s>cf["min_size"]).dropna()

    s.to_csv(snakemake.output.sequestration_potential)