基于Python的地理空间分析（六）：GeoDB数据可视化分析

引言

前面已经介绍了如何安装PostGIS、创建数据表、添加数据并执行基本的空间查询。我们将在此基础上进一步学习如何使用地理空间数据库进行分析和地图可视化。我们将采用美国新墨西哥州Albuquerque的犯罪数据作为例子完成一系列的分析和图表报告。

具体而言，在这篇Notebook里，我们打算创建一个关于犯罪的仪表盘（crime dashboard），大概可以分为三步：

1. 创建数据库
2. 创建数据分析工具
3. 可视化分析结果

创建犯罪案件数据库

如前所述，我们采用的是美国新墨西哥州Albuquerque市的开放数据（http://www.cabq.gov/abq-data/ ），该市提供了犯罪事件、区域指挥（Area commands）和巡逻区（beats）等相关数据集，我们将通过联合使用这三个数据集进行初步的分析。

注：关于数据集中Layer不同字段的说明见 http://coagisweb.cabq.gov/arcgis/rest/services/public/APD_Incidents/MapServer/0 。

创建数据表

分别为犯罪事件、区域指挥和巡逻区创建相应的数据表，导入所需要的库：

import psycopg2 # connect to PostGIS
import requests # grab data
import datetime # convert geogson
from shapely.geometry import Point, Polygon, MultiPolygon, mapping

# create a connection
connection = psycopg2.connect(database='pythonspatial',
                             user='postgres', password='postgres')

cursor = connection.cursor()

proxies = {
  'http': 'socks5://127.0.0.1:13579',
  'https': 'socks5://127.0.0.1:13579',
}

# create the areacommand table with name and geometry field
# the area command filed has a lenth of 20 in the ArcServer service
cursor.execute('create table areacommand (id SERIAL PRIMARY KEY, name VARCHAR(20),'
               'geom GEOMETRY)')

# create the beats table
cursor.execute('CREATE TABLE beats (id SERIAL PRIMARY KEY, beat VARCHAR(6),'
               'agency VARCHAR(3), areacomm VARCHAR(15), geom GEOMETRY)')

# create the incidents table
cursor.execute('CREATE TABLE incidents (id SERIAL PRIMARY KEY, address VARCHAR(72), '
               'crimetype VARCHAR(255), date DATE, geom GEOMETRY)')

connection.commit()

填充数据

# grab the area commands and insert them into table
url='http://coagisweb.cabq.gov/arcgis/rest/services/public/adminboundaries/MapServer/8/query'
params={"where":"1=1","outFields":"*","outSR":"4326","f":"json"}

# query the url passing parameters
r=requests.get(url,params=params, proxies=proxies)
data=r.json()

# insert the data
for acmd in data['features']:
    polys=[]

    # loop through the rings
    # A polygon contains an array of rings and a spatialReference.
    # Each ring is represented as an array of points.
    for ring in acmd['geometry']['rings']:
        polys.append(Polygon(ring))
    p=MultiPolygon(polys)
    name=acmd['attributes']['Area_Command']

    # uses Shapely to convert the coordinates to WKT (p.wkt).
    cursor.execute("INSERT INTO areacommand (name, geom) VALUES ('{}',ST_GeomFromText('{}'))".format(name, p.wkt))

connection.commit()

上面SQL查询的关键是INSERT INTO table (field, field) VALUES (value,value)命令，其中潜逃了ST_GeometryFromText()命令。ST_GeometryFromText和ST_GeomFromText一样，根据WKT描述的对象返回一个具体的geometry对象，WKT可以是POLYGON、MULTIPOLYGON、LINESTRING等。

此外，关于ESRI ArcServer查询的参数命令，可以查询： http://coagisweb.cabq.gov/arcgis/sdk/rest/index.html#/Query_Map_Service_Layer/02ss0000000r000000/

采用同样的套路，向beats表写入数据：

url='http://coagisweb.cabq.gov/arcgis/rest/services/public/adminboundaries/MapServer/9/query'
params={"where":"1=1","outFields":"*","outSR":"4326","f":"json"}
r=requests.get(url,params=params, proxies=proxies)
data=r.json()

for acmd in data['features']:
    polys=[]
    for ring in acmd['geometry']['rings']:
        polys.append(Polygon(ring))

    p=MultiPolygon(polys)

    beat = acmd['attributes']['BEAT']
    agency = acmd['attributes']['AGENCY']
    areacomm = acmd['attributes']['AREA_COMMA']

    cursor.execute("INSERT INTO beats (beat, agency,areacomm,geom) VALUES('{}','{}','{}', ST_GeomFromText('{}'))".format(beat,agency,areacomm,p.wkt))

connection.commit()

最后向incidents表写入数据：

url='http://coagisweb.cabq.gov/arcgis/rest/services/public/APD_Incidents/MapServer/0/query'
params={"where":"1=1","outFields":"*","outSR":"4326","f":"json"}
r=requests.get(url,params=params, proxies=proxies)
data=r.json()
for a in data["features"]:
    address=a["attributes"]["BlockAddress"]
    crimetype=a["attributes"]["IncidentType"]
    if a['attributes']['date'] is None:
        pass
    else:
        date = datetime.datetime.fromtimestamp(a['attributes']['date'] /1e3).date()

    try:
        p=Point(float(a["geometry"]["x"]),float(a["geometry"]["y"]))
        cursor.execute(("INSERT INTO incidents (address,crimetype,date, geom) VALUES('{}','{}','{}',"
                        "ST_GeomFromText('{}'))").format(address,crimetype,str(date), p.wkt))
    except KeyError:
        pass

connection.commit()

与areacommand和beats表相比，在incidents表写入数据过程中，我们做了一定程度的数据清洗，因为可能有些记录缺少date（直接pass）或者坐标系信息（用try catch语句去捕获异常）。

现在我们已经完成了数据的导入工作，接下来我们就可以在PostgreSQL数据库中进行数据的查询并进一步分析和呈现数据。

数据查询结果可视化

之前，我们已经从地理空间数据库中查询数据，并返回WKT格式的数据，但是仅仅是一系列坐标点却并不是想象中的地图。接下来我们使用ipyleaflet和Jupyter实现查询结果的地图可视化。

ipyleaflet实现了在Jupyter中交互式操作地图，具体的代码、用例，可以访问：https://github.com/ellisonbg/ipyleaflet

准备工作

# Install ipyleaflet
!pip install ipyleaflet

# enable the extension for notbook
# can be skipped for notebook 5.3 and above
!jupyter nbextension enable --py --sys-prefix ipyleaflet

# enable widgetsnbextension for notebook
!jupyter nbextension enable --py --sys-prefix widgetsnbextension

# for jupyter lab
!jupyter labextension install jupyter-leaflet

!jupyter labextension install @jupyter-widgets/jupyterlab-manager

注意：以上命令实际应当在命令行中运行（去掉感叹号！）

重启过Jupyter之后就可以使用了。

ipyleaflet交互式操作地图

import psycopg2
from shapely.geometry import Point, Polygon, MultiPolygon
from shapely.wkb import loads
from shapely.wkt import dumps, loads
import datetime
import json
import ipyleaflet as ipl

connection = psycopg2.connect(database='pythonspatial', user='postgres', password='postgres')
cursor = connection.cursor()
cursor.execute('SELECT name, ST_AsGeoJSON(geom) from areacommand')
c = cursor.fetchall()
c[0][0]

'FOOTHILLS'

之前我们用ST_AsText()命令查询PostgreSQL数据库返回geometry数据，由于GeoJSON更加便于地图呈现，所以我们采用ST_AsGeoJSON()命令。以上代码返回了areacommand表中的所有记录，每条记录包括name和geometry两个字段。

ST_AsText和ST_AsGeoJSON只是PostGIS获得geometry的17种方法里的两种，更多方法参见：https://postgis.net/docs/reference.html#Geometry_Accessors

接下来就该创建地图并显示获得的GeoJSON数据了：

center = [35.106196, -106.629515]
# the zoom level—the higher the number, the closer the zoom.
zoom = 9
map = ipl.Map(center=center, zoom=zoom)

for x in c:
    layer = json.loads(x[1])
    layer_geojson = ipl.GeoJSON(data=layer)
    map.add_layer(layer_geojson)
map

ipyleaflet控件在转换为markdown文件时无法显示，所以我用folium来作图。

import folium
map = folium.Map(center, zoom_start=9,control_scale=True)
for x in c:
    layer = json.loads(x[1])
    folium.GeoJson(layer).add_to(map)
map

# map the beats in a similar way
cursor.execute('SELECT beat, ST_AsGeoJSON(geom) from beats')
c = cursor.fetchall()
for x in c:
    layer = json.loads(x[1])
    layer_geojson = ipl.GeoJSON(data=layer)
    map.add_layer(layer_geojson)

同样地，我们也可以在地图显示incidents的数据，但是由于incidents的数目太多，全显示出来反而让人无所适从。因此，我们将通过条件的限定，优化我们的显示结果。

根据日期显示incidents

center = [35.106196, -106.629515]
# the zoom level—the higher the number, the closer the zoom.
zoom = 9
map = folium.Map(center, zoom_start=zoom)
date = datetime.datetime.strptime('20181031', '%Y%m%d').date()
print(date)
cursor.execute("SELECT address, crimetype, date, ST_AsGeoJSON(geom) from incidents where date = '{}'".format(str(date)))
incidents_date = cursor.fetchall()

for x in incidents_date:
    layer = json.loads(x[3])
    folium.Marker([layer['coordinates'][1], layer['coordinates'][0]]).add_to(map)
map.save('crimes_by_date.html')
map

2018-10-31

此外，我们也可以查询某天之后或者某天之前一段时间内的incidents，例如查询最近10天以内的犯罪事件：

cursor.execute("select * from incidents where date >= Now() - interval '10 day'")

根据位置查询incidents

center = [35.106196, -106.629515]
# the zoom level—the higher the number, the closer the zoom.
zoom = 10
map = folium.Map(center, zoom_start=zoom)

cursor.execute("SELECT ST_AsGeoJSON(i.geom) FROM incidents i JOIN areacommand acmd ON ST_Intersects(acmd.geom, i.geom)"
               " WHERE acmd.name like'FOOTHILLS' and date >= NOW() - interval '12 day';")
crime = cursor.fetchall()
for x in crime:
    layer = json.loads(x[0])
    folium.Marker([layer['coordinates'][1], layer['coordinates'][0]]).add_to(map)


cursor.execute("SELECT name, ST_AsGeoJSON(geom) from areacommand WHERE name like 'FOOTHILLS'")
foothills= cursor.fetchall()

layer = json.loads(foothills[0][1])
folium.GeoJson(layer).add_to(map)
map.save('crime_by_loc.html')
map

{'type': 'Point', 'coordinates': [-106.509618200809, 35.0694532721882]}

上面查询了在Foothills区域最近10天的事件。在具体的操作上，需要使用JOIN命令实现多个表的合并查询，考虑到空间合并查询的特殊性，PostGIS提供了ST_Intersects, ST_Contains, and ST_DWithin等多种合并方式，其中我们用的是ST_Intersects。关于命令的使用方法参见：

• http://postgis.net/workshops/postgis-intro/joins.html
• https://postgis.net/docs/ST_Intersects.html

类似地，我们也可以查询不同巡逻区的犯罪事件：

center = [35.106196, -106.629515]
# the zoom level—the higher the number, the closer the zoom.
zoom = 11
map = folium.Map(center, zoom_start=zoom)

cursor.execute("SELECT ST_AsGeoJSON(i.geom) FROM incidents i JOIN beats b ON "
               "ST_Intersects(b.geom, i.geom) WHERE b.beat in ('336','523','117','226','638','636') "
               "and date >= NOW() - interval '20 day';")
crime = cursor.fetchall()
for x in crime:
    layer = json.loads(x[0])
    folium.Marker([layer['coordinates'][1], layer['coordinates'][0]]).add_to(map)

cursor.execute("SELECT beat, ST_AsGeoJSON(geom) from beats WHERE beat in ('336','523','117','226','638','636')")
beats = cursor.fetchall()
for b in beats:
    layer = json.loads(b[1])
    folium.GeoJson(layer).add_to(map)
map.save('crimes_by_beats.html')
map

前面我们展示的结果都是直接从数据库中读取的数据，接下来我们实现地理空间分析结果的可视化。

Buffer的使用

首先用Shapely创建一个点对象，然后用mapping()转换为GeoJSON。将这个GeoJSON作为输入发送给PostGIS，使用PostGIS的空间分析功能返回所需要的结果，而这些并不要求是数据库表中的数据。

center = [35.106196, -106.629515]
# the zoom level—the higher the number, the closer the zoom.
zoom = 12
map = folium.Map(center, zoom_start=zoom)

from shapely.geometry import mapping
p = Point([-106.578677,35.062485])
# pgeojson = mapping(p)
# map.add_layer(player)
folium.Marker([35.062485, -106.578677]).add_to(map)

cursor.execute("SELECT ST_AsGeoJSON(ST_Buffer(ST_GeomFromText('{}')::geography,1500));".format(p.wkt))
buff = cursor.fetchall()
buffer = json.loads(buff[0][0])
folium.GeoJson(buffer).add_to(map)
map.save('point_buffer.html')
map

那么，为了查询在buffer内的犯罪事件：

• 首先获得buffer区域的wkt（用ST_GeomFromTexxt()）；
• 然后使用loads()命令创建一个shapely的多边形对象；
• ST_Intersects()查找多边形内的incidents

center = [35.062485, -106.578677]
# the zoom level—the higher the number, the closer the zoom.
zoom = 12
map = folium.Map(center, zoom_start=zoom)
p = Point([-106.578677,35.062485])
folium.Marker([35.062485, -106.578677]).add_to(map)

cursor.execute("SELECT ST_AsText(ST_Buffer(ST_GeomFromText('{}')"
               "::geography, 1500));".format(p.wkt))
bufferwkt = cursor.fetchall()
b = loads(bufferwkt[0][0])

buffer_json = mapping(b)
folium.GeoJson(buffer_json).add_to(map)


cursor.execute("SELECT ST_AsGeoJSON(incidents.geom) FROM incidents where "
               "ST_Intersects(ST_GeomFromText('{}'), incidents.geom) and date >= "
               "Now() - interval '12 day';".format(b.wkt))
crime = cursor.fetchall()
for x in crime:
    layer = json.loads(x[0])
    folium.Marker([layer['coordinates'][1], layer['coordinates'][0]]).add_to(map)
map.save('crimes_in_buffer.html')
map

k近邻分析

对于某一个点，如果想查询其周围指定范围内的对象，用buffer非常方便；但是如果想查询最近的n个对象，则应该使用k近邻分析(k nearest neighbors)。例如查找p周围最近的10个事件，可以分为以下步骤：

• 用Shapely创建查询点的对象；
• 用ST_Distnace()计算距离
• 利用Order by、<->和limit查找k近邻。

center = [35.062485, -106.578677]
# the zoom level—the higher the number, the closer the zoom.
zoom = 12
map = folium.Map(center, zoom_start=zoom)

p = Point([-106.578677,35.062485])
cursor.execute("SELECT ST_AsGeoJSON(incidents.geom),ST_Distance("
               "incidents.geom::geography,ST_GeometryFromText('{}')::geography)"
               "from incidents ORDER BY incidents.geom<->ST_GeometryFromText('{}') "
               "LIMIT 5".format(p.wkt,p.wkt))

c = cursor.fetchall()
for x in c:
    layer=json.loads(x[0])
    folium.Marker([layer['coordinates'][1], layer['coordinates'][0]]).add_to(map)
    print(layer)
map.save('crime_knn.html')
map

{'type': 'Point', 'coordinates': [-106.577602528136, 35.0619088546598]}
{'type': 'Point', 'coordinates': [-106.578702884409, 35.0637239657475]}
{'type': 'Point', 'coordinates': [-106.578702884409, 35.0637239657475]}
{'type': 'Point', 'coordinates': [-106.578702884409, 35.0637239657475]}
{'type': 'Point', 'coordinates': [-106.578702884409, 35.0637239657475]}

交互式控件

在Jupyter当中可以通过引入控件提供交互分析。例如，引入DatePicker，我们使用DataPicker作为装饰器，当DataPicker发生变化时，theDate会被调用，打印DataPicker的日期。借助这个想法，我们可以将日期作为SQL查询以及地图可视化的参数，即当DataPicker发生变化时，进行新的查询，并刷新显示结果：

from ipywidgets import interact, interactive, fixed, interact_manual, DatePicker

@widgets.interact(x=DatePicker())
def map_by_date(x):
    if x:
        for l in map.layers[1:]:
            map.remove_layer(l)
        cursor.execute("SELECT ST_AsGeoJSON(geom) from incidents where "
                       "date = '{}'".format(str(x)))
        c = cursor.fetchall()
        for x in c:
            layer = json.loads(x[0])
            layer_geojson = ipl.GeoJSON(data=layer)
            map.add_layer(layer_geojson)
        return len(c)
    else:
        pass

center = [35.106196, -106.629515]
# the zoom level—the higher the number, the closer the zoom.
zoom = 10
map = ipl.Map(center=center, zoom=zoom)
map

@widgets.interact(x='None')
def map_by_area(x):
    if x:
        for l in map.layers[1:]:
            map.remove_layer(l)
        cursor.execute("SELECT ST_AsGeoJSON(i.geom) FROM incidents i JOIN "
                       "areacommand acmd ON ST_Intersects(acmd.geom, i.geom)"
                       " WHERE acmd.name like '{}' and date >= NOW() - interval "
                       "'10 day';".format(x))
        c = cursor.fetchall()

        for x in c:
            layer = json.loads(x[0])
            layer_geojson = ipl.GeoJSON(data=layer)
            map.add_layer(layer_geojson)
        return c
    else:
        pass
center = [35.106196, -106.629515]
# the zoom level—the higher the number, the closer the zoom.
zoom = 10
map = ipl.Map(center=center, zoom=zoom)
map

基于图表的可视化

除了直接在地图上呈现数据以外，还可以用一些常用的图表从不同维度展示和分析数据。为了方便作图，我们首先从PostgreSQL读取数据存储为pandas.DataFrame格式。

import pandas as pd
d = datetime.datetime.strptime('20180101', '%Y%m%d').date()
cursor.execute("SELECT date, count(date) from incidents where date >'{}'"
               "group by date".format(str(d)))
df = pd.DataFrame(cursor.fetchall(), columns=['date', 'counts'])

df.date = pd.to_datetime(df.date)
print(df.shape)
df.head()

(91, 2)

	date	counts
0	2018-09-25	185
1	2018-10-07	696
2	2019-02-16	12
3	2018-12-23	24
4	2018-09-03	151

import plotly.graph_objs as go
import plotly.offline as py_offline

trace = go.Bar(x=df.date, y=df.counts)
margin = go.layout.Margin(
    l=40,
    r=40,
    t=60,
    b=80,
    pad=20)
layout = go.Layout(margin=margin, template='ggplot2')
data = [trace]
fig = go.Figure(data=data, layout=layout)
py_offline.iplot(fig)
py_offline.plot(fig, include_plotlyjs="cdn", output_type='div', filename='crime_bar.html')

从图上我们可以看可以发现些有意思的现象：

1. 犯罪事件的发生存在一定的周期性（似乎以周为单位）；
2. 从2018年11月中旬以来，犯罪事件明显下降。

我们再来看一下空间上的情况：

cursor.execute("SELECT beats.beat, beats.agency, count(incidents.geom) "
               "as crimes from beats left join incidents on ST_Contains("
               "beats.geom, incidents.geom) group by beats.beat, beats.agency")
area = pd.DataFrame(cursor.fetchall(), columns=['area', 'agency',
                                               'crimes'])
print(area.shape)
print(area.head())
trace = go.Bar(x=area.area, y=area.crimes)
margin = go.layout.Margin(
    l=40,
    r=40,
    t=60,
    b=80,
    pad=20)
layout = go.Layout(margin=margin, template='ggplot2')
data = [trace]
fig = go.Figure(data=data, layout=layout)
py_offline.iplot(fig)
py_offline.plot(fig, include_plotlyjs="cdn", output_type='div')

(141, 3)
  area agency  crimes
0         APD       0
1         AVI       0
2  100    APD       0
3  111    APD     102
4  112    APD     235

在上面的代码中，根据beat的区域对incidents进行归类，其中用到`left join`，从而可以获得incidents数目为0的beats。从结果来看，似乎犯罪情况在区域分布上具有明显的集聚性，而且存在某些高犯罪率的区域。

数据库触发器（Triggers）

在数据库中，触发器是特定事件出现的时候，自动执行的代码块。触发器可以视作一种特殊的存储过程，但是它的执行不是由程序调用，也不需要手动操作，而是由事件来触发，即由对表进行增删改操作所触发的。当对一个数据库或表进行增删改（ Insert，Delete，Update）的时就会激活触发器。

空间数据库也支持触发器，我们可以在Python中创建触发器。例如，我们可以创建一个触发器，防止不完整记录添加到数据库当中(使用PL/pgSQL语言），首先创建一个函数：

query=('CREATE FUNCTION newcrime()\n'
       'RETURNS trigger\n'
       'AS $newcrime$\n'
       'BEGIN\n'
       'IF NEW.crimetype IS NULL THEN\n'
       'RAISE EXCEPTION' +" '% Must Include Crime Type', NEW.address;"+'\n'
       'END IF;\n'
       'RETURN NEW;\n'
       'END;\n'
       '$newcrime$\n'
       'LANGUAGE \'plpgsql\';'
)
print(query)
cursor.execute(query)

CREATE FUNCTION newcrime()
RETURNS trigger
AS $newcrime$
BEGIN
IF NEW.crimetype IS NULL THEN
RAISE EXCEPTION '% Must Include Crime Type', NEW.address;
END IF;
RETURN NEW;
END;
$newcrime$
LANGUAGE 'plpgsql';



---------------------------------------------------------------------------

ProgrammingError                          Traceback (most recent call last)

<ipython-input-271-36e72f7bf30a> in <module>
     12 )
     13 print(query)
---> 14 cursor.execute(query)


ProgrammingError: 错误:  带相同参数类型的函数 "newcrime" 已经存在

然后创建一个触发器调用该函数：

query=('CREATE TRIGGER newcrime BEFORE INSERT OR UPDATE ON incidents FOR EACH '
       'ROW EXECUTE PROCEDURE newcrime()')
cursor.execute(query)
connection.commit()

我们插入一条记录，作为测试：

p=Point([-106,35])
address="123 Sesame St"
cursor.execute("INSERT INTO incidents (address, geom) VALUES ('{}',"
               "ST_GeomFromText('{}'))".format(address, p.wkt))

---------------------------------------------------------------------------

InternalError                             Traceback (most recent call last)

<ipython-input-274-67d485332cb4> in <module>
      2 address="123 Sesame St"
      3 cursor.execute("INSERT INTO incidents (address, geom) VALUES ('{}',"
----> 4                "ST_GeomFromText('{}'))".format(address, p.wkt))
      5 cursor.execute("SELECT * FROM incidents WHERE address LIKE '123 Sesame St'")
      6 cursor.fetchall()


InternalError: 错误:  123 Sesame St Must Include Crime Type
CONTEXT:  在RAISE的第4行的PL/pgSQL函数newcrime()

cursor.execute("SELECT * FROM incidents WHERE address LIKE '123 Sesame St'")
cursor.fetchall()

除了前面用到的PL/pgSQL语言以外，还可以使用Stack Builder安装其他语言的触发器，即选择安装’Categories’->‘Add-ons, tools and utilities’ -> ‘EDB Language Pack’。

小结

本文主要介绍了如何使用地理数据库和地图数据可视化进行查询、分析和结果的呈现，接下来我们将介绍如何使用QGIS进行地理空间分析。