Cartopy调用天地图作为底图

MeteoAI

发布于 2019-09-17 15:13:13

3.9K0

发布于 2019-09-17 15:13:13

文章被收录于专栏：MeteoAIMeteoAI

概述

在捍卫祖国领土从每一张地图开始，Python绘制气象实用地图[Code+Data](续)中我们介绍了cartopy这个库，专门用于绘制地图和地理空间信息分析的python库，但是cartopy中的底图都是国外资源，一个是未经审核，另外调用时加载速度也会受影响。本期文章将会介绍如何在cartopy中使用天地图提供的影像，矢量和地形的图层服务。

cartopy调用天地图图层

cartopy自带的底图是有Google、MapQuest、Stamen、Mapbox、Quadtree等多家图层服务，提供影像、矢量和地形图，可以在img_tiles.py文件中查看具体的实现形式。在国内，调用天地图的服务是最具有保证和方便的形式，我们可以通过继承cartopy提供的GoogleWTS类，自定义天地图图层的类，从而使用天地图的服务。

import matplotlib.pyplot as plt
import cartopy.io.img_tiles as cimgt
import cartopy.crs as ccrs
from cartopy.io import shapereader
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
class TDT(cimgt.GoogleWTS):
    def _image_url(self, tile):
        x, y, z = tile
        url = 'http://t3.tianditu.gov.cn/DataServer?T=vec_w&x=%s&y=%s&l=%s&tk=a4ee5c551598a1889adfabff55a5fc27'% (x, y, z)
        return url
class TDT_ter(cimgt.GoogleWTS):
    def _image_url(self, tile):
        x, y, z = tile
        url = 'http://t3.tianditu.gov.cn/DataServer?T=ter_w&x=%s&y=%s&l=%s&tk=a4ee5c551598a1889adfabff55a5fc27'% (x, y, z)
        return url
class TDT_img(cimgt.GoogleWTS):
    def _image_url(self, tile):
        x, y, z = tile
        url = 'http://t3.tianditu.gov.cn/DataServer?T=img_w&x=%s&y=%s&l=%s&tk=a4ee5c551598a1889adfabff55a5fc27'% (x, y, z)
        return url
def make_map(projection=ccrs.PlateCarree()):
    fig, ax = plt.subplots(figsize=(9, 13),                          subplot_kw=dict(projection=projection))
    gl = ax.gridlines(draw_labels=True)
    gl.xlabels_top = gl.ylabels_right = False
    gl.xformatter = LONGITUDE_FORMATTER
    gl.yformatter = LATITUDE_FORMATTER
    return fig, ax
extent = [115, 118, 39, 41.5] #北京市区图
request = cimgt.TDT() #矢量图层
#request = cimgt.TDT_img() #影像
#request = cimgt.TDT_ter() #地形
fig, ax = make_map(projection=request.crs)
ax.set_extent(extent)
ax.add_image(request, 10)# level=10 缩放等级 
plt.show()

矢量图层

影像图层

地形图层

cartopy以天地图为底图画利奇马台风

前段时间的利奇马台风对我国沿海造成了巨大的破坏，我们从中国台风网[1]爬取了利奇马台风的途径数据，利用catopy，以天地图为底图，对利奇马的路径和风力等级进行展示：爬取利奇马台风的数据：

import requests
import json
import time
import pandas as pd
year = "2019" # 输入台风的年份
typhoon_id = "1909" # 输入利奇马台风的编号 1909
url = "http://d1.weather.com.cn/typhoon/typhoon_data/"+year+"/"+typhoon_id+".json?callback=getData&_="+str(int(round((time.time()*1000))))
headers = {
        "User-Agent": "Mozilla/5.0 (Windows NT 10.0; WOW64) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/55.0.2883.87 Safari/537.36",
        "Referer": "http://typhoon.weather.com.cn/gis/typhoon_p.shtml",
    }
r = requests.get(url,headers=headers)
a = json.loads(r.text[8:-1])#解析json文件
with open("./test2.json",'w',encoding='utf-8') as json_file:
    json.dump(a["typhoon"][8],json_file,ensure_ascii=False)
#台风轨迹点
start_time = [] #台风到达时间
lon = [] #台风到达地经度
lat = [] #台风到达地纬度
central_pressure = [] #台风中心气压
wind = [] #台风风速风力（米/秒）
direction = [] #未来移向
feature_speed = [] #未来移速
for i in range(len(a["typhoon"][8])):
    start_time.append(a["typhoon"][8][i][1])
    lon.append(a["typhoon"][8][i][4])
    lat.append(a["typhoon"][8][i][5])
    central_pressure.append(a["typhoon"][8][i][6])
    wind.append(a["typhoon"][8][i][7])
    direction.append(a["typhoon"][8][i][8])
    feature_speed.append(a["typhoon"][8][i][9])
typhoon_data=pd.DataFrame({"start_time":start_time,"lon":lon,"lat":lat,"central_pressure":central_pressure,
                            "wind":wind,"direction":direction,"feature_speed":feature_speed})
typhoon_data["start_time"]=pd.to_datetime(typhoon_data["start_time"],format='%Y-%m-%d %H时')
typhoon_data.to_csv('%s.csv'%typhoon_id,index=False)

获取利奇马的数据文件1909.csv。对台风路径进行绘图展示：

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import cartopy.io.img_tiles as cimgt
import cartopy.feature as cfeat
from cartopy.io.shapereader import Reader
import pandas as pd
import os
from PIL import Image
class TDT(cimgt.GoogleWTS):
    def _image_url(self, tile):
        x, y, z = tile
        url = 'http://t3.tianditu.gov.cn/DataServer?T=vec_w&x=%s&y=%s&l=%s&tk=a4ee5c551598a1889adfabff55a5fc27'% (x, y, z)
        return url
class TDT_ter(cimgt.GoogleWTS):
    def _image_url(self, tile):
        x, y, z = tile
        url = 'http://t3.tianditu.gov.cn/DataServer?T=ter_w&x=%s&y=%s&l=%s&tk=a4ee5c551598a1889adfabff55a5fc27'% (x, y, z)
        return url
class TDT_img(cimgt.GoogleWTS):
    def _image_url(self, tile):
        x, y, z = tile
        url = 'http://t3.tianditu.gov.cn/DataServer?T=img_w&x=%s&y=%s&l=%s&tk=a4ee5c551598a1889adfabff55a5fc27'% (x, y, z)
        return url
def make_map(WTS,level):
    extent = [72, 135, 15, 54]
    shp= 'cartopy/shp/Province_9/Province_9.shp'
    proj = ccrs.PlateCarree()
    fig, ax = plt.subplots(figsize=(10, 6),
                           subplot_kw=dict(projection=WTS.crs))
    gl = ax.gridlines(draw_labels=True)
    gl.xlabels_top = gl.ylabels_right = False
    gl.xformatter = LONGITUDE_FORMATTER
    gl.yformatter = LATITUDE_FORMATTER
    ax.set_extent(extent)
    ax.add_image(WTS, level) 
    reader = Reader(shp)
    provinces = cfeat.ShapelyFeature(reader.geometries(),proj, edgecolor='white', facecolor='none')
    ax.add_feature(provinces, linewidth=0.5) 
    sub_extent = [105, 125, 0, 25]
    sub_ax = fig.add_axes([0.66, 0.11, 0.14, 0.155],
                      projection=WTS.crs)
    sub_ax.set_extent(sub_extent)
    sub_ax.add_image(WTS, level) 
    sub_ax.add_feature(provinces, linewidth=0.3)
    return fig,ax,sub_ax
if __name__ == "__main__":
    typhoon=pd.read_csv("1909.csv")
    typhoon.dropna(subset=['wind'],inplace=True)
    typhoon['class']=''
    typhoon.loc[(typhoon['wind']>=10.8) & (typhoon['wind']<=17.1),'class']='TROPICAL DEPRESSION'
    typhoon.loc[(typhoon['wind']>17.2) & (typhoon['wind']<=24.4),'class']='TROPICAL STORM'
    typhoon.loc[(typhoon['wind']>24.5) & (typhoon['wind']<=32.6),'class']="SEVERE TROPICAL STORM"
    typhoon.loc[(typhoon['wind']>32.7) & (typhoon['wind']<=41.4),'class']="TYPHOON"
    typhoon.loc[(typhoon['wind']>41.5) & (typhoon['wind']<=50.9),'class']="SEVERE TYPHOON"
    typhoon.loc[(typhoon['wind']>51),'class']="SUPER TYPHOON"
    color_dict = {'TROPICAL DEPRESSION': '#7fffd4', 'TROPICAL STORM': '#008000', "SEVERE TROPICAL STORM": '#0000ff', 
    "TYPHOON": '#ffff00', "SEVERE TYPHOON": '#ffa500', "SUPER TYPHOON": '#ff4500'}
    color_dict = {'TROPICAL DEPRESSION': 'lime', 'TROPICAL STORM': 'blue', "SEVERE TROPICAL STORM": 'yellow', 
    "TYPHOON": 'orange', "SEVERE TYPHOON": 'red', "SUPER TYPHOON": 'darkred'}
    wts = TDT_img()  
    for j in typhoon.index[154::1]:
        fig,ax,sub_ax=make_map(wts,5)
        for i in typhoon.index[:j]:
            a0,b0 = typhoon.loc[i,"lon"], typhoon.loc[i,"lat"]
            a1,b1 = typhoon.loc[i+1,"lon"], typhoon.loc[i+1,"lat"]
            ax.plot((a0,a1),(b0,b1), linewidth=2.5,
                    color=color_dict[typhoon.loc[i+1,"class"]], transform=ccrs.Geodetic())
            ax.scatter(a0, b0, s=10, c=color_dict[typhoon.loc[i,"class"]], alpha=0.8, transform=ccrs.Geodetic())
            sub_ax.plot((a0,a1),(b0,b1), linewidth=2.5,
                    color=color_dict[typhoon.loc[i+1,"class"]], transform=ccrs.Geodetic())
            sub_ax.scatter(a0, b0, s=10, c=color_dict[typhoon.loc[i,"class"]], alpha=0.8, transform=ccrs.Geodetic())
        ax.annotate('Typhoon: %s\n Date: %s\n Data: http://typhoon.weather.com.cn/\n' % ('Lekima(ID-1909)', "2019-08-04 14:00:00"),
            xy=(0, 1), xytext=(12, -12), va='top', ha='left',xycoords='axes fraction', textcoords='offset points')
    # plt.show()
        plt.savefig("./pic/%d.png"%j)
    os.chdir("./pic/")
    imgFiles = [fn for fn in os.listdir('.') if fn.endswith('.png')]
    imgFiles.sort(key=lambda x:int(x[:-4]))
    print(imgFiles)
    images = [Image.open(fn) for fn in imgFiles]
    im = images[0]
    filename = 'test.gif'
    im.save(fp=filename, format='gif', save_all=True, append_images=images[1:], duration=500)

利奇马完整的路径

gif动态图

参考文章： python如何调用天地图来绘制底图[2] 飓风“桑迪”路径图的制作[3] 基于Python的GIS分析之台风路径可视化（三）[4]

References

[1] 中国台风网: http://typhoon.weather.com.cn [2] python如何调用天地图来绘制底图: http://bbs.06climate.com/forum.php?mod=viewthread&tid=89165 [3] 飓风“桑迪”路径图的制作: https://www.cnblogs.com/vamei/archive/2012/11/07/2758006.html [4] 基于Python的GIS分析之台风路径可视化（三）: https://beiyuan.me/tytrack-detailed/

本文参与腾讯云自媒体分享计划，分享自微信公众号。

原始发表：2019-09-13，如有侵权请联系 cloudcommunity@tencent.com 删除

python