在 for 循环内对 for 循环中的数据帧值求和

Question

我有一个large polygon文件，small polygon文件，还是points文件。我在这里所做的是循环遍历大多边形以查找哪些小多边形相交。然后计算大多边形内每个小多边形的面积。然后我循环遍历小多边形以查找每个多边形中的点统计数据。

我发现了number_of_somethin每个小多边形的价值。问题是如何对number_of_somethin大多边形内的所有小多边形值求和并将结果存储在原始large_polygon文件中作为新列，比方说large_polygon['smth_sum']？

我df_res_2.loc[idx, 'smth'] = number_of_somethin在大多边形内的每个小多边形中得到了 number_of_somethin 值。现在我需要将它们相加large_polygon['smth_sum']

注： FID为大多边形的id，ID为小多边形的id

import geopandas as gpd

small_polygon = gpd.read_file(r'R:\...\small.shp')
large_polygon = gpd.read_file(r'R:\...\large.shp')
points = gpd.read_file(r'R:\...\points.shp')

SmallJoin =gpd.sjoin(small_polygon, large_polygon)[['FID', 'ID', 'someValue','geometry']]

for i in large_polygon.index:
    df_i = SmallJoin[SmallJoin['FID'] == i]

    # i do something here, f.e. calculate small polgyon area
    df_res = gpd.overlay(large_polygon, df_i, how='intersection')
    df_res['area'] = round((df_res.apply(lambda row: row.geometry.area, axis=1)), 4)

    # now i know area for each small polygon within large polygon
    df_res_2 = df_res[df_res['FID_1'] == i]

    # now point statistics in small polygons
    PointsJoin =gpd.sjoin(points, df_res)[['ID','someAttribute', 'someAttribute2','geometry']]

    for idx, val in df_res_2['ID'].items():
        df_idx = PointsJoin[PointsJoin['ID'] == val]
        number_of_somethin = df_idx ['someAttribute'] + 121 + df_idx['someAttribute2']
        df_res_2.loc[idx, 'smth'] = number_of_somethin

我有一些想法如何做到这一点，但没有一个是行不通的，所以我认为还有其他方法。

large_polygon.loc[i, 'smth_sum'] = df_res_2['smth']
large_polygon.loc[i, 'smth_sum'] = df_res_2['smth'].sum()

large_polygon['smth_sum'] = large_polygon[large_polygon['FID'] == df_res_2['FID_1'].sum()]

Answer 1

• 你描述了三个GeoDataFrame

  1. 大 - 为此使用了国家边界
  2. 小 - 为此使用了 UTM 区域边界
  3. 点 - 使用随机生成的点，大部分重叠 2

• 您定义每个大型几何图形需要两个输出（此处为国家/地区）

  • 面积-每个小几何图形的交叉 面积之和
  • value -空间连接到大几何图形的小几何图形内的点的值之和

• 上述所有内容都可以通过空间连接和pandas 来实现merge()groupby()

• 为了使这一点更清楚 - 还包括一种可视化所有这些的方法

import geopandas as gpd
import shapely.geometry
import requests
import numpy as np
import plotly.express as px

# get some sample data....
# fmt: off
gdf_utm = gpd.GeoDataFrame.from_features(requests.get("https://opendata.arcgis.com/datasets/b294795270aa4fb3bd25286bf09edc51_0.geojson").json()).set_crs("EPSG:4326")
gdf_countries = gpd.read_file(gpd.datasets.get_path("naturalearth_lowres"))

large_polygon = gdf_countries.loc[lambda d: d["iso_a3"].isin(["BEL", "LUX", "NLD", "DEU", "AUT"])]
# large_polygon.boundary.plot()

small_polygon = gpd.sjoin(gdf_utm, large_polygon).loc[:, gdf_utm.columns].groupby(["FID", "ZONE"]).first().reset_index()
# fmt: on

# some points within geometry of small_polygon
b = small_polygon.total_bounds
POINTS = 10
points = gpd.GeoDataFrame(
    geometry=[
        shapely.geometry.Point(x, y)
        for x, y in zip(
            np.random.uniform(*b[[0, 2]], POINTS),
            np.random.uniform(*b[[1, 3]], POINTS),
        )
    ],
    data={"value": np.arange(POINTS)},
    crs="EPSG:4326",
)

# spatial small to large with geometry from large
SmallJoin = gpd.sjoin(small_polygon, large_polygon).merge(
    large_polygon["geometry"],
    left_on="index_right",
    right_index=True,
    suffixes=("", "_large"),
)
SmallJoin["area"] = SmallJoin.intersection(gpd.GeoSeries(SmallJoin["geometry_large"])).area

# get sums of area of overlap and sum of values from points
Final = (
    SmallJoin.rename(columns={"index_right": "index_large"})
    .sjoin(points)
    .groupby("index_large")
    .agg({"area": "sum", "value": "sum", "geometry_large": "first"})
)

输出

索引大	区域	价值
114	24.6382	25
121	90.3565	45
128	0.603031	20
129	7.65999	20
130	10.5284	20

想象它

px.choropleth_mapbox(
    Final,
    geojson=gpd.GeoSeries(Final["geometry_large"]),
    locations=Final.index,
    color="value",
    hover_data=["area"],
).add_traces(
    px.scatter_mapbox(
        points,
        lat=points.geometry.y,
        lon=points.geometry.x,
        color="value",
    )
    .update_traces(marker_coloraxis="coloraxis2", marker_size=10)
    .data
).update_layout(
    mapbox={
        "style": "carto-positron",
        "center": {"lon": sum(b[[0, 2]]) / 2, "lat": sum(b[[1, 3]]) / 2},
        "zoom": 3,
        "layers": [{"source": small_polygon.__geo_interface__, "type": "line"}],
    },
    coloraxis2={
        "colorbar": {"x": -0.1, "title": "scatter"},
        "colorscale": [[0, "blue"], [1, "blue"]],
    },
    coloraxis={"colorscale": [[0, "white"], [1, "green"]]},
    margin={"l": 0, "r": 0, "t": 0, "b": 0},
)