空间连接后基于栅格面积权重计算各国PDSI空间平均值的方法咨询

我正在使用下方代码从坐标数据中提取所属国家信息，相关代码说明可参考公开技术社区的过往教程。我使用的核心变量是NOAA公开数据集的月度平均PDSI数值，代码生成的可视化结果局部中，阴影方块代表PDSI数值对应的空间区域，叠加在世界shapefile图层上。

以比利时为例，与比利时陆域相交的4个栅格方块同时也覆盖其他国家，若直接将这些栅格的PDSI值全部计入比利时的均值计算，会导致结果偏高，尤其是下方两个栅格仅极少面积覆盖比利时，计算均值时这部分数值的权重应大幅降低。请问是否可以实现以栅格落在对应国家境内的面积为权重的加权平均计算？同时我希望该方法可以标准化适配所有国家的计算需求。

原有代码如下：

import geopandas as gpd
import numpy as np
import plotly.express as px
import requests
from pathlib import Path
from zipfile import ZipFile
import urllib
import shapely.geometry
import xarray as xr

# download NetCDF data...
# fmt: off
url = "https://psl.noaa.gov/repository/entry/get/pdsi.mon.mean.selfcalibrated.nc?entryid=synth%3Ae570c8f9-ec09-4e89-93b4-babd5651e7a9%3AL2RhaV9wZHNpL3Bkc2kubW9uLm1lYW4uc2VsZmNhbGlicmF0ZWQubmM%3D"
f = Path.cwd().joinpath(Path(urllib.parse.urlparse(url).path).name)
# fmt: on

if not f.exists():
    r = requests.get(url, stream=True, headers={"User-Agent": "XY"})
    with open(f, "wb") as fd:
        for chunk in r.iter_content(chunk_size=128):
            fd.write(chunk)
ds = xr.open_dataset(f)
pdsi = ds.to_dataframe()
pdsi = pdsi.reset_index().dropna()  # don't care about places in oceans...

# use subset for testing... last 5 times...
pdsim = pdsi.loc[pdsi["time"].isin(pdsi.groupby("time").size().index[-5:])]

# create geopandas dataframe
gdf = gpd.GeoDataFrame(
    pdsim, geometry=pdsim.loc[:, ["lon", "lat"]].apply(shapely.geometry.Point, axis=1)
)

# make sure that data supports using a buffer...
assert (
    gdf["lat"].diff().loc[lambda s: s.ne(0)].mode()
    == gdf["lon"].diff().loc[lambda s: s.ne(0)].mode()
).all()
# how big should the square buffer be around the point??
buffer = gdf["lat"].diff().loc[lambda s: s.ne(0)].mode().values[0] / 2
gdf["geometry"] = gdf["geometry"].buffer(buffer, cap_style=3)

# Import shapefile from geopandas
path_to_data = gpd.datasets.get_path("naturalearth_lowres")
world_shp = gpd.read_file(path_to_data)

# the solution... spatial join buffered polygons to countries
# comma separate associated countries
gdf = gdf.join(
    world_shp.sjoin(gdf.set_crs("EPSG:4326"))
    .groupby("index_right")["name"]
    .agg(",".join)
)
gdf["time_a"] = gdf["time"].dt.strftime("%Y-%b-%d")

# simplest way to test is visualise...
px.choropleth_mapbox(
    gdf,
    geojson=gdf.geometry,
    locations=gdf.index,
    color="pdsi",
    hover_data=["name"],
    animation_frame="time_a",
    opacity=.3
).update_layout(
    mapbox={"style": "carto-positron", "zoom": 1},
    margin={"l": 0, "r": 0, "t": 0, "b": 0},
)

解决方案 #

该需求可以直接实现，且天然支持所有国家的标准化计算，你只需要将原有代码中空间连接后聚合国家名的逻辑，替换为以下面积加权计算逻辑即可：

# 为栅格数据设置坐标系
gdf = gdf.set_crs("EPSG:4326")
# 空间连接保留所有相交的栅格-国家配对
sjoined = gpd.sjoin(gdf, world_shp, how='inner', predicate='intersects')

# 转换为等面积投影计算面积，避免经纬度投影的高纬度面积失真
equal_area_crs = "+proj=cea"
# 计算每个栅格与对应国家的交集面积
sjoined["intersection_area"] = sjoined.apply(
    lambda row: row["geometry"].intersection(world_shp.loc[row["index_right"], "geometry"]).to_crs(equal_area_crs).area,
    axis=1
)
# 计算单个栅格的总面积
sjoined["cell_total_area"] = sjoined["geometry"].to_crs(equal_area_crs).area
# 计算权重：交集面积 / 栅格总面积
sjoined["weight"] = sjoined["intersection_area"] / sjoined["cell_total_area"]

# 按国家、时间分组计算加权平均PDSI
country_monthly_pdsi = sjoined.groupby(["name", "time"]).apply(
    lambda group: np.average(group["pdsi"], weights=group["weight"])
).reset_index(name="weighted_avg_pdsi")

逻辑说明 #

• 所有相交的栅格和国家配对都会被保留，不会遗漏任何覆盖关系
• 采用等面积投影计算面积，避免了WGS84经纬度坐标系下不同纬度面积换算的误差
• 权重完全由栅格落在一国境内的面积占比自动决定，覆盖面积越小的栅格对该国均值的影响越低，完全解决你提到的小面积覆盖栅格干扰计算结果的问题
• 无需针对单个国家做任何定制化调整，所有国家的计算逻辑完全统一

内容的提问来源于stack exchange，提问作者ggmaster