用10行代码调用网络地图服务

点击图片上方蓝色字体“慧天地”即可订阅

文章转载自微信公众号数读城事，作者，陈玮，排版，数读菌，校阅，数读菌，版权归原作者及刊载媒体所有。

前言：

Cartopy是一个Python包，用于地理空间数据处理，以便生成地图和其他地理空间数据分析。

01Why

Reason1: Cartopy调用在线地图服务的教程很多，例如公众号气海无涯的教程，洋流大佬以天地图为例，详细介绍了如何用Cartopy调用自定义网络地图服务并绘制，但没有对Cartopy各方法的关键参数进行讲解，也没有介绍如何在底图的基础上叠加路网和标注图层，所以我在洋流大佬教程的基础上做一点补充。大佬原文：
Python气象数据处理与绘图：Cartopy调用天地图在线服务
Reason 2: 天地图提供的地形图渲染简直太丑了，下面两图分别用天地图地形和Stamen绘制地形，Stamen不仅在有水系、聚居区，地形渲染的颜色和山体阴影也非常合理。

Reason 3：由于国外网络地图一般在中国边界上都存在不同的问题，且标注信息多为英文，所以只能通过自己叠加国内在线地图服务添加边界、路网、地名标注等信息。

02

How

先上示意图，使用Staman无标注地形图、卫星影像和高德地图路网、边界、地名标注绘制。

图1｜地形图

                                               图2｜卫星影像

代码解释：

# 导入相关库
from matplotlib import pyplot as plt
import gdal
from pylab import mpl
import cartopy.crs as ccrs
import cartopy.io.img_tiles as cimgt

Cartopy的cimgt专门用来调用在线地图服务并绘制在Matplotlib中。cimgt中虽然自带了诸如谷歌地图、谷歌卫星影像、Stamen等常用在线地图服务，但我们也可以通过修改img_tiles.py文件添加天地图、高德地图等自定义网络地图服务。修改方法见洋流大佬的教程。

# 首先定义四个方向的边界，即要绘制的WGS84投影经纬度坐标
west = 72
east = 136
south = 2.5
north = 54

# 获取地形图层，此处采用Stamen的无标注地形
terrain = cimgt.Stamen('terrain-background')
# 获取地形图层，此处采用高德路网和边界，需要自行添加
label = cimgt.GD_road_l()

fig = plt.figure(figsize=(20, 20))
# 此处将Axes投影系设置为网络地图使用的莫卡托投影，否则网络地图会因投影不正确而变形
ax = fig.add_subplot(projection=ccrs.GOOGLE_MERCATOR)

# 将Axes绘图边界设置为我们开始定义的东、西、南、北坐标
# 此处crs=ccrs.PlateCarree()参数表示原始输入数据投影为WGS84
# Cartopy将自动转换为之前设置的GOOGLE_MERCATOR投影
ax.set_extent([west, east, south, north], crs=ccrs.PlateCarree())

# 在Axes上添加获取到的网络地图图层
# 此处第二个数字参数表示网络地图的缩放级别，一般为0-18
ax.add_image(terrain, 5)
ax.add_image(label, 5)

# 保存图像，bbox_inches='tight'参数可以去除多余白边
plt.savefig('China.png', bbox_inches='tight')

由于Cartopy的cimgt中没有高德地图，所以这里需要仿照洋流大佬的教程，在img_tiles.py文件中的添加高德路网在线地图服务。

Tips：在Pycharm中按住Ctrl键并单击代码中的cimgt即可打开img_tiles.py文件

# 高德路网（含注记）
class GD_road_l(GoogleWTS):
def _image_url(self, tile):
        x, y, z = tile
        url = 'https://wprd02.is.autonavi.com/appmaptile?x=%s&y=%s&z=%s&lang=zh_cn&size=1&scl=1&style=8' % (x, y, z)
return url

# 高德路网（不含注记）
class GD_road(GoogleWTS):
def _image_url(self, tile):
        x, y, z = tile
        url = 'https://wprd02.is.autonavi.com/appmaptile?x=%s&y=%s&z=%s&lang=zh_cn&size=1&scl=2&style=8' % (x, y, z)
return url

由于边界、路网以及标注图层是透明图层，所以我们需要将img_tiles.py文件中第55行desired_tile_form='RGB'字段修改为desired_tile_form='RGBA'，添加Alpha通道，这样才能使cimgt支持透明图层。

class GoogleWTS(six.with_metaclass(ABCMeta, object)):
"""
    Implement web tile retrieval using the Google WTS coordinate system.

    A "tile" in this class refers to the coordinates (x, y, z).

    """
    _MAX_THREADS = 24

def __init__(self,desired_tile_form='RGBA',
user_agent='CartoPy/' + cartopy.__version__):
self.imgs =[]
self.crs = ccrs.Mercator.GOOGLE
self.desired_tile_form =desired_tile_form
self.user_agent =user_agent
# some providers like osm need a user_agent in the request issue #1341
# osm may reject requests if there are too many of them, in which case
# a change of user_agent may fix the issue.

03Showcase
下面展示通过上述方法绘制的地图，其中包含shp、geojson矢量图形和tif栅格数据，并将地名标注绘制在最上层。

内容转载、商务活动、投稿等合作请联系

微信号：huitiandi321

邮箱：[e-mail]

欢迎关注慧天地同名新浪微博：

ID：慧天地_geomaticser

往期精彩推荐

《慧天地》敬告

《慧天地》公众号聚焦国内外时空信息科技前沿、行业发展动态、跨界融合趋势，探索企业核心竞争力，传播测绘地理信息文化，为测绘、地信、遥感等相关专业的同学提供日常学习、考研就业一站式服务，旨在打造政产学研用精准对接的平台。《慧天地》高度重视版权，对于原创、委托发布的稿件，会烦请作者、委托方亲自审核通过后才正式推发；对于来自网站、期刊、书籍、微博、微信公众号等媒介的稿件，会在作者栏或者文章开头显著标明出处，以表达对作者和推文引用平台版权的充分尊重和感谢；对于来源于网络作者不明的优质作品，转载时如出现侵权，请后台留言，我们会及时删除。感谢大家一直以来对《慧天地》的关注和支持！——《慧天地》运营团队

编辑：朱奔宇 审核：马冉冉

指导：万剑华教授