Cartopy调用天地图作为底图
概述
在 捍卫祖国领土从每一张地图开始 ， 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]