在平常的python使用中，有些时候需要基于gis的地理数据绘制相关的数据图形，如上图所示，python中的matplotlib绘图包当然能够胜任这个工作，但实际操作中国确实有很多细节需要注意。
首先来分析下绘图的整体思路，我们所有的图形操作都是在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()