matplotlib basemap 绘制多边形区域曲线

1.简介

在平常的python使用中，有些时候需要基于gis的地理数据绘制相关的数据图形，如上图所示，python中的matplotlib绘图包当然能够胜任这个工作，但实际操作中国确实有很多细节需要注意。

2.运行环境

python 3.5 :笔者使用python3.x，当然2.7也是一样可以的

matplotlib 1.5.3

Basemap 1.1.0:basemap的安装包比较大，快捷安装可以参考笔者的另一篇博客文章

pyshp 1.2.12:python中处理shapefile的包，程序中需要读入基于shapefile的gis地理数据作为地形边界

json 2.5.1(可选)：json包用于读取数据文件，读者可以根据自己的数据类型选取相应的数据处理包

numpy 1.11.1:python中用的相当多的包，对数组的封装堪称完美

shapely 1.6.0:python中专门处理多边形的包，程序中用于判断数据点是否在指定区域内

3.代码流程

首先来分析下绘图的整体思路，我们所有的图形操作都是在Basemap上的操作，其余的包和数据都是辅助。步骤是先绘制指定方形区域的数据图，然后在将指定多边形区域外的内容去掉。

导入相关包

import numpy as np
from mpl_toolkits.basemap import Basemap
from matplotlib.path import Path
from matplotlib.patches import PathPatch
import matplotlib.pyplot as plt
import shapefile as sf
from matplotlib.mlab import griddata
import json
from shapely.geometry import Polygon as ShapelyPolygon
from shapely.geometry import Point as ShapelyPoint

设定些基本常数

# 绘图的经纬度范围，经度范围113.4-115.8,纬度范围37.3-38.9
PLOTLATLON = [113.4, 115.8, 37.3, 38.9]
# 地图的中心点经纬度坐标
SJZ_Center = [38.03, 114.48]
# 插值时设置的格点数，需要依据实际情况动态调整
PLOT_Interval = 100

读入离散点的json数据,插值到指定的格点上

with open('path/to/your/data/data.json') as data_file:
data = json.load(data_file)
# 经度数据list
x = []
# 纬度数据list
y = []
# 实际数据list,这里以温度Temperature为例
Temperature = []
#数据预处理，剔除掉错误或缺测值，只读取有用的数据点
for station in data:
if not station['TEM'] or station['TEM'] == '999999':
# Temperature.append(np.NaN)
pass
else:
Temperature.append(float(station['TEM']))
x.append(float(station['Lon']))
y.append(float(station['Lat']))

# 定义经度和纬度两个方向上的插值格点，这里用到之前全局定义的插值数PLOT_Interval
xi = np.linspace(PLOTLATLON[0], PLOTLATLON[1], PLOT_Interval)
yi = np.linspace(PLOTLATLON[2], PLOTLATLON[3], PLOT_Interval)
# 使用matplotlib.mlab下的griddata函数来将离散的数据格点插值到固定网格上［xi,yi］，这里插值函数设置值为'linear'为线性插值，当然还有另外一种是临近插值'nn'，详细请参看https://matplotlib.org/api/mlab_api.html
zi = griddata(np.array(x), np.array(y), np.array(Temperature), xi, yi, interp='linear')

4.创建basemap地图

# matplotlib的基本制图初始化操作
fig = plt.figure()
ax = fig.add_subplot(111)
# 创建basemap，参数中设定图像左下角和右上角的经纬度，共4个值，basemap包含的地图投影有很多种，这里选取'lcc',详细参数设置参看http://basemaptutorial.readthedocs.io/en/latest/basemap.html
map = Basemap(llcrnrlon=PLOTLATLON[0],llcrnrlat=PLOTLATLON[2],\
urcrnrlon=PLOTLATLON[1],urcrnrlat=PLOTLATLON[3],\
resolution='i', projection='lcc', lat_0 = SJZ_Center[0], lon_0 = SJZ_Center[1])
# basemap当然也可以读取shape文件，调用下面这个readshapefile即可,参看http://basemaptutorial.readthedocs.io/en/latest/shapefile.html
map.readshapefile('path/to/your/shapefile/path', 'shapefile_alias')

5.绘制数据图形

# 调整格点投影坐标,这是basemap一个值得注意的细节，因为投影格式的区别，需要将经纬度的数值［xi,yi］投影到实际的制图坐标
x1, y1 = map(xi, yi)
# 网格化经纬度网格：调用numpy的meshgrid函数来生成后面需要的网格化的数据格点xx,yy
xx, yy = np.meshgrid(x1, y1)
# 绘图等值线:这里使用basemap的两种绘制数据的方法，pcolor和等值线contour，详细参看http://basemaptutorial.readthedocs.io/en/latest/plotting_data.html
# 等值线的相关参数api地址：https://matplotlib.org/api/pyplot_api.html#matplotlib.pyplot.contour
PCM = map.pcolor(xx, yy, zi, alpha=0.8)
CS = map.contour(xx, yy, zi,\
alpha=0.8,
linestyles = 'dashed',
levels = np.arange(np.min(Temperature),np.max(Temperature),1)
)
# matplotlib中通过调用plt.clabel给等值线上标上相关的文本数据，相关的参数设置参看：https://matplotlib.org/devdocs/api/_as_gen/matplotlib.axes.Axes.clabel.html
CS_label = plt.clabel(CS, inline=True, inline_space=10, fontsize=8, fmt='%2.0f', colors='k')

6.裁剪边缘

# 调用shapefile模块读入shape文件，shapefile的官方模块路径为：https://pypi.python.org/pypi/pyshp
sjz = sf.Reader('path/to/your/shapefile/path')
for shape_rec in sjz.shapeRecords():
vertices = []
codes = []
pts = shape_rec.shape.points
prt = list(shape_rec.shape.parts) + [len(pts)]
for i in range(len(prt) - 1):
for j in range(prt[i], prt[i+1]):
vertices.append(map(pts[j][0], pts[j][1]))
codes += [Path.MOVETO]
codes += [Path.LINETO] * (prt[i+1] - prt[i] -2)
codes += [Path.CLOSEPOLY]
clip = Path(vertices, codes)
clip = PathPatch(clip, transform=ax.transData)

＃等值线的导出类型CS中包含collections，可以直接通过其中的set_clip_path函数来切除边缘，不过只限于contour等值线层，详细参见http://basemaptutorial.readthedocs.io/en/latest/clip.html
for contour in CS.collections:
contour.set_clip_path(clip)

# 下面通过调用shapely的多边形来将之前用clabel标上的文本数字进行过滤删除边界外的数据，参见https://stackoverflow.com/questions/42929534/set-clip-path-on-a-matplotlib-clabel-not-clipping-properly
clip_map_shapely = ShapelyPolygon(vertices)
for text_object in CS_label:
if not clip_map_shapely.contains(ShapelyPoint(text_object.get_position())):
text_object.set_visible(False)
# 因为pcolor返回的是一个polygon的对象，因此可以直接调用set_clip_path来切除该图层的边缘
PCM.set_clip_path(clip)

plt.show()

5.总结

图虽看着简单，其中的细节处理确实值得深入研究，笔者也未完全熟练掌握，其中用了pcolor和contour两个函数来描绘数据图层，是因为基于笔者的数据一个contourf绘出结果奇怪。当然，如果只有一个contourf图层则可以用官方推荐的边界删除方法直接去除。http://basemaptutorial.readthedocs.io/en/latest/clip.html