将纬度/经度点转换为 GeoPandas 中的网格多边形

问题描述

我正在尝试将文本文件（按纬度、经度和花粉通量值组织）中的数据绘制为 Python 中的栅格网格。我正在使用 https://autogis-site.readthedocs.io/en/latest/notebooks/L5/02_interactive-map-folium.html 上的 Choropleth Map 的代码来尝试显示数据。我的 GeoPandas 地理数据框有点几何；然而，看起来教程中点的几何形状已经是多边形，我假设它是网格中的正方形。如何将我的数据（假设每个纬度/经度点是网格中像素的中心）转换为网格化的 geopandas（geodataframe）数据？我将使用的投影是 LAmbert Conformal Conic 投影。

为了澄清我的地理数据框的样子，在执行 gdf.head(10).to_dict() 时，它看起来像这样

    {'geoid': {0: '0',1: '1',2: '2',3: '3',4: '4',5: '5',6: '6',7: '7',8: '8',9: '9'},'geometry': {0: <shapely.geometry.point.Point at 0x7fa3e7feee90>,1: <shapely.geometry.point.Point at 0x7fa3e7Feed10>,2: <shapely.geometry.point.Point at 0x7fa3e7feef90>,3: <shapely.geometry.point.Point at 0x7fa3e7fe4f90>,4: <shapely.geometry.point.Point at 0x7fa3e7fe4e50>,5: <shapely.geometry.point.Point at 0x7fa3e7fe4bd0>,6: <shapely.geometry.point.Point at 0x7fa3e7fe4ed0>,7: <shapely.geometry.point.Point at 0x7fa3e7fe4c90>,8: <shapely.geometry.point.Point at 0x7fa3e7fe4d50>,9: <shapely.geometry.point.Point at 0x7fa3e7fe4c10>},'pollenflux': {0: 0.0,1: 0.0,2: 0.0,3: 0.0,4: 0.0,5: 0.0,6: 0.0,7: 0.0,8: 0.0,9: 0.0}}

什么时候应该这样格式化：

    {'geoid': {0: '0','geometry': {0: <shapely.geometry.multipolygon.Multipolygon at 0x7fa3e9363f50>,1: <shapely.geometry.multipolygon.Multipolygon at 0x7fa3e9363c90>,2: <shapely.geometry.multipolygon.Multipolygon at 0x7fa3e93631d0>,3: <shapely.geometry.multipolygon.Multipolygon at 0x7fa3e9363f10>,4: <shapely.geometry.multipolygon.Multipolygon at 0x7fa3e9363410>,5: <shapely.geometry.multipolygon.Multipolygon at 0x7fa3e9363a90>,6: <shapely.geometry.multipolygon.Multipolygon at 0x7fa3e9363d90>,7: <shapely.geometry.multipolygon.Multipolygon at 0x7fa3e9363d10>,8: <shapely.geometry.multipolygon.Multipolygon at 0x7fa3e9363390>,9: <shapely.geometry.multipolygon.Multipolygon at 0x7fa3e9363190>},'pop18': {0: 108,1: 273,2: 239,3: 202,4: 261,5: 236,6: 121,7: 196,8: 397,9: 230}}

解决方法

我相信您遇到的问题是代码需要带有多边形或多多边形的 GeoDataFrame，但您的只有点。

这里有一种快速的方法来生成一个新的 GeoDataFrame，在您的点周围有正方形：

import shapely
import numpy 

def get_square_around_point(point_geom,delta_size=0.0005):
    
    point_coords = np.array(point_geom.coords[0])

    c1 = point_coords + [-delta_size,-delta_size]
    c2 = point_coords + [-delta_size,+delta_size]
    c3 = point_coords + [+delta_size,+delta_size]
    c4 = point_coords + [+delta_size,-delta_size]
    
    square_geom = shapely.geometry.Polygon([c1,c2,c3,c4])
    
    return square_geom

def get_gdf_with_squares(gdf_with_points,delta_size=0.0005):
    gdf_squares = gdf_with_points.copy()
    gdf_squares['geometry'] = (gdf_with_points['geometry']
                               .apply(get_square_around_point,delta_size))
    
    return gdf_squares

# This last command actually executes the two functions above. 
gdf_squares = get_gdf_with_squares(gdf,delta_size=0.0005)

请注意，delta_size 参数表示正方形角点坐标与中心点之间的距离。当您的原始数据位于 WGS84 坐标 (EPSG 4326) 中时，如果您的数据位于德克萨斯州中部，则使用 0.0005 的增量将导致大约 100 米的平方。

查看您的输入数据，找到它正在使用的 CRS 并尝试估计一个好的 delta 值，该值将生成足够大的正方形但不会相互重叠。

希望能让其余代码正常工作。