计算空间连接后每个国家/地区的空间平均值

Posted 2023-03-23

【中文标题】计算空间连接后每个国家/地区的空间平均值

【英文标题】：Calculating spatial averages for each country after spatial join

【发布时间】：2021-12-31 05:24:57

【问题描述】：

您好，我在底部使用以下代码从坐标中提取国家/地区。请参阅以下网址，该网址提供了更详细的代码说明：Extracting countries from NetCDF data using geopandas。

我的主要变量/值是来自https://psl.noaa.gov/data/gridded/data.pdsi.html 的月平均 pdsi 值。下图表示由下面的代码创建的可视化的一部分。阴影方块表示 pdsi 值的空间区域，与世界的 shapefile 重叠。

从比利时的图片可以看出，与比利时国土面积相接的4个方格也与其他国家相接。如果我将基值归因于比利时，我相信这会高估平均 pdsi 值。尤其是考虑到底部的两个方块几乎不接触比利时时，这些值在计算平均值时的权重应该显着降低。因此，有没有一种方法可以结合某种加权平均值，其中一个国家内每个正方形的面积可以用作调整每个 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,
)

【问题讨论】：

这听起来像是一个科学问题，而不是一个编码问题。 SO 适用于您陈述您的方法并且您需要帮助编码该方法的问题

【参考方案1】：

使用https://geopandas.org/en/stable/docs/reference/api/geopandas.GeoSeries.intersection.html，您可以获得与国家多边形相交的部分网格 使用面积，可以计算重叠比例 由此我生成了两个可视化

显示网格重叠的国家/地区以及重叠程度 使用加权平均值汇总到国家/地区并计算可用于提高透明度的其他指标

我不知道以这种方式聚合 PDSI（平均值或加权平均值）在数学/科学上是否合理。这确实演示了如何获得您的问题请求的结果。

# the solution... spatial join buffered polygons to countries
# plus work out overlap between PDSI grid and country.  Area of each grid is constant...
gdf_c = (
    world_shp.sjoin(gdf.set_crs("EPSG:4326"))
    .merge(
        gdf.loc[:, "geometry"],
        left_on="index_right",
        right_index=True,
        suffixes=("", "_pdsi"),
    )
    .assign(
        overlap=lambda d: (
            d["geometry"]
            .intersection(gpd.GeoSeries(d["geometry_pdsi"], crs="EPSG:4326"))
            .area
            / (buffer * 2) ** 2
        ).round(3)
    )
)

# comma separate associated countries and a list of overlaps
gdf_pdsi = gdf.loc[:, ["geometry", "time", "pdsi"]].join(
    gdf_c.groupby("index_right").agg("name": ",".join, "overlap": list)
)
gdf_pdsi["time_a"] = gdf_pdsi["time"].dt.strftime("%Y-%b-%d")

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

fig_buf

import pandas as pd

# prepare data to plot by country
df_pdsi = (
    gdf_c.groupby(["name", "time"])
    .apply(
        lambda d: pd.Series(
            
                "weighted_pdsi": (d["pdsi"] * d["overlap"]).sum() / d["overlap"].sum(),
                "unweighted_pdsi": d["pdsi"].mean(),
                "min_pdsi": d["pdsi"].min(),
                "max_pdsi": d["pdsi"].max(),
                "min_overlap": d["overlap"].min(),
                "max_overlap": d["overlap"].max(),
                "size_pdsi": len(d["pdsi"]),
                # "pdsi_list":[round(v,2) for v in d["pdsi"]]
            
        )
    )
    .reset_index()
)
df_pdsi["time_a"] = df_pdsi["time"].dt.strftime("%Y-%b-%d")
fig = px.choropleth_mapbox(
    df_pdsi,
    geojson=world_shp.set_index("name").loc[:, "geometry"],
    locations="name",
    color="weighted_pdsi",
    hover_data=df_pdsi.columns,
    animation_frame="time_a",
    opacity=0.3,
).update_layout(
    mapbox="style": "carto-positron", "zoom": 1,
    margin="l": 0, "r": 0, "t": 0, "b": 0,
)
fig

【讨论】：

嗨@RobRaymond，非常感谢您再次提供帮助。您的代码成功运行，但是当我尝试将 df_pdsi 数据框导出到 csv 文件时，仅列出了 2014 年。我需要完整的时间序列数据，但我无法弄清楚如何添加完整的年份。你能帮我解决这个问题吗？

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

抱歉，我不确定在哪里包含此代码？我确实在开头包含了这段代码并运行了你发布的新代码。在运行 df_pdsi 代码之前，我应该第二次使用相同的代码吗？

看看你在问题中发布的代码......你会看到早期有一行明确评论使用子集进行测试......最后5次。 .. 删除/评论这一行，它将适用于所有数据。这是一个熊猫过滤器pandas.pydata.org/docs/getting_started/intro_tutorials/…。我认为这些图表不适用于所有日期，情节数据过多/动画中的日期过多