使用 geopandas (Python) 从空间数据框中进行空间分箱

Question

我想做一个空间分箱（使用中值作为聚合函数） 从包含在经纬度位置测量的污染物值的 CSV 文件开始。
生成的地图应该是这样的：

但对于应用于城市范围的数据。 在这方面，我发现这个tutorial 接近我想要做的，但我无法得到想要的结果。 我认为我缺少有关如何正确使用dissolve 并绘制结果数据的内容（最好使用Folium） 有什么有用的示例代码吗？

Answer 1

• 您没有提供样本数据。所以我使用全球地震作为加利福尼亚的点和几何范围/范围
• 使用shapely.geometry.box() 创建网格很简单
• 我展示了使用 median 和另一个 aggfunc 来演示可以计算多个指标
• 已使用 folium 进行绘图。此功能是 geopandas 0.10.0 https://geopandas.org/en/stable/docs/user_guide/interactive_mapping.html 中的新功能

import geopandas as gpd
import shapely.geometry
import numpy as np

# equivalent of CSV, all earthquake points globally
gdf_e = gpd.read_file(
    "https://earthquake.usgs.gov/earthquakes/feed/v1.0/summary/all_month.geojson"
)

# get geometry of bounding area.  Have selected a state rather than a city
gdf_CA = gpd.read_file(
    "https://raw.githubusercontent.com/glynnbird/usstatesgeojson/master/california.geojson"
).loc[:, ["geometry"]]

BOXES = 50
a, b, c, d = gdf_CA.total_bounds

# create a grid for Califormia, could be a city
gdf_grid = gpd.GeoDataFrame(
    geometry=[
        shapely.geometry.box(minx, miny, maxx, maxy)
        for minx, maxx in zip(np.linspace(a, c, BOXES), np.linspace(a, c, BOXES)[1:])
        for miny, maxy in zip(np.linspace(b, d, BOXES), np.linspace(b, d, BOXES)[1:])
    ],
    crs="epsg:4326",
)

# remove grid boxes created outside actual geometry
gdf_grid = gdf_grid.sjoin(gdf_CA).drop(columns="index_right")

# get earthquakes that have occured within one of the grid geometries
gdf_e_CA = gdf_e.loc[:, ["geometry", "mag"]].sjoin(gdf_grid)
# get median magnitude of eargquakes in grid
gdf_grid = gdf_grid.join(
    gdf_e_CA.dissolve(by="index_right", aggfunc="median").drop(columns="geometry")
)
# how many earthquakes in the grid
gdf_grid = gdf_grid.join(
    gdf_e_CA.dissolve(by="index_right", aggfunc=lambda d: len(d))
    .drop(columns="geometry")
    .rename(columns={"mag": "number"})
)

# drop grids geometries that have no measures and create folium map
m = gdf_grid.dropna().explore(column="mag")
# for good measure - boundary on map too
gdf_CA["geometry"].apply(lambda g: shapely.geometry.MultiLineString([p.exterior for p in g.geoms])).explore(m=m)

Answer 2

我想将 pandas DataFrame 转换为启用空间的 geopandas 之一：

df=pd.read_csv('../Desktop/test_esri.csv')
df.head()

然后使用转换：

gdf = geopandas.GeoDataFrame(
    df, geometry=geopandas.points_from_xy(df.long, df.lat))
from pyproj import crs
crs_epsg = crs.CRS.from_epsg(4326)
gdf=gdf.set_crs('epsg:4326')

然后我想将空间网格过度叠加为：

import numpy as np
import shapely
from pyproj import crs
# total area for the grid
xmin, ymin, xmax, ymax= gdf.total_bounds
# how many cells across and down
n_cells=30
cell_size = (xmax-xmin)/n_cells
# projection of the grid
# crs = "+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs"
# create the cells in a loop
grid_cells = []
for x0 in np.arange(xmin, xmax+cell_size, cell_size ):
    for y0 in np.arange(ymin, ymax+cell_size, cell_size):
        # bounds
        x1 = x0-cell_size
        y1 = y0+cell_size
        grid_cells.append( shapely.geometry.box(x0, y0, x1, y1)  )
cell = geopandas.GeoDataFrame(grid_cells, columns=['geometry'], 
                                 crs=crs.CRS('epsg:4326'))

然后将网格与 geopandas 数据框合并：

merged = geopandas.sjoin(gdf, cell, how='left', predicate='within')

最终在“溶解”中计算所需的指标：

# Compute stats per grid cell -- aggregate fires to grid cells with dissolve
dissolve = merged.dissolve(by="index_right", aggfunc="median")

但我认为我在“单元格”网格上做错了，我无法弄清楚！ 在here 找到使用的 csv 文件的提取。

Answer 3

感谢@Rob Raymond，
终于用下面的代码解决了：

import pandas as pd
import geopandas as gpd
import pyproj
import matplotlib.pyplot as plt
import numpy as np
import shapely
from folium import plugins

df=pd.read_csv('../Desktop/test_esri.csv')
gdf_monica = gpd.GeoDataFrame(
    df, geometry=gpd.points_from_xy(df.long, df.lat))
gdf_monica=gdf_monica.set_crs('epsg:4326')

gdf_area = gpd.read_file('https://raw.githubusercontent.com/openpolis/geojson-italy/master/geojson/limits_IT_municipalities.geojson')#.loc[:, ["geometry"]]

gdf_area =gdf_area[gdf_area['name']=='Portici'].loc[:,['geometry']]

BOXES = 50
a, b, c, d = gdf_area.total_bounds


gdf_grid = gpd.GeoDataFrame(
    geometry=[
        shapely.geometry.box(minx, miny, maxx, maxy)
        for minx, maxx in zip(np.linspace(a, c, BOXES), np.linspace(a, c, BOXES)[1:])
        for miny, maxy in zip(np.linspace(b, d, BOXES), np.linspace(b, d, BOXES)[1:])
    ],
    crs="epsg:4326",
)

# remove grid boxes created outside actual geometry
gdf_grid = gdf_grid.sjoin(gdf_area).drop(columns="index_right")


gdf_monica_binned = gdf_monica.loc[:, ["geometry", "CO"]].sjoin(gdf_grid)
# get median magnitude of CO pollutant
gdf_grid = gdf_grid.join(
    gdf_monica_binned.dissolve(by="index_right", aggfunc="median").drop(columns="geometry")
)
# how many earthquakes in the grid
gdf_grid = gdf_grid.join(
    gdf_monica_binned.dissolve(by="index_right", aggfunc=lambda d: len(d))
    .drop(columns="geometry")
    .rename(columns={"CO": "number"})
)

# drop grids geometries that have no measures and create folium map
m = gdf_grid.dropna().explore(column="CO")
# for good measure - boundary on map too
gdf_area["geometry"].apply(lambda g: shapely.geometry.MultiLineString([p.exterior for p in g.geoms])).explore(m=m)

产生：

如您所知，我对空间分析知之甚少或一无所知。如果不使用描述兴趣点所在几何的 geojson 数据，我无法获得正确的结果。 如果有人可以添加更多见解...谢谢！