python地理空间_GIS Experience (十):Python地理空间数据分析

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

前言

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

1 入门级

1.1 geopandas

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

'''读写数据'''

# 读取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()

#

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的范围

#

# print(shapes[1].parts)

# print(len(shapes)) # 输出要素数量

# print(file.numRecords) # 输出要素数量

# print(file.records()) # 输出所有属性表

#

'''

字段类型:此列索引处的数据类型。类型可以是:

“C”:字符,文字。

“N”:数字,带或不带小数。

“F”:浮动(与“N”相同)。

“L”:逻辑,表示布尔值True / False值。

“D”:日期。

“M”:备忘录,在GIS中没有意义,而是xbase规范的一部分。

'''

# fields = file.fields

# print(fields)

#

#

for index in range(len(shapes)):

geometry = shapes[index]

# print(geometry.shapeType)

# print(geometry.points)

#

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 参考资料

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值