GIS Experience （十）：Python地理空间数据分析

目录

• 前言
• 1 入门级
  • 1.1 geopandas
  • 1.2 pyshp
  • 1.4 GDAL & OGR
  • 1.5 pysal
    • 1.5.1 explore
    • 1.5.2 viz
    • 1.5.3 model
    • 1.5.4 lib
• 2 进阶级
  • 2.1 矢量创建
    • 2.1.1 点转线
• 3 参考资料

备注：虽然在Pycharm中借助SciView可以很好地进行分屏显示，但地图数据一般数据量较大，所以用python进行地图可视化需要先行关闭。

前言

GISer都知道在常用的桌面端GIS应用ArcGIS和QGIS工具中都大量的使用了Python语言，考虑到当前python也被大量应用到机器学习和人工智能领域，在进行空间处理时候直接通过编写代码也可以使得工作更为高效。

空间处理库	链接地址	描述
geopandas	https://pypi.org/project/geopandas/#files	扩展pandas对地理数据的支持
pyproj	https://pypi.org/project/pyproj/#files	坐标转换
Shapely	https://pypi.org/project/Shapely/#files	2D地理对象操作及分析
Fiona	https://pypi.org/project/Fiona/#files	空间数据读写
pyshp	https://pypi.org/project/pyshp/#files	ESRI系列空间数据读写
GDAL	https://pypi.org/project/GDAL/#files	gdal用于复杂地理栅格数据处理，ogr用于复杂地理矢量数据处理
pysal	https://pypi.org/project/pyshp/#files	空间分析
geoplot	https://pypi.org/project/geoplot/	空间可视化
geopy	https://pypi.org/project/geopy/	网页地图服务-空间定位

1 入门级

1.1 geopandas

GeoPandas主要目的是使得在Python环境下更方便的处理地理空间数据，其扩展了pandas的数据类型，允许其在几何类型上进行空间操作。几何操作由 shapely执行，fiona进行文件存取，并借助descartes和matplotlib 进行绘图，详见geopandas官方文档。

备注：pip install geopandas时候可能存在无法安装的情况，可以前往国内镜像网站直接下载并安装对应的whl文件。
1）http://mirrors.aliyun.com/pypi/simple/ 阿里云
2）https://pypi.mirrors.ustc.edu.cn/simple/ 中国科技大学
3）https://pypi.tuna.tsinghua.edu.cn/simple/ 清华大学
4）http://pypi.mirrors.ustc.edu.cn/simple/ 中国科学技术大学
5）http://pypi.sdutlinux.org/ 山东理工大学
6）http://pypi.hustunique.com/ 华中理工大学
7）https://www.lfd.uci.edu/~gohlke/pythonlibs/ 加利福利亚大学
8）http://pypi.douban.com/ 豆瓣

'''读写数据'''
# 读取shp数据
df = gpd.GeoDataFrame.from_file("Point.shp")

# 写入为shp数据
df.to_file(filename="New_Point.shp", driver="ESRI Shapefile", schema=None)
# 写入为geojson数据
df.to_file("Point.geojson", driver='GeoJSON')
# 写入为geopackage数据
df.to_file("Point.gpkg", layer='Point', driver="GPKG")

'''坐标转换'''
# way1 proj4写法，详见[spatialreference](https://spatialreference.org/)
df = df.to_crs(crs="+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext  +no_defs")

# way2 EPSGcode，详见[epsg](http://epsg.io/)
df = df.to_crs({'init': 'epsg:3857'})

1.2 pyshp

fiona与pyshp都可以读写ESRI shapefiles，但pyshp支持的效果更好，相比geopandas将几何地理对象读写为数据框，pyshp是直接读写为几何对象，故pyshp进行坐标和几何信息提取也稍显便捷。

# 导入shapely
from shapely import wkt, geometry
 
wktPoly = "POLYGON((0 0, 4 0, 4 4, 0 4, 0 0))"
poly = wkt.loads(wktPoly)
 
# 打印多边形的面积
print(poly.area)
 
# 建立缓冲区并计算面积
buf = poly.buffer(5.0)
print(buf.area)
 
# 计算面积差异
print(buf.difference(poly).area)

import shapefile  # 使用pyshp库
 
file = shapefile.Reader("data\市界.shp")
shapes = file.shapes()
 
# <editor-fold desc="读取元数据">
print(file.shapeType)  # 输出shp类型
'''
NULL = 0
POINT = 1
POLYLINE = 3
POLYGON = 5
MULTIPOINT = 8
POINTZ = 11
POLYLINEZ = 13
POLYGONZ = 15
MULTIPOINTZ = 18
POINTM = 21
POLYLINEM = 23
POLYGONM = 25
MULTIPOINTM = 28
MULTIPATCH = 31
'''
print(file.bbox)  # 输出shp的范围
# </editor-fold>
# print(shapes[1].parts)
# print(len(shapes))  # 输出要素数量
# print(file.numRecords)  # 输出要素数量
# print(file.records())  # 输出所有属性表
 
# <editor-fold desc="输出字段名称和字段类型">
'''
字段类型：此列索引处的数据类型。类型可以是：
“C”：字符，文字。
“N”：数字，带或不带小数。
“F”：浮动（与“N”相同）。
“L”：逻辑，表示布尔值True / False值。
“D”：日期。
“M”：备忘录，在GIS中没有意义，而是xbase规范的一部分。
'''
# fields = file.fields
# print(fields)
# </editor-fold>
 
 
# <editor-fold desc="输出几何信息">
for index in range(len(shapes)):
    geometry = shapes[index]
    # print(geometry.shapeType)
    # print(geometry.points)
# </editor-fold>

1.4 GDAL & OGR

GDAL和OGR是开源空间信息基金会(Open Source Geospatial Foundation，简称OSGeo)推出的开源空间处理库，两者分别被应用于栅格和矢量数据处理。

备注：安装QGIS后本地会存在GDAL和OGR，可以将其复制到pycharm或者是Anaconda第三方库文件，此时包引入规则from osgeo import gdal, ogr

1.5 pysal

Pysal是一个面向地理空间数据科学的开源跨平台库，重点是用python编写的地理空间矢量数据。它支持空间分析高级应用程序的开发（详见官方文档或者osgeo译制文档）。此外，考虑下载速度原因，采用国内清华镜像下载安装pip install -i https://pypi.tuna.tsinghua.edu.cn/simple pysal，中途若遇到安装Rtree时提示 OSError: could not find or load spatialindex_c-64.dll，请自行前往加利福利亚大学镜像站获取whl文件。

1.5.1 explore

用于对空间和时空数据进行探索性分析的模块，包括对点、网络和多边形格的统计测试。还包括空间不等式和分布动力学的方法。

1.5.2 viz

可视化空间数据中的模式，以检测集群、异常值和热点。

1.5.3 model

用各种线性、广义线性、广义加性和非线性模型对数据中的空间关系进行建模。

1.5.4 lib

解决各种各样的计算几何问题：

• 从多边形格、线和点构建图形。
• 空间权重矩阵与图形的构建与交互编辑
• α形状、空间指数和空间拓扑关系的计算
• 读写稀疏图形数据，以及纯python空间矢量数据阅读器。

2 进阶级

2.1 矢量创建

2.1.1 点转线

1）由点生成闭合线

    def PointToClosedLine():
    	# 读取shp 地图格式数据，读取后的数据结构为DataFrame格式
    	df = gpd.GeoDataFrame.from_file("折点.shp")
        # 按原始面分组
        dataGroup = df.groupby('ORIG_FID')
        # 打印分组单元的前五行记录
        # print(dataGroup.head(5))

        # 构造数据
        tb = []
        geomList = []
        for name, group in dataGroup:
            # 分离出属性信息，取每组的第1行前8列作为数据属性
            tb.append(group.iloc[0, :8])
            # 借助列表推导式，把同一组的点打包到一个list中
            xyList = [xy for xy in zip(group.POINT_X, group.POINT_Y)]
            # 取同一组所有点生成闭合图形
            line = LineString(xyList)
            geomList.append(line)
        # print(tb)
        # 点转线
        geoLine_DataFrame = gpd.GeoDataFrame(tb, geometry=geomList)
        geoLine_DataFrame.plot()
        plot,show()

2）由点生成隔离线段

    def PointToIsolatediLine():
    	# 读取shp 地图格式数据，读取后的数据结构为DataFrame格式
    	df = gpd.GeoDataFrame.from_file("折点.shp")
        # 按原始面分组
        dataGroup = df.groupby('ORIG_FID')
        # 打印分组单元的前五行记录
        # print(dataGroup.head(5))

        # 构造数据
        tb = []
        geomList = []
        for name, group in dataGroup:
            # 分离出属性信息，取每组的第1行前8列作为数据属性
            tb.append(group.iloc[0, :8])
            # 借助列表推导式，把同一组的点打包到一个list中
            xyList = [xy for xy in zip(group.POINT_X, group.POINT_Y)]
            for near in range(len(xyList)-1):
                # 分离出属性信息，取每组的第near行前8列作为数据属性
                tb.append(group.iloc[near, :8])
                # 取同一组所有点生成分离图形
                line = LineString([xyList[near], xyList[near+1]])
                geomList.append(line)
        # print(tb)
        # 点转线
        geoLine_DataFrame = gpd.GeoDataFrame(tb, geometry=geomList)
        geoLine_DataFrame.plot()
        plot,show()

3 参考资料

1. GeoPandas 0.6.0