Python空间绘图绘图——Cartopy 进阶
双节活动获奖的书籍会在下周一(10月12日寄出)，获奖的同学请耐心等待哦！奖品收到后记得拍照发我下哦！
Cartopy进阶——自由的接口
一、复习回顾
在前面一节中，我们已经介绍了cartopy的大致用法——全球地图的绘制、范围的设定以及更改地理信息的精度。但是，有时候这并不能满足我们的需求，比如我作为某地级市的预报员，绘制该市降水图时，为使图片整洁，一般不希望多出其他市县。还有进行地区级别的研究，比如青藏高原地理区划将包含尼泊尔与不丹，cartopy的基础地理信息添加暂时无法做到，但是该库包已经准备了额外的接口以满足这种需求，并且比NCL更加灵活。
二、数据读取接口
Cartopy提供了一个基于pyshp的接口以实现对地理文件的简单读取和操作：
from cartopy.io.shapereader.Reader import Readerreader（读取器）就是用来读取你想要读取的shp文件。如何使用呢，下面通过两个例子来解释。
例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 mticker然后，以shp_path=r"文件路径"存储好你存放shp文件的地址。为什么需要加一个r在前面呢？这是因为电脑比较笨，不给绝对地址，他可能找不到。
shp_path=r'E:\enshi\恩施.shp'#确定shp文件地址接着，按照前面教的绘图流程应该添加画布，增加子图，准备绘制。
proj= ccrs.PlateCarree() # 简写投影
fig = plt.figure(figsize=(4, 4), dpi=400) # 创建画布
ax = fig.subplots(1, 1, subplot_kw={'projection': proj}) # 创建子图准备工作做好之后，就可以使用Reader来读取你的shp文件，并通过cartopy.feature中的ShapelyFeature添加shp特征：
extent=[108.2,110.8,29.1,31.401]#限定绘图范围
reader = Reader(shp_path)
enshicity = cfeat.ShapelyFeature(reader.geometries(), proj, edgecolor='k', facecolor='none')
ax.add_feature(enshicity, linewidth=0.7)#添加市界细节
ax.set_extent(extent, crs=proj)通过plt.show()语句展示绘制出来的图像：
在添加地理信息时，还有两个参数——edgecolor和facecolor，这两个参数直接控制展示出来的图形框线和填充颜色：
enshicity = cfeat.ShapelyFeature(reader.geometries(), proj, edgecolor='g', facecolor='r')
按照之前的讲述的规则，线条的粗细也是可以改变的：
ax.add_feature(enshicity, linewidth=3)#添加市界细节
最后，通过gridliner语句，完善经纬度信息：
gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True, linewidth=0.5, color='r', 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]+0.5, 0.4))
gl.ylocator = mticker.FixedLocator(np.arange(extent[2], extent[3]+0.5, 0.4))
gl.xlabel_style={'size':7}
gl.ylabel_style={'size':7}
不过，因为当时我是单独绘制的，所以市县信息是查出来，单独存放在字典里，然后通过scatter和text命令将站点打印出来：
nameandstation={"恩施":[109.5,30.2],"利川":[109,30.3],"巴东":[110.34,31.04],"建始":[109.72,30.6],"宣恩":[109.49,29.987],"来凤":[109.407,29.493],"咸丰":[109.14,29.665],"鹤峰":[110.034,29.89]}
for key,value in nameandstation.items():
ax.scatter(value[0] , value[1] , marker='.' , s=60 , color = "k" , zorder = 3)
ax.text(value[0]-0.07 , value[1]+0.03 , key , fontsize = 7 , color = "k")
我是通过key存放市县名，value存放经纬度信息。然后通过for in 遍历字典绘制站点。这算是我在两个月之前刚学习时想出的笨办法，如果读者有更方便的办法，可在后台留言交流。
例2
通过add_geometries()添加地理信息：
shp_path=r'E:\shp\行政边界.shp'
a_shapes=list(Reader(shp_path).geometries())
ax.add_geometries(a_shapes[:],crs=proj,edgecolor='k',facecolor='',lw=0.75)
这种绘图方式有什么用处呢？关键在list列表里，现在我们更改列表的切片范围：
ax.add_geometries(a_shapes[2:9],crs=proj,edgecolor='k',facecolor='',lw=0.75)上一步是[:]表示从头到尾全部取完，现在我们取[2:9]
amazing!!!!显然，这样使我们的绘制灵活度上升了不少。我们可以看看[2:9]切片应该有多少县呢？从索引2开始，2、3、4、5、6、7、8，应该有七个县，绘制的县有多少呢？也是七个。这样即明白地展示其原理。我们还可以只绘制一个县：
ax.add_geometries(a_shapes[:1],crs=proj,edgecolor='k',facecolor='',lw=0.6)
如何知道每个县对应的列表索引呢？在几何图形比较少的情况下（<10），大可以逐个实验，对列表单独切片。另外的利器有meteoinfo，专门的气象地图软件上查看，具体如何操作呢？下面以恩施州分县地图来说明。
shp_path=r'E:\shp\恩施土家族苗族自治州_行政边界\恩施土家族苗族自治州_行政边界.shp'
a_shapes=list(Reader(shp_path).geometries())
ax.add_geometries(a_shapes[:],crs=proj,edgecolor='k',facecolor='',lw=0.6)
现在是从头至尾全部绘制，然后我们按照在Python气象绘图教程特刊（一）中的方法，查出图层属性：
我们可以看出，第一个是鹤峰县，那么我们使切片变为[:1]并绘制：
ax.add_geometries(a_shapes[:1],crs=proj,edgecolor='k',facecolor='',lw=0.6)
的确是鹤峰县被单独绘制出来了，现在试着取来凤县和利川市[1:3]
ax.add_geometries(a_shapes[1:3],crs=proj,edgecolor='k',facecolor='',lw=0.6)
这样即可灵活实现地图的绘制，满足我们日常的需求。
腾讯云开发者
