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