基于Python的地理空间分析（四）：矢量数据分析

发表于 分类于 阅读次数： Disqus：

引言

本文主要介绍矢量数据的处理和地理空间分析，主要用到Shapely、OGR和GeoPandas三个库，展示如何用这三个库执行基本的地理空间分析，以及更加复杂的数据科学任务，主要内容包括：

  • 读写、创建、操作矢量数据
  • 矢量数据的地图可视化
  • 使用地图投影、执行空间操作
  • 矢量几何和属性数据分析

OGR简单特征库

OGR简单特征库（OGR Simple Features Library）是GDAL的一部分，提供了一套处理矢量数据的工具。尽管OGR（矢量部分）和GDA（山各部分）L已经集成在一起，但是两者在使用上还是比较清晰的。虽然OGR是C编写的（文档也是C的），但是Python bindings使得我们可以通过Python访问所有GDAL的功能。接下来，介绍三部分内容：

  • OGR批处理命令与Python脚本
  • OGR库与Python bindings

OGR批处理命令

OGR提供了一系列的批处理命令，包括我们前面已经用过的ogrinfo和ogr2ogr命令，这些命令支持SQL查询的语法。此外，ogrindex和ogr2vrt可以用来创建矢量数据文件，ogrmerge能够将多个数据文件合并，并且创建指定格式的数据文件。

GDAL在安装的过程中，已经提供了一组用于地理空间分析的Python脚本，这些脚本的路径类似于：C:\ProgramData\Anaconda3\pkgs\gdal-2.4.0-py37hdf5ee75_1002\Scripts ，可以直接从Jupyter Notebook（使用魔法命令%run）或者终端运行。不过，从这个文件夹中的Python脚本列表中可以看到，几乎所有这些脚本都是针对GDAL而不是OGR。
需要注意的是，如果要运行这些脚本，需要正确引用这些脚本的路径，或者添加到系统变量当中。

OGR库与Python bindings

OGR库与其Python绑定相结合，构成了在Python中处理矢量数据的最重要部分。利用OGR库，我们可以创建点，线和多边形，并对这些元素执行空间操作。 例如，您可以计算几何体的面积，将不同的数据叠加在一起，并使用接近工具（如缓冲区）。 此外，与ogrinfo和ogr2ogr一样，OGR库提供了读取矢量数据文件，迭代单个元素以及选择和重新投影数据的工具。

OGR主要模块和类

OGR主要包括两个模块：ogr模块主要处理矢量几何，和osr模块，这两个模块都是osgeo的子模块。ogr模块主要处理矢量几何，osr则是处理投影。OGR库提供了以下七个类：

  • Geometry
  • Spatial Reference
  • Feature
  • Feature Class Definition
  • Layer
  • Dataset
  • Drivers

这些库的名字直观地说明了类的功能以及OGR的结构。OGR的模块、类和函数能够在GDAL网站(www.gdal.org/python )上找到文档，但是没有Python的代码示例。好在其他的Python库例如Fiona、GeoPandas提供了更加用户友好的方式。

创建几何多边形

我们可以用OGR编写点、线、多边形等矢量几何，并且可以通过投影给出坐标或者米等单位。并且可以用well-known binary(WKB)进行编码，并将最终的段变形转换为well-known text(WKT)：

1
2
3
4
5
6
7
8
9
10
11
from osgeo import ogr

r= ogr.Geometry(ogr.wkbLinearRing)
r.AddPoint(1,1)
r.AddPoint(5,1)
r.AddPoint(5,5)
r.AddPoint(1,5)
r.AddPoint(1,1)
poly = ogr.Geometry(ogr.wkbPolygon)
poly .AddGeometry(r)
print(poly.ExportToWkt())
POLYGON ((1 1 0,5 1 0,5 5 0,1 5 0,1 1 0))

从GeoJSON创建多边形

除了直接创建多边形以外，我们还可以通过传入GeoJSON创建，具体代码如下：

1
2
3
4
5
6
geojson = """{"type":"Polygon","coordinates":
[[[1,1],[5,1],[5,5],[1,5], [1,1]]]}"""
# geojson本身的易读性有点差

polygon = ogr.CreateGeometryFromJson(geojson)
print(polygon)
POLYGON ((1 1,5 1,5 5,1 5,1 1))

基本几何操作

下面是一些对多边形的基本操作，包括创建区域，获得中心、边界、凸包、缓冲（buffer），以及判断是否有某个点：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
# 1. 创建区域
print(f'The area of our polygon is {polygon.Area()}')

# 2. 计算中心点
print(f'The centroid of our polygon is {polygon.Centroid()}')

# 3. 获得边界
print(f'The boundary of our polygon is {polygon.GetBoundary()}')

# 4. 获得凸包
print(f'As our polygon is a square, the convex hull is the same:{polygon.ConvexHull()}')

# 5. 获得缓冲
print(f'Set a buffer of 0 for our polygon:{polygon.Buffer(0)}, which is the same.')

# 6. 判断点是否在多边形内
point = ogr.Geometry(ogr.wkbPoint)
point.AddPoint(6, 6)
print(polygon.Contains(point))
The area of our polygon is 16.0
The centroid of our polygon is POINT (3 3)
The boundary of our polygon is LINESTRING (1 1,5 1,5 5,1 5,1 1)
As our polygon is a square, the convex hull is the same:POLYGON ((1 1,1 5,5 5,5 1,1 1))
Set a buffer of 0 for our polygon:POLYGON ((1 1,1 5,5 5,5 1,1 1)), which is the same.
False

创建shapefile并写入数据

利用OGR还可以将刚才创建的polygon数据写入到shapefile中，需要用到ogr和osr两个库，具体过程如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
from osgeo import osr
# 1. 设置空间参考坐标系
spatialReference = osr.SpatialReference()
spatialReference.ImportFromProj4('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs')

# 2. 创建新的shapefile
driver = ogr.GetDriverByName('ESRI Shapefile')
shapeData = driver.CreateDataSource('my_polygon.shp')

# 3. 添加图层
layer = shapeData.CreateLayer('polygon_layer', spatialReference, ogr.wkbPolygon)
layerDefinition = layer.GetLayerDefn()

# 4. 添加geometry，创建feature
featureIndex = 0
feature = ogr.Feature(layerDefinition)
feature.SetGeometry(polygon)
feature.SetFID(featureIndex)

# 5. 将feature添加到图层
layer.CreateFeature(feature)
0

接下来，可以用ogrinfo查看刚才创建的shapefile：

1
!ogrinfo my_polygon.shp
INFO: Open of `my_polygon.shp'
      using driver `ESRI Shapefile' successful.
1: my_polygon (Polygon)

使用filter选择features

我们还是使用前面用过的Natural Earth数据集，使用经纬度坐标以边界框的形式创建filter，通过选择边界框以内的数据，达到选择数据子集的效果。我们将使用OGR的SpatialFilterRec方法，该方法采用四个值-minx，miny，maxx和maxy来创建边界框。在实例中，选择我们的边界框中的城市（北京周边）。 为了进一步过滤我们的结果，我们只想要中国的城市。 这意味着我们必须在for循环中使用额外的if / else语句过滤搜索结果。 具体代码如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
import os

# reference the shapefile and specify driver type
shapefile = r"data/10m_cultural/ne_10m_populated_places.shp"

driver = ogr.GetDriverByName('ESRI Shapefile')

# open the datasource with driver, zero means open in read-only mode
dataSource = driver.Open(shapefile, 0)

# use the GetLayer() function for referencing the layer that holds the data
layer = dataSource.GetLayer()

# pass in the coordinates for the dataframe to the SetSpatialFilterRect() function,
# which creates a rectangular extent and selects the features inside the extent.
layer.SetSpatialFilterRect(115, 35, 125, 45)

for feature in layer:
    # select only the cities inside of China
    # we can do this through a SQL query
    # skip the cities that are not in CHina,
    # and print the names of the cities
    if feature.GetField("ADM0NAME") != 'China':
        pass
    else:
        print(feature.GetField('NAME'))

Zhuozhou
Hengshui
Xuanhua
Beipiao
Wafangdian
...

Shapely和Fiona库

之前已经介绍过Shapeply库和Fiona库，之所以要放在一起，是因为Shapely依赖其他库实现文件读写，而Fiona刚好可以胜任。接下来，Fiona读取文件，然后由Shapely进行分析。

Shapely对象和类

Shapely库用于创建和操作2D矢量数据，而无需空间数据库，此外它也不关心投影和数据格式，而只专注于几何形状本身。 Shapely的优势在于它使用易于阅读的语法来创建可用于几何操作的各种几何形状。借助其他Python包，可以将这些几何和几何操作的结果写入矢量文件格式并在必要时进行投影。我们将介绍将pyproj和Fiona与Shapely的功能结合起来的示例。例如，我们可以使用Fiona从shapefile中读取矢量几何图形，然后使用Shapely来简化或清理现有几何图形（类似数据清洗）。清理后的几何图形可用作其他工作流程的输入，例如，用于创建专题图或执行数据科学任务。

Shapely库使用一组类，这些类是三种基本类型的几何对象（点，曲线和曲面）的实现，具体包括：

地理对象 Shapely类名
点 Point
线 LineString, Linear Ring
面 Polygon
点集合 MultiPoint
线集合 MultiLineString
面集合 MultiPolygon

Shapely的地理空间分析方法

在Shapley中，拓扑关系（例如包含、相邻等）作为集合对象的方法加以实现，Shapely还提供产生新的几何对象的方法（例如交叉、联合等）的实现。通过使用buffering方法可以使得形状更加清晰。此外，借助WKB和WKT、数组、Bumpy以及Python Geo接口，Shapely实现了与其他软件的互操作。

Fiona数据模型

尽管Fiona是OGR是Python的封装，但FIona与OGR的数据模型并不同。OGR使用数据源（data source）、图层（layer）和要素（feature），Fiona基于GeoJSON使用条目记录（term record）访问向量数据中的地理要素。一条记录具有ID、几何形状和属性键，可以用Python字典对象通过键引用记录。

创建几何形状

如同OGR一样，Shapely也可以用来创建形状，并且可以绘制出来，而不需要其他命令。

1
2
3
4
5
from shapely.geometry import Polygon

p1 = Polygon(((1, 2), (5, 3), (5, 7), (1, 9), (1, 2)))
p2 = Polygon(((6,6), (7,6), (10,4), (11,8), (6,6)))
p1

1
p2

1
2
3
4
# Point takes tuples as well as positional coordinate values
from shapely.geometry import Point
point = Point(2.0, 2.0)
point

1
2
3
4
# line geometry
from shapely.geometry import LineString
line = LineString([(0, 0), (10, 10)])
line

1
2
3
4
# linear rings
from shapely.geometry.polygon import LinearRing
ring = LinearRing([(0,0), (3,3), (3,0)])
ring

1
2
3
4
# collections of points
from shapely.geometry import MultiPoint
points = MultiPoint([(0,0), (0,3), (3,0), (3,3)])
points

1
2
3
4
# collection of lines
from shapely.geometry import MultiLineString
coords = MultiLineString([((0, 0), (1, 1)), ((-1, 0), (1, 0))])
coords

1
2
3
4
# cooection of polygons
from shapely.geometry import MultiPolygon
polygons = MultiPolygon([p1, p2])
polygons

使用Shapely的几何方法

1
dir(p1)[-20:]
['length',
 'minimum_rotated_rectangle',
 'overlaps',
 'project',
 'relate',
 'relate_pattern',
 'representative_point',
 'simplify',
 'svg',
 'symmetric_difference',
 'to_wkb',
 'to_wkt',
 'touches',
 'type',
 'union',
 'within',
 'wkb',
 'wkb_hex',
 'wkt',
 'xy']

1
2
3
print(p1.area)
print(p1.bounds)
print(p1.to_wkt())
22.0
(1.0, 2.0, 5.0, 9.0)
POLYGON ((1.0000000000000000 2.0000000000000000, 5.0000000000000000 3.0000000000000000, 5.0000000000000000 7.0000000000000000, 1.0000000000000000 9.0000000000000000, 1.0000000000000000 2.0000000000000000))

读取JSON数据

尽管Shapely不能读写数据文件，我们可以用Shapely读取json数据：

1
2
3
4
5
6
7
8
import json
from shapely.geometry import mapping, shape

p = shape(json.loads('{"type": "Polygon", "coordinates": [[[1,1], [1,3], [3,3]]]}'))
# the mapping command returns a new, independent geometry
# with coordinates copied from the context
print(json.dumps(mapping(p)))
p.area
{"type": "Polygon", "coordinates": [[[1.0, 1.0], [1.0, 3.0], [3.0, 3.0], [1.0, 1.0]]]}

2.0

Fiona库的使用

读取数据

1
2
3
4
import fiona
c = fiona.open(r"data/10m_cultural/ne_10m_admin_1_states_provinces.shp")
rec = next(iter(c))
rec.keys()
dict_keys(['type', 'id', 'properties', 'geometry'])

使用Python的ppprint(pretty-print)，打印数据集中第一个要素的各个键值如下：

1
2
3
4
5
import pprint
pprint.pprint(rec['type'])
pprint.pprint(rec['id'])
pprint.pprint(rec['properties'])
pprint.pprint(rec['geometry'].keys)
'Feature'
'0'
OrderedDict([('featurecla', 'Admin-1 scale rank'),
             ('scalerank', 3),
             ('adm1_code', 'ARG-1309'),
             ('diss_me', 1309),
             ('iso_3166_2', 'AR-E'),
             ('wikipedia', None),
             ('iso_a2', 'AR'),
             ('adm0_sr', 1),
             ('name', 'Entre RÃ\xados'),
             ('name_alt', 'Entre-Rios'),
             ('name_local', None),
             ('type', 'Provincia'),
             ('type_en', 'Province'),
             ('code_local', None),
             ('code_hasc', 'AR.ER'),
             ('note', None),
             ('hasc_maybe', None),
             ('region', None),
             ('region_cod', None),
             ('provnum_ne', 10),
             ('gadm_level', 1),
             ('check_me', 20),
             ('datarank', 3),
             ...
<built-in method keys of dict object at 0x0000025406F075E8>

1
2
3
4
5
6
7
8
# print total amount of features
print(len(c))

# print driver name
print(c.driver)

# print coordinate reference system of data
print(c.crs)
4594
ESRI Shapefile
{'init': 'epsg:4326'}

联合使用Shapely和Fiona

前面已经提到，我们可以用Fiona打开shapefile，然后用Shapely操作分析其中的要素数据。

1
2
3
src = fiona.open(r"data/10m_cultural"
                "/ne_10m_admin_1_states_provinces.shp")
print(src[1000]['properties']['name'])
Heilongjiang

从上面结果可知，第1001个feature的名字是Heilongjiang，那我们可以用shapely的shape方法进行操作：

1
2
heilongjiang = shape(src[1000]['geometry'])
heilongjiang

Shapefile作图

正如我们前面所做的，引用单独的地理要素并使用Python作图并不是那么简单。幸运的是，已经有许多人为我们造好了许多轮子，如果不将shapefile转为GeoJSON，而是直接做图的话，可以采用一下方式：

  • 使用Numpy 数组和matplotlib：将shapefile的所有坐标装入到数组当中，然后可以用各种包作图。
  • 使用Shapely，以使用地理区域名称作为键，地理要素数据作为值，从已有的shapefile创建字典，然后用shapefly作图。
  • 使用pyshp和matplotlib：类似地，pyshp库也是用来读取地理信息。
  • 使用GeoPandas和matplotlib：GeoPandas的好处在于可以将属性表存储为pandas的dataframe。

GeoPandas的使用

GeoPandas的好处是创建类似pandas的数据结构，同时还可以实现pandas所不能的地理空间功能（几何运算、数据叠加、地理编码等）。接下来，我们看一下如何使用GeoPandas实现数据处理、地理空间分析以及绘图等功能。

选择几何数据并作图

1
2
3
4
5
import geopandas as gpd
%matplotlib inline
df = gpd.read_file(r'data/10m_cultural/'
                  'ne_10m_admin_1_states_provinces.shp')
print(df[df['iso_a2']=='CN'].name)
17             Xinjiang
18               Xizang
118        Inner Mongol
124               Gansu
261              Yunnan
1000       Heilongjiang
1005              Jilin
1010           Liaoning
1270            Guangxi
1670          Guangdong
2014             Hainan
2023             Fujian
2113           Zhejiang
2114           Shanghai
2116            Jiangsu
2220           Shandong
2221              Hebei
2222            Tianjin
3075    Paracel Islands
3786            Beijing
3789            Sichuan
3790          Chongqing
3791            Guizhou
3792              Hunan
3793            Ningxia
3794            Shaanxi
3795            Qinghai
3796             Shanxi
3797            Jiangxi
3798              Henan
3799              Hubei
3800              Anhui
Name: name, dtype: object

如上所示，'iso_a2’为’CN’的共有32个，是中国大陆的省、直辖市及自治区（不含港澳台）。

1
2
beijing = df.loc[df['name']=='Shandong']
beijing.plot(figsize=(15,12))

output_70_1

1
2
china = df.loc[df['iso_a2']=='CN']
china.plot(cmap='Set2', figsize=(15,12))

output_71_1

上面绘制的是中国大陆的省级行政区（注意：西藏等地区的边界线并非我国主张的边境线）。此外，我们还可以绘制指定经纬度范围内的地图：

1
2
exp = df.cx[80:130, 15:50]
exp.plot(cmap='Set2', figsize=(15,12))

output_73_1

绘制美国山火（wildfire）地图

近两年来，美国尤其是加州山火肆虐、经久不熄，给当地造成了严重的人身和财产损失，也导致了保险行业的巨大损失，接下来我们就用Geopandas对美国的山火情况进行简单的分析。

我们所使用的数据是MTBS提供的数据，MTBS，全称Monitoring Trends in Burn Severity，是一个由多个机构联合实施的项目，旨在为1984年以来美国的火灾严重程度和波及范围提供一致性的度量。

Top20_Destruction

1
2
3
4
5
6
import matplotlib

# import the shapefile with all of the state boundaries
states = df[df['iso_a2']=='US']
p = states.plot(figsize=(15,12))
p.set_xlim((-180,-60))
(-180, -60)

output_76_1

1
2
3
# read wildfire data
fires = gpd.read_file(r'data/mtbs_fod_pts_data/mtbs_fod_pts_DD.shp')
fires.plot(markersize=1, figsize=(15,12))

output_77_1

1
2
3
print(states.crs)
print(fires.crs)
print(fires.columns)
{'init': 'epsg:4326'}
{'init': 'epsg:4269'}
Index(['Fire_ID', 'Fire_Name', 'Asmnt_Type', 'Pre_ID', 'Post_ID', 'Fire_Type',
       'ND_T', 'IG_T', 'Low_T', 'Mod_T', 'High_T', 'Ig_Date', 'Lat', 'Long',
       'Acres', 'geometry'],
      dtype='object')

通过以上代码，我们发现了两个问题：

  1. states与fires采用的是不同的投影坐标系；
  2. fires中的山火信息不包括其发生在哪个州。

针对第一个问题，我们可以通过重新投影解决：

1
fires = fires.to_crs({'init': 'epsg:4326'})

针对第二个问题，利用sjoin方法，通过空间联合，指示fires的经纬度是否位于某一个州的经纬度范围内：

1
2
3
state_fires = gpd.sjoin(fires, states[['name', 'geometry']].copy(),
                       op='within')
state_fires[['Fire_Name', 'Ig_Date', 'Lat', 'Long','name']].head()

接下来，我们就可以按照州为单位统计发生的山火次数：

1
2
counts_per_state = state_fires.groupby('name').size()
counts_per_state.sort_values(axis=0, ascending=False)
name
Florida           3710
California        1630
Kansas            1422
Idaho             1316
...
dtype: int64

可以看到，Florida、California和Kansas在1984-2016年发生的山火次数高居全美前三位。但是，数字毕竟没有图像直观，接下来，我们将各州山火情况可视化：

首先，我们可以将计算的数据作为原来shapefile的一个属性：

1
states = states.merge(counts_per_state.reset_index(name='number_of_fires'))

然后画出等值图（chorpleth map）:

1
2
3
4
5
6
7
8
import matplotlib.pyplot as plt

fig, ax = plt.subplots(1, figsize=(18,6))
ax = states.plot(column='number_of_fires',
                 cmap='OrRd', legend=True, ax=ax)
ax.set_xlim((-180, -60))
fig.suptitle('US Wildfire count per state in 1984-2015')
ax.set_axis_off()

output_89_0.png

1
2
3
4
5
6
7
import plotly.graph_objs
import plotly.offline as py_offline

scl = [[0.0, 'rgb(242,240,247)'],[0.2, 'rgb(218,218,235)'],[0.4, 'rgb(188,189,220)'],
       [0.6, 'rgb(158,154,200)'],[0.8, 'rgb(117,107,177)'],[1.0, 'rgb(84,39,143)']]
states['text'] = states['name'] + '<br>' + 'Fires:' + states['number_of_fires'].astype(str)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
data = [dict(
    type='choropleth',
    colorscale=scl,
    autocolorscale=False,
    locations=[x[-2:] for x in states['iso_3166_2']],
    z=states.number_of_fires,
    locationmode='USA-states',
    text = states.text,
    marker = dict(
        line = dict(
            color = 'rgb(255,255,255)',
            width = 2
        )
    ),
    colorbar=dict(title='Number of fires')
)]

layout = dict(
    title = 'Number of fires per state<br>(1984-2016)',
    geo = dict(
        scope='usa',
        projection=dict( type='albers usa' ),
        showlakes = True,
        lakecolor = 'rgb(255, 255, 255)',
    ),
)

fig = dict(data=data, layout=layout)

url = py_offline.iplot(fig, filename='wildfire-map.html')
py_offline.plot(fig, include_plotlyjs="cdn", output_type='div')
## 数据检查的重要性

在我们真正进入分析之前，首先最好对数据进行检查，有一个初步的认识 例如，列出有关数据集的统计信息，以显示有多少元素，以及是否存在任何缺失值。 在进行分析之前，通常需要清理数据。 以我们前面使用的山火数据为例，统计山火的发生次数与各州山火的总数，发现分州汇总之后的数据少61个。这可能意味着：

  • 有些山火的数据是缺失（无效的）；
  • 山火的经纬度坐标超出了Natrual Earth所给的坐标范围；
  • …

这就需要我们重新去检查数据的有效性，并做相应的处理。

1
2
3
4
5
6
7
print(len(fires.index))
print(counts_per_state.sum())
print(fires.empty)
print(fires.geometry.empty)
for index, row in fires.iterrows():
    if row['geometry'] =='NULL':
        print(row)
21673
21612
False
False

小结

本文介绍了三个用于处理矢量数据-OGR，Shapely和GeoPandas的Python库以及如何使用这三者进行地理空间分析和处理。 每个库都分别包含其类，方法，数据结构和常用功能。接下来，我们将探讨栅格数据的分析处理以及GDAL和Rasterio库的具体使用。