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)