GeoPandas入门实战:空间数据分析的Pandas式操作

大家好,我是小白鸽,本期为大家分享如何用GeoPandas进行空间数据分析。文章首先对GeoPandas进行了简要介绍,其次对GeoPandas的数据结构做了介绍,然后就GeoPandas的主要功能使用展示了代码示例,最后结合案例对GeoPandas功能做了补充介绍。

   GeoPandas简介

GeoPandas 是一个开源的Python 库,用于处理和分析地理空间数据。它建立在 Pandas 库之上,并添加了对几何形状和空间操作的支持。这使得 GeoPandas 成为处理地理空间数据的重要工具。

   GeoPandas 主要功能

  1. 地理数据类型: GeoPandas 支持多种几何数据类型,包括点、线、多边形和几何集合。

  2. 空间操作: GeoPandas 提供了丰富的空间操作函数,例如缓冲区生成、裁剪、合并、分割、拓扑运算等。

  3. 地图绘制: GeoPandas 可以与matplotlib 等绘图库集成,轻松创建地图和空间数据可视化。

  4. 数据分析: GeoPandas 可以进行各种空间数据分析,例如计算几何图形的面积、长度、质心等属性,以及进行空间统计分析。

  5. 数据导入和导出: GeoPandas 支持多种地理数据格式,例如 Shapefile、GeoJSON、GML、KML 等,方便进行数据交换和共享。

   GeoPandas数据结构

  • GeoSeries

    GeoSeries 是GeoPandas 的基本数据结构,它是一个类似 Pandas 的 Series,但是专门用于存储几何对象。每个元素都是一个几何形状,可以是点、线、多边形等。

主要属性和方法:

geometry: 返回一个包含所有几何对象的列表。

crs: 表示坐标参考系(Coordinate Reference System, CRS)的信息。

area: 计算每个几何对象的面积。

boundary: 返回每个几何对象的边界。

centroid: 返回每个几何对象的重心(质心)。

contains: 判断几何对象是否包含另一个几何对象。

distance: 计算两个几何对象之间的最短距离。

intersection: 计算两个几何对象的交集。

buffer: 为几何对象生成缓冲区。

创建 GeoSeries 示例:

import geopandas as gpd
from shapely.geometry import Point
# 创建一个包含点的GeoSeries
points = gpd.GeoSeries([Point(1, 1), Point(2, 2), Point(3, 3)])
print(points)
GeoDataFrame
  • GeoDataFrame

    GeoDataFrame 是GeoPandas 的另一个核心数据结构,它是 DataFrame 的扩展,专门用于处理包含几何数据的表格数据。GeoDataFrame 有一个特殊的列,通常命名为 geometry,用于存储几何对象。

主要属性和方法:

geometry: 特殊的列,用于存储几何数据。

crs: 表示整个GeoDataFrame 的坐标参考系。

to_crs: 将GeoDataFrame 的坐标参考系转换为新系统。

plot: 绘制GeoDataFrame 中的几何对象。

geom_type: 返回几何类型。

buffer: 为GeoDataFrame 中的所有几何对象生成缓冲区。

centroid: 返回GeoDataFrame 中所有几何对象的重心。

创建 GeoDataFrame 示例:

import geopandas as gpd
from shapely.geometry import Polygon
# 创建一个多边形
polygon = Polygon([(0, 0), (1, 0), (1, 1), (0, 1)])
# 创建一个 GeoDataFrame
gdf = gpd.GeoDataFrame({'geometry': [polygon]})
print(gdf)

   功能示例

  • 矢量数据读写
    import geopandas as gpd
    from shapely.geometry import Point
    # 生成虚拟POI数据(咖啡店)
    poi_data = {
        'name': ['Starbucks', 'Luckin Coffee', 'Costa'],
        'geometry': [Point(116.3, 39.9), Point(116.31, 39.91), Point(116.29, 39.88)]
    }
    poi_gdf = gpd.GeoDataFrame(poi_data, crs="EPSG:4326")  # WGS84坐标系
    # 生成虚拟地铁站数据
    metro_data = {
        'station': ['Station_A', 'Station_B'],
        'geometry': [Point(116.305, 39.905), Point(116.31, 39.895)]
    }
    metro_gdf = gpd.GeoDataFrame(metro_data, crs="EPSG:4326")
    # 保存为Shapefile和GeoJSON
    poi_gdf.to_file("poi_data.shp")  # 写入Shapefile
    metro_gdf.to_file("metro_data.geojson", driver='GeoJSON')  # 写入GeoJSON
    # 从文件读取数据
    poi = gpd.read_file("poi_data.shp")
    metro = gpd.read_file("metro_data.geojson")
  • 空间连接
    # 创建地铁站500米缓冲区(转换为投影坐标系计算距离)
    metro_projected = metro.to_crs(epsg=32650)  # UTM 50N
    metro_projected['buffer'] = metro_projected.geometry.buffer(500)  # 500米缓冲区
    # 将POI数据也转换到同一坐标系
    poi_projected = poi.to_crs(metro_projected.crs)
    # 空间连接:找出位于缓冲区内的POI
    joined = gpd.sjoin(
        poi_projected, 
        metro_projected[['station', 'buffer']].set_geometry('buffer'),
        how='inner', 
        predicate='within'
    )
    print("地铁站500米范围内的POI:")
    print(joined[['name', 'station']])
  • 缓冲区分析
    import matplotlib.pyplot as plt
    # 创建可视化画布
    fig, ax = plt.subplots(figsize=(10, 8))
    # 绘制原始地铁站
    metro.plot(ax=ax, color='red', markersize=100, label='Metro Stations')
    # 绘制缓冲区
    metro_projected.set_geometry('buffer').plot(
        ax=ax, 
        alpha=0.3, 
        edgecolor='blue', 
        facecolor='none',
        label='500m Buffer'
    )
    # 绘制POI点
    poi.plot(ax=ax, color='green', markersize=50, label='POI')
    # 添加图例和标题
    plt.legend()
    plt.title("Metro Station Buffers and POI Distribution")
    plt.show()
  • 坐标转换
    # 查看原始坐标系
    print("POI数据坐标系:", poi.crs)  # EPSG:4326 (WGS84)
    # 转换为Web墨卡托投影(EPSG:3857)
    poi_web_mercator = poi.to_crs(epsg=3857)
    print("转换后坐标系:", poi_web_mercator.crs)
    # 坐标系转换前后坐标示例
    print("\n原始坐标(WGS84):")
    print(poi.geometry.head())
    print("\n转换后坐标(Web墨卡托):")
    print(poi_web_mercator.geometry.head())

   案例分享

  • 背景

    假设现在有地铁站点数据和Poi数据,需要分析地铁站与Poi的关系,并且求出每个地铁站1km范围内的Poi数量。

  • 代码示例
    import contextily as ctx
    import matplotlib.pyplot as plt
    import geopandas as gpd
    from shapely.geometry import Point
    # 生成模拟数据(100个随机POI和10个地铁站)
    import numpy as np
    np.random.seed(42)
    # 生成随机POI(北京范围)
    poi_points = [Point(lon, lat) for lon, lat in 
                 zip(np.random.uniform(116.2, 116.4, 100),
                     np.random.uniform(39.8, 40.0, 100))]
    poi_gdf = gpd.GeoDataFrame(
        geometry=poi_points,
        data={'category': np.random.choice(['餐饮', '零售', '办公'], 100)},
        crs="EPSG:4326"
    )
    # 生成地铁站
    metro_points = [Point(lon, lat) for lon, lat in 
                   zip(np.random.uniform(116.25, 116.35, 10),
                       np.random.uniform(39.85, 39.95, 10))]
    metro_gdf = gpd.GeoDataFrame(
        geometry=metro_points,
        data={'station': [f'Station_{i}' for i in range(10)]},
        crs="EPSG:4326"
    )
    # 分析每个地铁站1km范围内的POI数量
    metro_proj = metro_gdf.to_crs(epsg=32650)
    metro_proj['buffer'] = metro_proj.geometry.buffer(1000)
    poi_proj = poi_gdf.to_crs(metro_proj.crs)
    joined = gpd.sjoin(
        poi_proj, 
        metro_proj.set_geometry('buffer'),
        how='inner',
        predicate='within'
    )
    # 统计结果
    result = joined.groupby('station').agg(
        total_poi=('category', 'count'),
        restaurant_count=('category', lambda x: (x=='餐饮').sum())
    ).reset_index()
    print("\n各地铁站周边POI统计:")
    print(result)
    # 可视化
    fig, ax = plt.subplots(figsize=(12, 10))
    # 绘制底图
    metro_proj.to_crs(epsg=3857).plot(
        ax=ax, 
        color='red', 
        markersize=50, 
        label='Metro Stations'
    )
    metro_proj.set_geometry('buffer').to_crs(epsg=3857).plot(
        ax=ax, 
        alpha=0.2, 
        color='blue'
    )
    poi_proj.to_crs(epsg=3857).plot(
        ax=ax, 
        column='category', 
        legend=True,
        markersize=20,
        categorical=True,
        cmap='tab20'
    )
    # 添加底图瓦片
    #OpenStreetMap的底图不一定能连接上,建议换自己的底图或者删除本行
    ctx.add_basemap(ax, source=ctx.providers.OpenStreetMap.Mapnik)
    plt.title("地铁站周边POI分布热力图")
    plt.axis('off')
    plt.show()
  • 结果

图片

   更多内容

Deepseek与规划协同:填平技术门槛,利用Deepseek进行空间句法分析

网页自动化神器|可视化爬虫|automa使用教程

【高德路径规划API】使用指南

用Deepseek做了个文献总结工具,让你写基金飞起来

更多内容可以关注微信公众号“小白鸽gis”

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值