python 畫廣東省等壓線圖

最近開發時要實現一個業務邏輯：

• 調用中國氣象數據網API接口獲取廣東省實時氣象數據
• 根據數據，基於廣東省地圖渲染等壓線圖

最終效果圖是這樣的：

• 首先是獲取實時氣壓數據，由於中國氣象數據網每次只能獲得30個站點的氣象數據，而廣東省共有86個氣象站點，所以分成3批獲取，存入數組。
  獲取到的數據格式是[{Station_Id_C,Year,Mon,Day,Hour,PES}]，然后遍歷氣象站點列表，把對應的經緯度(longitude、latitude)存入數組中，得到了這樣的三元組數組(longitude、latitude、value)，用pandas處理一下：

x = df['longitude'].values
y = df['latitude'].values
z = df[request.GET.get("c")].values

• 然后通過numpy庫中的 meshgrid 函數把x、y兩個一維數組矢量化為兩個二維數組X , Y = np.meshgrid(x,y)
• 用zip函數取出所有點points = [[a, b] for a, b in zip(x, y)]
• 使用scipy庫中的 griddata 函數插值數據集，便於形成N條等壓線Z = griddata(points, z, (X, Y))
• 為了與前端采用的geojson作為地圖數據統一，選用Geopands庫作為地圖渲染庫，首先打開廣東省的geojson文件df = gpd.read_file(os.path.dirname(__file__) + '\\gd.json')，然后將地圖渲染到matplotlib的畫布上ax = df.plot(figsize=(10, 10), alpha=0.2, edgecolor='k')
• 利用 contour 函數將等壓線渲染到畫布上C = ax.contour(X, Y, Z)，再標注等壓線的數值plt.clabel(C, inline=True, fontsize=10)
• 去除x軸、y軸plt.xticks(()) plt.yticks(())，調用show函數查看圖像
• 實際業務中的代碼如下，為了用戶下載圖片浪費服務器資源，選擇將靜態圖片保存至服務器端：

def get_contour(request):
    index_url = "http://api.data.cma.cn:8090/api?"
    get_params = {
        "dataFormat": "json",
        "interfaceId": "getSurfEleByTimeRangeAndStaID",
        "dataCode": "SURF_CHN_MUL_HOR",
        "timeRange": "[" + datetime.strptime(request.GET.get('f'), "%Y-%m-%dT%H:%M").strftime(
            "%Y%m%d%H%M%S") + "," + datetime.strptime(request.GET.get('g'), "%Y-%m-%dT%H:%M").strftime(
            "%Y%m%d%H%M%S") + "]",
        "staIDs": int(request.GET.get('e')),
        "elements": 'Station_Id_C,Year,Mon,Day,Hour,' + request.GET.get("c")
    }

    session = requests.Session()
    f = session.get(index_url + parse.urlencode(get_params))
    s = json.loads(f.text)
    l = s['DS'][len(s['DS']) - 1]
    date = datetime(year=int(l['Year']), month=int(l['Mon']),
                    day=int(l['Day']), hour=int(l['Hour']))
    s = list(DCmaStation.objects.all().values())
    d = []
    for i in range(int(len(s) / 30) + 1):
        n = [j['id'] for j in s[i * 30: (i + 1) * 30]]
        get_params = {
            "dataFormat": "json",
            "interfaceId": "getSurfEleByTimeRangeAndStaID",
            "dataCode": "SURF_CHN_MUL_HOR",
            "timeRange": "[" + date.strftime("%Y%m%d%H%M%S") + "," + date.strftime("%Y%m%d%H%M%S") + "]",
            "staIDs": str(n),
            "elements": 'Station_Id_C,Year,Mon,Day,Hour,' + request.GET.get("c")
        }
        f = session.get(index_url + parse.urlencode(get_params))
        d.extend(json.loads(f.text)['DS'])
    for i in d:
        for j in s:
            if int(i['Station_Id_C']) == int(j['id']):
                i['longitude'] = j['longitude']
                i['latitude'] = j['latitude']
                break
    df = pd.DataFrame(d)
    x = df['longitude'].values
    y = df['latitude'].values
    z = df[request.GET.get("c")].values

    def plot_contour(x, y, z, resolution=50, contour_method='linear'):
        resolution = str(resolution) + 'j'
        X, Y = np.mgrid[min(x):max(x):complex(resolution), min(y):max(y):complex(resolution)]
        points = [[a, b] for a, b in zip(x, y)]
        Z = griddata(points, z, (X, Y), method=contour_method)
        return X, Y, Z

    X, Y, Z = plot_contour(x, y, z, resolution=50, contour_method='linear')
    locale.setlocale(locale.LC_CTYPE, 'chinese')
    plt.rcParams['font.sans-serif'] = ['SimHei']
    df = gpd.read_file(os.path.dirname(__file__) + '\\gd.json')
    ax = df.plot(figsize=(10, 10), alpha=0.2, edgecolor='k')
    C = ax.contour(X, Y, Z)
    plt.clabel(C, inline=True, fontsize=10)
    plt.xticks(())
    plt.yticks(())
    plt.title(
        "廣東省" + date.strftime("%Y年%m月%d日%H時") + "等" + DCmaDict.objects.filter(key=request.GET.get('c')).first().value+"圖")

    filename = datetime.now().strftime("%Y%m%d%H%M%S") + ".png"
    plt.savefig(os.path.dirname(__file__) + '\\static\\img\\' + filename, bbox_inches='tight')
    return HttpResponse(filename)