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上的操作，其餘的包和資料都是輔助。步驟是先繪製指定方形區域的資料圖，然後在將指定多邊形區域外的內容去掉。

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
from shapely.geometry import Polygon as ShapelyPolygon
from shapely.geometry 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'為線性插值，當然還有另外一種是臨近插值'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