Cartopy调用天地图作为底图
- 2019 年 10 月 4 日
- 筆記
概述
在捍卫祖国领土从每一张地图开始，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: %sn Date: %sn 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/
