matplotlibbasemap绘制多边形区域曲线
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中⽤的相当多的包，对数组的封装堪称完美
matplotlib中subplotshapely 1.6.0:python中专门处理多边形的包，程序中⽤于判断数据点是否在指定区域内
3.代码流程
⾸先来分析下绘图的整体思路，我们所有的图形操作都是在Basemap上的操作，其余的包和数据都是辅助。步骤是先绘制指定⽅形区域的数据图，然后在将指定多边形区域外的内容去掉。
1. 导⼊相关包
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
ry import Polygon as ShapelyPolygon
ry import Point as ShapelyPoint
1. 设定些基本常数
# 绘图的经纬度范围，经度范围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
1. 读⼊离散点的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'为线性插值，当然还有另外⼀种是临近    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'，详细参数设置参看basemaptuto 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即可，参看adthedocs.io/en/latest/shapefile.html
5.绘制数据图形
# 调整格点投影坐标，这是basemap⼀个值得注意的细节，因为投影格式的区别，需要将经纬度的数值［xi,yi］投影到实际的制图坐标
x1, y1 = map(xi, yi)
# ⽹格化经纬度⽹格：调⽤numpy的meshgrid函数来⽣成后⾯需要的⽹格化的数据格点xx,yy
xx, yy = np.meshgrid(x1, y1)
# 绘图等值线：这⾥使⽤basemap的两种绘制数据的⽅法，pcolor和等值线contour，详细参看adthedocs.io/en/latest/plotting_data.html # 等值线的相关参数api地址：/api/pyplot_api.html#ur
PCM = map.pcolor(xx, yy, zi, alpha=0.8)
CS = ur(xx, yy, zi,\
alpha=0.8,
linestyles = 'dashed',
levels = np.arange(np.min(Temperature),np.max(Temperature),1)
)
# matplotlib中通过调⽤plt.clabel给等值线上标上相关的⽂本数据，相关的参数设置参看：/devdocs/api/_as_gen/matplotlib.axes.Axes.clab CS_label = plt.clabel(CS, inline=True, inline_space=10, fontsize=8, fmt='%2.0f', colors='k')
6.裁剪边缘
# 调⽤shapefile模块读⼊shape⽂件，shapefile的官⽅模块路径为：/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, ansData)
＃等值线的导出类型CS中包含collections，可以直接通过其中的set_clip_path函数来切除边缘，不过只限于contour等值线层，详细参见basemaptutorial.r for contour llections:
contour.set_clip_path(clip)
# 下⾯通过调⽤shapely的多边形来将之前⽤clabel标上的⽂本数字进⾏过滤删除边界外的数据，参见stackoverflow/questions/42929534/set-clip-pa clip_map_shapely = ShapelyPolygon(vertices)
for text_object in CS_label:
if not clip_ains(ShapelyPoint(_position())):
text_object.set_visible(False)
# 因为pcolor返回的是⼀个polygon的对象，因此可以直接调⽤set_clip_path来切除该图层的边缘
PCM.set_clip_path(clip)
plt.show()
5.总结
图虽看着简单，其中的细节处理确实值得深⼊研究，笔者也未完全熟练掌握，其中⽤了pcolor和contour两个函数来描绘数据图层，是因为基于笔者的数据⼀个c ontourf绘出结果奇怪。当然，如果只
有⼀个contourf图层则可以⽤官⽅推荐的边界删除⽅法直接去除。adthedocs.io/en/latest/clip.html

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。