Geopandas

最新推荐文章于 2025-11-14 16:21:16 发布

原创最新推荐文章于 2025-11-14 16:21:16 发布 · 1.9k 阅读

2 ·

CC 4.0 BY-SA版权

文章标签：

#python

Z类分类_python包的使用专栏收录该内容

23 篇文章

订阅专栏

本文介绍如何使用GeoPandas和Shapely进行空间数据操作与可视化，包括多边形创建、空间叠加分析、图形转换及画图技巧。通过示例展示了矩形多边形的交集计算，并提供了将多边形转换为列表的方法，以及如何绘制不同类型的边界图形。

from shapely.geometry import Polygon, mapping
import geopandas as gpd
from matplotlib import pyplot as plt
import  pandas as pd


polys1 = gpd.GeoSeries([Polygon([(0,0), (2,0), (2,2), (0,2)]),
         Polygon([(2,2), (4,2), (4,4), (2,4)])])

polys2 = gpd.GeoSeries([Polygon([(1,1), (3,1), (3,3), (1,3)]),
         Polygon([(3,3), (5,3), (5,5), (3,5)])])

df1 = gpd.GeoDataFrame({'geometry': polys1, 'df1':[1,2]})
df2 = gpd.GeoDataFrame({'geometry': polys2, 'df2':[1,2]})

res_union = gpd.overlay(df1, df2, how='intersection')
# >>> res_union
# >>>    df1  df2                             geometry
# >>> 0    1    1  POLYGON ((1 2, 2 2, 2 1, 1 1, 1 2))
# >>> 1    2    1  POLYGON ((2 2, 2 3, 3 3, 3 2, 2 2))
# >>> 2    2    2  POLYGON ((3 4, 4 4, 4 3, 3 3, 3 4))

x,y = res_union['geometry'][0].exterior.coords.xy
# >>>x,y
# >>>(array('d', [1.0, 2.0, 2.0, 1.0, 1.0]), array('d', [2.0, 2.0, 1.0, 1.0, 2.0]))
coordinates = mapping(res_union['geometry'][0])['coordinates'][0]
# >>> ((1.0, 2.0), (2.0, 2.0), (2.0, 1.0), (1.0, 1.0), (1.0, 2.0))

ax = df1.plot(color='red')
df2.plot(ax=ax, color='green', alpha=0.5)
plt.show()

res_union.plot(ax=ax,color='yellow',alpha=0.5)
plt.show()

Geopandas 依赖库
从 http://www.lfd.uci.edu/~gohlke/pythonlibs 下载 Fiona , GDAl , pyproj , Shapely,Rtree，descartes

这是两组矩形图，重叠部分可以从颜色上看出来
在这里插入图片描述
黄色是后续计算出重叠部分，然后再重新画上颜色

20200327增加

将polygon转成list的小功能

def polygon2list(polys):
    '''
    shapely空间函数数据点提取
    :param polys: shapely.geometry.polygon.Polygon
    :return: list
    '''
    temp_str = polys.to_wkt()
    temp_str = temp_str.replace('POLYGON ','').replace('(','').replace(')','')
    temp_str = temp_str.split(',')
    temp_str = [i.strip().split(' ') for i in temp_str]
    temp_str = [[round(float(i[0]),6),round(float(i[1]),6)] for i in temp_str]
    return temp_str

额外的画图

    # 合并画图
    polys = gpd.GeoSeries(source['polygon'].tolist())
    b = polys.boundary
    c = polys.convex_hull
    d = gpd.GeoSeries(polys.unary_union)
    e = gpd.GeoSeries(d).boundary
    b.plot(alpha=0.5)
    plt.savefig(f'pic/all_boundary.jpg')
    c.plot(alpha=0.5)
    plt.savefig(f'pic/all_convex_hull.jpg')
    d.plot(alpha=0.5)
    plt.savefig(f'pic/all_unary_union.jpg')
    e.plot(alpha=0.5)
    plt.savefig(f'pic/all_unary_union_boundary.jpg')

boundary
在这里插入图片描述
convex_hull

unary_union

unary_union_boundary

    # 20200408补充分块画图
    rssult_df['border_polygon'] = rssult_df['border'].apply(lambda x:Polygon(x))
    polys = gpd.GeoSeries(rssult_df['border_polygon'].tolist())
    b = polys.boundary
    c = polys.convex_hull
    d = gpd.GeoSeries(polys.unary_union)
    e = gpd.GeoSeries(d).boundary

    b.plot(alpha=0.5,figsize=(30,20))
    plt.savefig(f'pic/all_boundary_1.jpg')

在这里插入图片描述

20200706增加

使用得是shapely来做计算

from shapely.geometry import Point,MultiPoint,Polygon
def createArea(lon_lat:list,area_tyoe = 'rectangle', radius = 0.009):
    '''
    创建区域(标准方形还是圆形)
    :param lon_lat: 经纬度
    :param area_tyoe: 创建类型:rectangle 方形 roundness 圆形
    :param radius: 方形为边长1半，圆形为半径
    :return: 图形边界数据
    '''
    if area_tyoe == 'rectangle':
        return Polygon([[lon_lat[0] + radius, lon_lat[1] + radius],
                        [lon_lat[0] + radius, lon_lat[1] - radius],
                        [lon_lat[0] - radius, lon_lat[1] - radius],
                        [lon_lat[0] - radius, lon_lat[1] + radius],
                        [lon_lat[0] + radius, lon_lat[1] + radius]]
                   )
    if area_tyoe == 'roundness':
        return Point(lon_lat).buffer(radius)

def roundness(points):
    '''
    最小外圆形:
    首先计算最小外凸多边形，然后计算质心点
    然后计算质心点与最小外凸多边形各个顶点
    的最远距离，即为对应的最小外围圆范围
    :param points: shapely.geometry.multipoint.MultiPoint
    :return: center_point:list, roundness:shapely.geometry.polygon.Polygon
    '''
    polys = points.convex_hull
    border = polygon2list(polys)  # 最小外凸多边形
    center_point = polys.centroid.to_wkt().replace('POINT (','').replace(')','').split(' ')   # 质心点
    center_point = [round(float(dot),6)for dot in center_point]
    buffer_length = [ Point(center_point).distance(Point(point)) for point in border]
    buffer_length = round(max(buffer_length),6)
    return center_point,Point(center_point).buffer(buffer_length)

def drawPic(points,polys,file_path):
    '''
    画图，作为开发中参考画出的图形的实际效果
    :param points: 点
    :param polys: 图形边界
    :param file_path: 存储文件位置
    :return: 无
    '''
    from descartes.patch import PolygonPatch
    from matplotlib import pyplot as plt
    plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签
    plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号
    fig = plt.figure()  # 设定图片大小
    ax = fig.add_subplot()
    patch = PolygonPatch(polys, facecolor='#6495ED', edgecolor='#6495ED', alpha=0.5, zorder=2)
    ax.add_patch(patch)
    for point in points:
        ax.plot(point.x, point.y, 'o', color='GRAY')
    # plt.show()
    plt.savefig(file_path)
    plt.close()


# **************************边界计算****************************
# 计算每个聚落的最小外凸多边形，最小外矩形，最小外圆形
splace_data_border = pd.DataFrame()
for label in splace_data_lables:
    # label = 60
    temp = splace_data[splace_data['label'] == label]
    points = MultiPoint(temp[['longitude', 'latitude']].values)
    # 存在数据质量问题，需要判断全列是否都是一个数
    if len(temp['longitude'].unique()) * len(temp['latitude'].unique()) == 1:
        # 当数据全一样时，直接以数据点作为原点，100米作为半径
        lon_lat = [list(temp['longitude'].unique())[0], list(temp['latitude'].unique())[0]]
        polys = createArea(lon_lat, area_tyoe='roundness')
        border = polygon2list(polys)
        center_point = lon_lat
    # 存在数据质量问题，需要判断全列是否都是成一条直线
    elif 'LINESTRING' in points.convex_hull.to_wkt():
        # 当数据是一条直线时，取其中点作为原点，其中点到一端得距离为半径
        temp_str = points.convex_hull.to_wkt().replace('LINESTRING ','').replace('(','').replace(')','')
        temp_str = temp_str.split(',')
        temp_str = [i.strip().split(' ') for i in temp_str]
        lon_lat = [(float(temp_str[0][0]) + float(temp_str[1][0]))/2,
                   (float(temp_str[0][1]) + float(temp_str[1][1]))/2]
        radius = round(Point(lon_lat).distance(Point([float(a) for a in temp_str[0]])),6)
        # 20200629添加对于半径过小导致精度不满足，计算结果为0情况修正
        if radius == 0:
            radius = 0.000009
        polys = createArea(lon_lat, area_tyoe='roundness',radius = radius)
        border = polygon2list(polys)
        center_point = lon_lat
    else:
        # 画图可以注释掉
        polys = points.convex_hull
        # drawPic(points, polys, f'result/pic/{c_class_big_id}_{year}_{month}_{label}_convex_hull.jpg')
        border_3 = polygon2list(polys)  # 最小外凸多边形
        polys = points.envelope
        # drawPic(points, polys, f'result/pic/{c_class_big_id}_{year}_{month}_{label}_envelope.jpg')
        border_1 = polygon2list(polys)  # 最小外矩形(边与坐标轴平行)
        polys = points.minimum_rotated_rectangle
        # drawPic(points, polys, f'result/pic/{c_class_big_id}_{year}_{month}_{label}_minimum_rotated_rectangle.jpg')
        border_2 = polygon2list(polys)  # 最小外矩形(边不限于与坐标轴平行)
        polys = roundness(points)
        # drawPic(points, polys, f'result/pic/{c_class_big_id}_{year}_{month}_{label}_roundness.jpg')
        border = polygon2list(polys)  # 最小外圆形

在这里插入图片描述