Python气象数据处理与绘图：Cartopy调用天地图在线服务

气海无涯 2021-07-28

2781

1、前言

我们常用的在线地图，例如搜狗、高德、腾讯等都是进行加密处理的，也就是其坐标系和实际坐标系之间存在偏差，这就导致我们制图的结果也是有偏差的，尤其是在小区域尺度下，那么有没有无偏移的在线地图呢？天地图就可以满足大家的需求，但是调用天地图的服务并不是十分容易，在此给出可行的办法供大家参考。需要注意的是，天地图采用的是CGCS2000的地理坐标系统，而我们的数据通常是WGS1984的地理坐标系统，但两者之间几乎没有差异。

网上也有相关帖子介绍了Cartopy调用天地图的办法，例如：
1、好奇心Log的推送：

1https://my.oschina.net/u/4581316/blog/4396806

2、气象家园论坛：

1http://bbs.06climate.com/forum.php?mod=viewthread&tid=89165

但是由于天地图服务升级，原先的方法都不再适用，这里给出的是最新的调用方法.

2、注册天地图

国家地理信息公共服务平台：

1https://www.tianditu.gov.cn/

注册并申请地图API的秘钥，选择服务器端的key。

3、修改img_tiles.py文件

在个人电脑里，找到

1Anaconda3\Lib\site-packages\cartopy\io\img_tiles.py

并在最后添加以下代码，填入之前申请到的key。

在KesciLab里，img_tiles.py 文件所在目录是

1/opt/conda/lib/python3.7/site-packages/cartopy/io/img_tiles.py

（可以通过 !cd /opt/conda/ && grep -rl cartopy . 语句查到）

可以通过新建一个 py 文件，比如 key.py，然后用 cat 追加到 img_tiles.py 里。代码是：

1!cat /home/kesci/work/key/key.py >> /opt/conda/lib/python3.7/site-packages/cartopy/io/img_tiles.py

推荐还是在本地的电脑环境进行配置！！！

 1# 天地图矢量
 2class TDT_vec(GoogleWTS):
 3    def _image_url(self, tile):
 4        x, y, z = tile
 5        key = 'your_key'
 6        url = 'http://t0.tianditu.gov.cn/DataServer?T=vec_w&x=%s&y=%s&l=%s&tk=%s' % (x, y, z, key)
 7        return url
 8
 9# 天地图遥感
10class TDT_img(GoogleWTS):
11    def _image_url(self, tile):
12        x, y, z = tile
13        key = 'your_key'
14        url = 'http://t0.tianditu.gov.cn/DataServer?T=img_w&x=%s&y=%s&l=%s&tk=%s' % (x, y, z, key)
15        return url
16
17# 天地图地形
18class TDT_ter(GoogleWTS):
19    def _image_url(self, tile):
20        x, y, z = tile
21        key = 'your_key'
22        url = 'http://t0.tianditu.gov.cn/DataServer?T=ter_w&x=%s&y=%s&l=%s&tk=%s' % (x, y, z, key)
23        return url

4、调用矢量底图、影像底图、地形底图

 1# 导入模块
 2import matplotlib.pyplot as plt
 3import cartopy.crs as ccrs
 4import cartopy.feature as cfeat
 5from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
 6from cartopy.mpl.ticker import LongitudeFormatter,LatitudeFormatter
 7import cartopy.io.shapereader as shpreader
 8import cartopy.io.img_tiles as cimgt
 9import numpy as np
10
11plt.rcParams['font.sans-serif'] = ['SimHei'] 
12plt.rcParams['axes.unicode_minus'] = False

绘图

 1fig = plt.figure(figsize=(18, 12))
 2
 3ax = fig.add_subplot(1, 3, 1, projection=ccrs.PlateCarree())
 4ax.set_extent([118, 122, 28, 32],crs=ccrs.PlateCarree())
 5request = cimgt.TDT_vec()
 6ax.add_image(request, 9)
 7ax.set_title('天地图矢量底图',fontsize=15)
 8gl = ax.gridlines(draw_labels=True, linewidth=1, color='k', alpha=0.5, linestyle='--')
 9gl.xlabels_top = gl.ylabels_right = False 
10gl.xformatter = LONGITUDE_FORMATTER 
11gl.yformatter = LATITUDE_FORMATTER
12
13ax = fig.add_subplot(1, 3, 2, projection=ccrs.PlateCarree())
14ax.set_extent([118, 122, 28, 32],crs=ccrs.PlateCarree())
15request = cimgt.TDT_img()
16ax.add_image(request, 9)
17ax.set_title('天地图影像底图',fontsize=15)
18gl = ax.gridlines(draw_labels=True, linewidth=1, color='k', alpha=0.5, linestyle='--')
19gl.xlabels_top = gl.ylabels_right = False 
20gl.xformatter = LONGITUDE_FORMATTER 
21gl.yformatter = LATITUDE_FORMATTER
22
23ax = fig.add_subplot(1, 3, 3, projection=ccrs.PlateCarree())
24ax.set_extent([118, 122, 28, 32],crs=ccrs.PlateCarree())
25request = cimgt.TDT_ter()
26ax.add_image(request, 9)
27ax.set_title('天地图地形底图',fontsize=15)
28gl = ax.gridlines(draw_labels=True, linewidth=1, color='k', alpha=0.5, linestyle='--')
29gl.xlabels_top = gl.ylabels_right = False 
30gl.xformatter = LONGITUDE_FORMATTER 
31gl.yformatter = LATITUDE_FORMATTER

5、比较不同地图层级

添加的在线地图是以瓦片形式存储的，不同的层级代表不同的分辨率，可以根据研究区域的大小选择合适的层级，注意过高的分辨率加载会比较慢。

 1fig = plt.figure(figsize=(18, 12))
 2for i in range(6,9):
 3    ax = fig.add_subplot(1, 3, i-5, projection=ccrs.PlateCarree())
 4    ax.set_extent([120.4, 122.1, 30.6, 32.1],crs=ccrs.PlateCarree())
 5    request = cimgt.TDT_img()
 6    ax.add_image(request, i)
 7    ax.set_title('Level='+str(i),fontsize=15)
 8    gl = ax.gridlines(xlocs=np.arange(120.5, 122.5, 0.5),
 9                      ylocs=np.arange(30, 32.5, 0.5),
10                      draw_labels=True,linewidth = 0.5,color='k',
11                      alpha=0.5,linestyle='--')
12    gl.xlabels_top = gl.ylabels_right = False

有问题可以到QQ群里进行讨论，我们在那边等大家。

QQ群号：854684131

python

文章转载自气海无涯，如果涉嫌侵权，请发送邮件至：contact@modb.pro进行举报，并提供相关证据，一经查实，墨天轮将立刻删除相关内容。