python三方库GIS用法总结
注:工作学习中运用到的python三方库GIS用法总结
第一章 Python GIS学习之geopandas的使用
文章目录
前言
本文大概对工作学习中geopandas的运用进行阶段总结及代码段功能作用总结。
一、GeoPandas 是什么?
Geopandas简介
顾名思义,GeoPandas 通过添加对地理空间数据的支持来扩展流行的数据科学库 pandas。如果您不熟悉 ,我们建议您在继续操作之前快速浏览其 入门文档。
geopandas.GeoDataFrame; pandas.DataFrame; geopandas.GeoSeries; pandas.Series; GeoDataFrame; pandas.Series; geopandas.GeoSeries;
GeoPandas 中的核心数据结构是 GeoDataFrame,它是 pandas 的子类,可以存储几何列并执行空间操作。GeoSeries 是 GeoDataFrame 中处理几何对象的子类。因此,您的 GeoDataFrame 是 pandas DataFrame 的结合,具有传统数据(数值、布尔、文本等),以及具有几何对象(点、多边形等)的 GeoSeries。您可以拥有任意数量的带有几何对象的列;对于桌面GIS软件来说,通常没有限制。
GeoSeries; GeoSeries.crs; GeoSeries; GeoDataFrame;
每个数组可以包含任何几何类型(甚至可以在单个数组中混合它们),并且都有一个属性,用于存储有关投影的信息(CRS代表坐标参考系统)。因此,每个几何对象都可以处于不同的投影中,这使得你可以拥有同一几何对象的多个版本(不同的投影)。
二、安装方法
1. 下载相关依赖包
geopandas包在安装之前,需要手动依次将GDAL , Fiona, Pyproj,Rtree,Shapely五个依赖包安装成功;
根据已安装的python版本下载对应的依赖包;
下载地址为 https://www.lfd.uci.edu/~gohlke/pythonlibs/#wordcloud 。
2. 读入数据
包的安装顺序依次为GDAL->Fiona->pyproj->Rtree->Shapely->geopandas
pip install GDAL-3.4.2-cp38-cp38-win_amd64.whl
具体安装教程可根据这个大佬的文章:https://blog.csdn.net/Aries1Chan/article/details/124659444
三、用法总结
对geopandas部分实际功能作用进行总结,以便复用;
1. 读取shp文件
polygon_path = './data/districts.shp'
gdf= gpd.read_file(polygon_path)
2. 判断是否相交
gdf['geometry'].iloc[i].intersects(gdf['geometry'].iloc[j]) # 判断元素间是否相交
gdf['geometry'].iloc[i].intersection(gdf['geometry'].iloc[j]) # 求出相交的范围
3. 根据行政区划切分路网并计算各市的路网总长度
def innerShp(polyline, polygon, outputFile):
"""
计算每个行政区域的总长度
:param polyline: 路网数据
:param polygon: 行政区划数据
:param outputFile: 导出文件的路径
:return:
"""
# 使用空间关系进行切割
roads_by_admin = gpd.sjoin(polyline, polygon, how="inner", op="intersects")
# 投影坐标系转换,UTM 50N,求出的单位是m
roads_by_admin = roads_by_admin.to_crs(epsg=32650)
# 计算每个行政区域的总长度
# roads_by_admin["length"] = roads_by_admin["geometry"].agg(lambda x: x.length)
roads_by_admin["length"] = roads_by_admin["geometry"].length
text = roads_by_admin.groupby('市')["length"].sum().reset_index()
# 保存结果为excel
text.to_excel(outputFile, index=False)
4. 三维面转二维,同时转换坐标系
读取的polygon类型的shp中,面的构成是由内环外环两条线,三维转二维是把z值去掉,重新构建一个二维面
import geopandas as gpd
from shapely.geometry import Polygon, MultiPolygon
import transbigdata as tbd
# 读取三维Shapefile
shapefile_path = './data/area_3d_wgs84.shp'
gdf = gpd.read_file(shapefile_path)
def drop_z(geom):
if geom.has_z: # 如果几何包含Z坐标
if isinstance(geom, Polygon):
# 对Polygon的外环和内环去掉Z坐标,可能是镂空的面,即环状
return Polygon([(x, y) for x, y, z in geom.exterior.coords],
[[(x, y) for x, y, z in interior.coords] for interior in geom.interiors])
elif isinstance(geom, MultiPolygon):
# 对MultiPolygon的每个Polygon都去掉Z坐标
return MultiPolygon([drop_z(poly) for poly in geom.geoms])
return geom # 如果没有Z坐标,直接返回
# 将所有几何从三维转换为二维
gdf['geometry'] = gdf['geometry'].apply(drop_z)
# 转换坐标系到gcj02
gdf = tbd.transform_shape(gdf, tbd.wgs84togcj02)
gdf.to_file('./output/area_2d_gcj02.shp')
5.
逐步更新中……