Python空间绘图--Cartopy简介
Python地理信息库包—— Cartopy
一、简介
在前面的教程中，我们已经讲解了常用的二维型数据的可视化方法。但是在日常研究中，由于大气科学属于地学系统，和地球地理信息的结合十分密切，大多数时间，需要在图形中添加地理信息。作为胶水语言，在Python中，目前还在使用的地理可视化库包尚有basemap、cartopy、geopandas等，但由于basemap是基于Python 2，而2已经不再维护，这意味着basemap也要为Python 2陪葬。而geopandas是基于pandas的，属于商务图表利器，但对于气象科研，显得力不从心。现在仅介绍basemap接班者cartopy。
二、Cartopy的安装
在前面已经推荐大家使用conda安装更新库包，从几个Python交流群反馈来看，请尽量不要混用pip与conda来安装库包，反复在pip和conda之间横跳最为致命。如果最开始选择conda，请尽量从一而终。安装Cartopy的命令主要有下面两种：
conda install cartopy
conda install -c conda-forge cartopy
###conda-forge是一个提供库包辅助社区按照之前的教程，请随意新建一个文档，输入import cartopy ,如果没有报错（无事发生*囧*），说明库包已正常安装。
三、初识Cartopy
由于地球是球体，而我们使用的地图是平面的，将球型展开为平面进行绘制时有距离、面积的失真。所以地图学家们提出了各种各样的投影方式，来尽量减小某方面的失真。Cartopy作为专业地理制图库包，提供了非常多的投影方式，能够满足气象业务的需求（import cartopy.crs as ccrs）。其中需要我们重点关注的是:
默认投影（PlateCarree）
兰勃脱投影（Lambert）
墨卡托投影（Mercator）
极投影（我对此类投影的统一称呼，当然学名不叫极投影）
其中，默认投影适合单独省份或者地级市的绘制，这种情况下其变形基本无法看出（内蒙古的meteoer请走开）;兰勃脱投影适合中纬度大范围绘制，比如绘制全中国大公鸡、东亚形势、西北太平洋等；墨卡托投影适合低纬度赤道附近的绘制，一般研究台风、纬向环流等（我也不是专门研究这个的，只能说也许、大概、差不离、估摸着是这样）；极投影适合研究寒潮的北极冷涡。我找到了一张图展示了墨卡托形变（南极离赤道最远，面积变化最大 ）：
当然，默认制图时，中心经线一般为本初子午线。这样中国就偏安东部了，cartopy提供了修改中心经线的命令：
cartopy.crs.PlateCarree(central_longitude=0.0)首先调出默认等距投影platecarree，在内部参数central_longitude处修改到你需要的中心经度。
四、实际操作
千读不如一练，Python气象绘图显然也是如此，下面通过简要的一幅小图，我们来直观感受cartopy的运作。
import matplotlib.pyplot as plt###引入库包
import cartopy.crs as ccrs
proj = ccrs.PlateCarree()
fig = plt.figure(figsize=(4, 4), dpi=550) # 创建画布
ax = fig.subplots(1, 1, subplot_kw={'projection': proj}) # 创建子图
ax.coastlines()##绘制默认海岸线
当然，此处只是绘制了海岸线（最粗糙的那种），一个命令即可绘制全球海岸线，是不是足够酷炫。然而，这样的地图显然无法展示cartopy的地理专业性，所以，我们来尝试添加更多的地理信息：
import cartpy.feature as cfeature
ax.add_feature(cfeature.LAND)####添加陆地######
ax.add_feature(cfeature.COASTLINE,lw=0.3)#####添加海岸线#########
ax.add_feature(cfeature.RIVERS,lw=0.25)#####添加河流######
ax.add_feature(cfeature.LAKES)######添加湖泊#####
ax.add_feature(cfeature.BORDERS, linestyle='-',lw=0.25)####不推荐，我国丢失了藏南、台湾等领土############
ax.add_feature(cfeature.OCEAN)######添加海洋########我们在引入库包时，多引用了一个模块feature(特征、特貌)，并把它缩写为cfeature，通过添加语句即可绘制大部分地理信息：
当然，和我们在前面matplotlib中了解到的一样，所有线条模式都可以改变他的粗细：
ax.add_feature(cfeature.LAND)####添加陆地######
ax.add_feature(cfeature.COASTLINE,lw=1)#####添加海岸线#########
ax.add_feature(cfeature.RIVERS,lw=1)#####添加河流######
ax.add_feature(cfeature.LAKES)######添加湖泊#####
ax.add_feature(cfeature.BORDERS, linestyle='-',lw=1)####不推荐，我国丢失了藏南、台湾等领土############
ax.add_feature(cfeature.OCEAN)######添加海洋########
同样的，所有线状、面状物都是可以赋予颜色的：
ax.add_feature(cfeature.LAND,color='y')####添加陆地######
ax.add_feature(cfeature.COASTLINE,lw=0.25)#####添加海岸线#########
ax.add_feature(cfeature.RIVERS,lw=0.25)#####添加河流######
ax.add_feature(cfeature.LAKES)######添加湖泊#####
ax.add_feature(cfeature.BORDERS, linestyle='-',lw=0.25)####不推荐，我国丢失了藏南、台湾等领土############
ax.add_feature(cfeature.OCEAN,color='g')######添加海洋########
根据前面的介绍，我们使用central_longitude=130，将我国移动至地图中心：
proj = ccrs.PlateCarree(central_longitude=130)##
fig = plt.figure(figsize=(3, 2), dpi=550) # 创建页面
ax = fig.subplots(1, 1, subplot_kw={'projection': proj}) # 创建子图
ax.add_feature(cfeature.LAND)####添加陆地######
ax.add_feature(cfeature.COASTLINE,lw=0.25)#####添加海岸线#########
ax.add_feature(cfeature.RIVERS,lw=0.25)#####添加河流######
ax.add_feature(cfeature.LAKES)######添加湖泊#####
ax.add_feature(cfeature.BORDERS, linestyle='-',lw=0.25)####不推荐，我国丢失了藏南、台湾等领土############
ax.add_feature(cfeature.OCEAN)######添加海洋########
!!!!!!!!!红色警告!!!!!!!!!!
由于该图包的默认命令的参数都是外国人输入的，在绘制国境线时，会有相当多的领土（比如藏南、阿克赛钦、台湾）可能不被画入我国，所以不推荐绘制国境线，必须绘制的情况下，也应规避这些地方或使用我国发布的有效地理信息。（如何绘制完美中国地图会在后续文章推送）
先等一等，这幅地图是不是缺点什么东西。没错，作为专业地图，竟然没有经纬度，这是对cartopy库包的侮辱。
方法一
通过下面一组命令（我是从气象学家上发现的），可以调节出经纬度：
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import matplotlib.ticker as mticker
extent=[-180,180,-90,90]##经纬度范围
gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True, linewidth=0.2, color='k', alpha=0.5, linestyle='--')
gl.xlabels_top = False ##关闭上侧坐标显示
gl.ylabels_right = False ##关闭右侧坐标显示
gl.xformatter = LONGITUDE_FORMATTER ##坐标刻度转换为经纬度样式
gl.yformatter = LATITUDE_FORMATTER
gl.xlocator = mticker.FixedLocator(np.arange(extent[0], extent[1], 30))
gl.ylocator = mticker.FixedLocator(np.arange(extent[2], extent[3], 30))
似乎出了点状况，不要担心，我们能解决的。首先解决刻度重叠的问题，在前面的文章中，我们已经指出——所有刻度类型的问题基本可以通过字典的方式解决，查阅官方文档，我们发现一条命令可以更改刻度字体大小：
gl.xlabel_style={'size':3.5}
gl.ylabel_style={'size':3.5}
经纬网格问题是怎么出现的呢？不知道到大家是否还记得Python切片操作时的一个特点。例如我们对a=[1,2,3,4,5]这个列表进行切片，b=a[0,1]，那么，b里面有几个数字呢？——事实上只有一个1。这条命令表示对a从索引0(事实上是第一个元素)开始取，一直取到索引1（事实上第二个元素），但是索引1是不在被取之列。那么坑爹的地方来了，我们在经纬度范围里设定的是-180~180，-90~90，那么，180和90是不在被取之列的。解决办法即在取值时添加一小段，确保将两端取入。
gl.xlocator = mticker.FixedLocator(np.arange(extent[0], extent[1]+10, 30))
gl.ylocator = mticker.FixedLocator(np.arange(extent[2], extent[3]+10, 30))我们将最后取值加上10个单位，这样就能取到东经190°，北纬100°（当然都是不存在的），又由于刻度间隔为30，这样缺失线就被绘制出来了。
方法二
这个的步骤就比较多了，属于精通matplotlib之后进行魔改。
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
import cartopy.crs as ccrs
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
dlon, dlat = 60, 30
xticks = np.arange(0, 360.1, dlon)
yticks = np.arange(-90, 90.1, dlat)
fig = plt.figure(figsize=(6,5))
ax = fig.add_subplot(1,1,1, projection=ccrs.PlateCarree(central_longitude=180))
ax.coastlines()
gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=False,
linewidth=1, linestyle=':', color='k', alpha=0.8)
gl.xlocator = mticker.FixedLocator(xticks)
gl.ylocator = mticker.FixedLocator(yticks)
ax.set_xticks(xticks, crs=ccrs.PlateCarree())
ax.set_yticks(yticks, crs=ccrs.PlateCarree())
ax.xaxis.set_major_formatter(LongitudeFormatter(zero_direction_label=True))
ax.yaxis.set_major_formatter(LatitudeFormatter())
fig_fname = "全球地图.png"
plt.savefig(fig_fname, dpi=500, bbox_inches='tight')
plt.show()
通过魔改matpltlib出来的地图比气象学家上推送文章更加精美，将南北90度也同时标记出来了。
五、Cartopy小步进阶
大多数时候，中国的气象研究者们可没时间随时随地盯着世界上的天气形势，虽然说全球气象一盘棋，上下游形势互相影响，但是该看本地就看本地，这样才能提高效率。那么，如何看本地呢？
Cartopy提供了extent(范围、长度)命令来调整地图的绘制：
extent=[100,140,30,55]##经纬度范围
ax.set_extent(extent)
完美，中国偏东北地区的区域地图就绘制出来了。不过似乎图形线条都比较僵硬，缺乏美感，也不符合一个专业制图软件的表现。
不要慌，库包开发者们已经提前准备好了，通过提高精度的方式来绘制更加完美的地图。目前cartopy提供了三种精度（10m,50m,110m），这其实是earthnatural的资源。
ax.add_feature(cfeature.OCEAN.with_scale('10m'))
ax.add_feature(cfeature.LAND.with_scale('10m'))
ax.add_feature(cfeature.RIVERS.with_scale('10m'),lw=0.6)
ax.add_feature(cfeature.LAKES.with_scale('10m'))
ax.add_feature(cfeature.BORDERS.with_scale('50m'), linestyle='-',lw=0.6)
ax.add_feature(cfeature.COASTLINE.with_scale('10m'),lw=0.5)
这样是不是就美观细致多了。不过这些资源都需要下载，所以在第一次绘制时候会报错（因为资源未下载），绘制高精度时出图特别慢（本质上是一根线一根线慢慢画，精度越高，转折弯曲越多，绘图越慢），但是绘制成功一次之后，再绘制时就非常快了（因为资源已经下载好了，并保存在本地）。有些教程打包了资源解压到指定文件夹加快绘图速度，我不建议这么做。一是更改目录是一种危险的行为（作为气象本专业的学生，我十分清楚大部分大气科学从业者的计算机水平是不专业的，这是大学专业学得杂，结果样样都会点，样样不精），二是即便是第一次绘制10m精度的图，也只需要大概五分钟即可全部完成，对初学者来说这点时间不值一提。
腾讯云开发者
