地理空间分析中常用的python模块

一、背景

本文内容主要参照《python地理空间分析》，同时加入一些自己的理解。python是一门非常强大的编程语言。对于遥感图像处理来说，其计算量非常大，例如一景Landsat 8 OLI影像大约1G大小，影像的每一个波段约100MB，幸运的是，python对于遥感图像处理提供了很多高效的解决方案。下面对一些在遥感图像处理中常用的python模块进行介绍。

二、运行环境

所有的程序都经过了我的测试，我的运行环境为：

anoconda3 （python3.6）

pycharm

win10

三、文件读写操作库

1.读写shapefile文件的模块——pyshp

pyshp是一个完全用python实现的软件包，其主要功能是用于读写shapefile文件。它不支持任何几何操作，只调用了python标准库，并且能够包含在一个单独的文件中和小型的嵌入式平台集成。

pyshp可以直接在pycharm中下载安装：

下面给出一个读写shapefile文件的示例：

# 导入pyshp
import shapefile

shp = shapefile.Reader("DEM.shp")
for feature in shp.shapeRecords():
    point = feature.shape.points[0]
    rec = feature.record[0]
    print(point[0], point[1], rec)

2.通用矢量库——OGR

OGR是一种通用矢量库，除了可以像pyshp一样读写shp，dbf等矢量文件外，还可以处理WKT文本字符串（WKT是一种表达集合图形和空间索引系统的一种简单文本格式，目前读写WKT的最佳方式是shapely库）。

anoconda中自带于osgeo库，ogr可以直接从osgeo中导入。

下面以读取坐标点并打印为例，说明OGR的简单使用：

# 导入OGR库
from osgeo import ogr

shp = ogr.Open("point.shp")
layer = shp.GetLayer()
for feature in layer:
    geometry = feature.GetGeometryRef()    
    print(geometry.GetX(), geometry.GetY(),
          feature.GetFiled("FIRST_FLD"))

3.通用几何库——shapely

shapely的主要用途是通用几何库，它是提供python风格GEOS库几何操作的库。事实上，shapely特意避免访问文件，它完全依赖其他模块导入或导出数据，功能聚焦于几何操作。

shapely库可以直接从pycharm中搜索下载。

下面以计算集合面积为例，说明shaply的简单使用：

# 导入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)

4.栅格数据处理——GDAL

GDAL库是处理栅格数据的主流地理空间库，在python中使用GDAL读写地理空间栅格数据非常高效。使用一些额外的小技巧能很好的执行影像的重投影操作。但是其真正价值在于能够和python的其他模块进行交互。

在anoconda中，GDAL库集成在osgeo中，可以直接导入。

以一个tif影像的读取操作为例，说明GDAL的简单使用：

# 导入GDAL
from osgeo import gdal

raster = gdal.Open("RSimage.tif")

print(raster.RasterCount)
print(raster.RasterXSize)
print(raster.RasterYsize)

5.图像编辑——PIL

PIL库本来是用于处理遥感图像，但目前在python中主要用于图像编辑，为了保证速度，PIL使用C语言实现，而且专门针对python做了一些优化。此外在图片创建和处理方面，它包含一个很有用的上个绘制模块。

anoconda中自带有PIL图像库，可以直接导入。

下面将会结合pyshp库和PIL库对shapefile文件进行栅格化，并将其保存为一张图片。示例代码如下：

# 导入PIL
from PIL import Image
from PIL import ImageDraw

# 导入pyshp
import shapefile

r = shapefile.Reader("DEM.shp")

# 确定缩放比例
xdist = r.bbox[2] - r.bbox[0]
ydist = r.bbox[3] - r.bbox[1]
iwidth = 400
iheight = 600
xratio = iwidth / xdist
yratio = iheight / ydist
pixels = []
for x,y in r.shape()[0].points:
    px = int(iwidth - (r.bbox[2] - x) * xratio)
    py = int((r.bbox[3] - y) * yratio)
    pixels.append((px, py))

img = Image.new("RGB", (iwidth, iheight), "white")
draw = ImageDraw.Draw(img)
draw.ploygon(pixels, outline="rgb(203, 196, 190)", fill="rgb(198, 204, 189)")
img.save("DEM.png")

6.海量数据表格处理——Geopandas

Geopandas是由shapely，Fiona，pyproj，matplotlib以及其他必须的库一同构建的pandas的地理空间扩展，它能够处理海里数据表格，有序化/无序化、标签化矩阵以及无标记的统计数据。当然它需要特定的数据库支持。

anoconda中带有Geopandas，也可以从pycharm中安装该库。

下面的示例功能是打开一个shapefile文件，将其转化为GeoJSON格式，然后使用matplotlib创建一张地图。代码如下：

# 导入geopandas
import geopandas
import matplotlib as plt

gdf = geopandas.GeoDataFrame
census = gdf.from_file("polygon.shp")
census.plot()
plt.show()

7.多维数组处理工具——numpy

numpy是基于C语言实现的，专门为科学计算设计的一款高效的，多维python数组处理工具。numpy可以和GDAL，shapely， PIL库进行数据交互。

anoconda中自带numpy，可以直接导入。

numpy是python中最常用的库之一，因此本例集合GDAL读取遥感图像，并获取第一个波段创建直方图，将结果保存为jpg。GDAL与numpy之间的接口是GDAL库中的gdal_array模块，该模块会讲numpy当作一个引用依赖。代码为：

from osgeo import gdal_array

srcArray = gdal_array.LoadFile("image.tif")
band1 = arcArray[0]
gdal_array.SaveArray(band1, "band1.jpg", format="JPEG")

四、文件下载操作库

1.通过URL访问文件——urllib

urllib库的主要用途是通过URL地址来访问任意文件。在urllib库中，调用open()方法能够创建一个资源的链接，但是不访问任何数据。调用urllib.request.retrieve()方法能够找到相应文件并下载到硬盘上。

anoconda中自带urllib库，可以直接导入。

例如下载USGS上的地震资料CSV文件，代码为：

import urllib.request
import urllib.error

url = "http://earthquake.usgs.gov/earthquakes/feed/v1.0/summary/all_hour.csv"

earthquake = urllib.request.urlopen(url)

# 读取CSV的单行数据
earthquake.readline()

# 逐行读取CSV数据
for record in earthquake: print(record)

2.FTP客户端下载数据——FTP

2004年之前，空间web服务流行起来之前，FTP曾经是最常用的发布地理空间数据的方式之一。ftplib能够通过FTP下载地理空间数据。

acoconda中带有ftplib，可直接导入使用。

以下示例访问NOAA中FTP中的一个文件，文件中包含监测世界各地海啸的网络浮标数据。代码如下：

import ftplib

server = "ftp.ngdc.noaa.gov"
dir = "hazards/DART/20070815_peru"
filename = "21415_from_20070727_08_55_15_tides.txt"

ftp = ftplib.FTP(server)
ftp.login()
ftp.cwd(dir)

out = open(filename, "wb")
ftp.retrbinary("RETR" + filename, out.write)
out.close()

dart = open(filename)
for line in dart:
    if "LAT," in line:
        print(line)
        break

3.解压缩文件——ZIP和TAR

很多地理空间数据问价是zip或者tar的压缩格式文件，python中专门读写ZIP和TAR格式的模块分别为zipfile和tarfile。

anoconda中带有zipfile和tarfile，可以直接导入使用。

下面以解压文件为例，说明zip和tar文件解压的简单操作：

# 解压zip文件
import zipfile

zip = open("image.zip", "rb")
zipimage = zipfile.ZipFile(zip)
for filename in zipimage.namelist():
    out = open(filename, "wb")
    out.write(zipimage.read(filename))
    out.close()

# 解压tar文件
import tarfile

tar = tarfile.open("image.tar.gz", "r:gz")
tar.extractall()
tar.close()

五、其他一些非官方库

这里的非官方库是指无法从pycharm中直接搜索到，需要自己手动下载安装的库，下面会介绍这些模块并给出其相应的下载地址。

1.专门处理dbf文件——dbfpy3

OGR库和pyshp库都能读写dbf文件，因为dbf文件都是shapefile文件规范的一部分。dbf文件包含shapefile文件的属性和字段信息，但是上面的两个库都只提供了对dbf文件的基本操作。若需要对dbf文件进行更复杂的操作，需要自己手动安装dbfpy3库。

dbfpy3不能从pycharm上直接下载，不过可以使用pip安装，这里给出该库的下载地址：

https://pypi.org/project/dbfpy/2.3.1/

2.轻量级的PIL——PNGCanvas

如果因为机器上没有安装python模块和使用C编译器的管理员权限而无法安装PIL，此时可以考虑安装PNGCanvas。该模块完全使用python实现，也具有对栅格数据处理的能力，是轻量级的PIL。

下面给出PNGCanvsa的下载地址：

https://pypi.org/project/pngcanvas/

3.创建PDF格式的地图——pyFPDF

pyFPDF库是完全使用python代码实现的，它是一种轻量级的创建包括地图等PDF文件的解决方案。

下面给出pyFPDF的下载地址：

https://pypi.org/project/fpdf2/

4.高光谱数据处理——SPy

SPy(Spectral Python)库是专门处理遥感应用的高级功能包，它在高光谱数据处理方面的功能非常强大。如果安装该库后并添加相关的依赖包，则SPy能够提供完全窗口化的处理环境。

SPy的官方地址（我这里需要翻墙才能打开）为：

http://www.spectralpython.net/