捍卫祖国领土从每一张地图开始
地图是表达国家版图最常用、最主要的形式。但在影视剧《亲爱的，热爱的》中出现了明显的错误，从上至下引起了极大的关注度。
地图错用的问题确实很普遍，不仅出现在影视剧中，还经常出现在科学论文中，甚至还会出现在世界顶级期刊中。主要原因是因为网上传播着大量的“问题地图”，大家在关注科学问题时而忽略了地图的重要性，致使大量“问题地图”出现在各大学术期刊中。
前几日自然资源部表示要强化国家版图意识的教育，引导大家使用正确的地图。MeteoAI作为又红又专的中国特色社会主义接班人决定积极响应号召，引导大家正确绘制中国地图。
MeteoAI的小伙伴从资源环境平台下载了中国的行政区划，将其制作成shp文件，并和权威机构的标准地图做了比对，吻合一致。详情见：中国地图的正确打开方式。这回我们就教大家如何用python的地图库Cartopy来加载自制shp文件来绘制出正确的中国地图。
Cartopy
Cartopy是用来绘制地图和地理空间信息分析的python库，主要介绍一下几个主要功能：
1.设置投影2.增加地图特征3.显示自定义shp
设置投影
Cartopy提供了大量的投影方式（Cartopy投影列表[1]），使用cartopy.crs可以调用各个投影，
cartopy.crs.LambertCylindrical #调用兰勃脱投影
cartopy.crs.Mercator #调用麦卡托投影
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
ax = plt.axes(projection=ccrs.PlateCarree()) # 结合matplotlib调用投影
增加地图特征
Cartopy会有一些现成的地理信息，比如：
1.河流 cartopy.feature.RIVERS2.湖泊 cartopy.feature.LAKES3.海岸线 cartopy.feature.COASTLINE4.陆地 cartopy.feature.LAND5.海洋cartopy.feature.OCEAN6.国界线 cartopy.feature.BORDERS
通过add_feature来增加以上的地图信息，这些信息带有三种分辨率，分别为110m，50m和10m。
import cartopy.feature as cfeat
import matplotlib.pyplot as plt
ax.add_feature(cfeat.RIVERS.with_scale('50m')) # 加载分辨率为50的河流
ax.add_feature(cfeat.LAKES.with_scale('110m')) # 加载分辨率为110的湖泊
这些地图特征有些是可以直接用的比如河流和湖泊，但是像国界线就是万万不能用的，因为其中包含的信息是不符合我们国家地图标准的。
显示自定义shp
使用cartopy.io.shapereader中的Reader可以读取shp文件。
from cartopy.io.shapereader import Reader
reader = Reader(your_shp)
再通过cartopy.feature中的ShapelyFeature可以加载自己的shp特征，并设置相关属性。
import cartopy.crs as ccrs
import cartopy.feature as cfeat
proj = ccrs.PlateCarree()
feature = cfeat.ShapelyFeature(reader.geometries(), proj,
edgecolor='k', facecolor=cfeat.COLORS['land'])
最后通过add_feature来增加以上的地图信息
ax.add_feature(feature, linewidth=1)
可视化标准中国地图
如果看完下面的介绍还不能画出正确的中国标准地图的童鞋请疯狂艾特我们。
对于正确的地图绘制我们是认真的，包教包会，一定要手把手教到你完全掌握。
因为中国一点都不能少！！！
import numpy as np
import cartopy.crs as ccrs
import cartopy.feature as cfeat
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
from cartopy.io.shapereader import Reader
import matplotlib.pyplot as plt
import matplotlib.ticker as mtickerimport warnings
def create_map():
extent = [70, 140, 0, 60]
shp_path = r'D:\programfile backup\china_shp\\'
# --创建画图空间
proj = ccrs.PlateCarree() # 创建坐标系
fig = plt.figure(figsize=(6, 8), dpi=350) # 创建页面
ax = fig.subplots(1, 1, subplot_kw={'projection': proj}) # 创建子图
# --设置地图属性
reader = Reader(shp_path + 'Province_9\Province_9.shp')
provinces = cfeat.ShapelyFeature(reader.geometries(), proj, edgecolor='k', facecolor='none')
ax.add_feature(provinces, linewidth=1)
ax.set_extent(extent, crs=proj)
# --增加低分辨率地形图（cartopy自带）
# ax.stock_img() # 增加地形图
# --增加高分辨率地形图（需自行下载）
fname = os.path.join(config["repo_data_dir"], 'raster', 'natural_earth', '10-natural-earth-1.tif')
ax.imshow(imread(fname), origin='upper', transform=proj, extent=[-180, 180, -90, 90])
# --设置网格点属性
gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True, linewidth=1.2, color='k', alpha=0.5, linestyle='--')
gl.xlabels_top = False # 关闭顶端的经纬度标签
gl.ylabels_right = False # 关闭右侧的经纬度标签
gl.xformatter = LONGITUDE_FORMATTER # x轴设为经度的格式
gl.yformatter = LATITUDE_FORMATTER # y轴设为纬度的格式
gl.xlocator = mticker.FixedLocator(np.arange(extent[0], extent[1]+10, 10))
gl.ylocator = mticker.FixedLocator(np.arange(extent[2], extent[3]+10, 10))
return ax
if __name__ == '__main__':
warnings.filterwarnings('ignore')
ax = create_map()
plt.show()
中国台湾岛、钓鱼岛、南海诸岛、藏南地区、阿克赛钦地区、九段线这些典型的易错易少区域都没有问题。
cartopy自带的地形图分辨率不太高，大家如果感兴趣可以去NaturalEarth[2]下载10m分辨率的图片。
最后放上奉上我们的shp文件：链接：https://pan.baidu.com/s/1l0tSO7x3m1V8qlF_21wAaQ 提取码：jadj
如果有小伙伴发现我们制作的shp有不严谨的地方欢迎随时联系我们，后续我们也会对shp文件不断更新，给大家提供更细致更精确的中国行政区划。
以上
References
[1] Cartopy投影列表: https://scitools.org.uk/cartopy/docs/latest/crs/projections.html [2] NaturalEarth: https://www.naturalearthdata.com/downloads/10m-raster-data/
往期回顾
交叉新趋势|采用神经网络与深度学习来预报降水、温度等案例(附代码/数据/文献)