双节活动 获奖的书籍会在下 周一(10月12日寄出) ，获奖的同学请耐心等待哦！奖品收到后记得 拍照发我下 哦！

Cartopy进阶——自由的接口

一、复习回顾

在前面一节中，我们已经介绍了cartopy的大致用法——全球地图的绘制、范围的设定以及更改地理信息的精度。但是，有时候这并不能满足我们的需求，比如我作为某地级市的预报员，绘制该市降水图时，为使图片整洁，一般不希望多出其他市县。还有进行地区级别的研究，比如青藏高原地理区划将包含尼泊尔与不丹，cartopy的基础地理信息添加暂时无法做到，但是该库包已经准备了额外的接口以满足这种需求，并且比NCL更加灵活。

二、数据读取接口

Cartopy提供了一个基于pyshp的接口以实现对地理文件的简单读取和操作：

from cartopy.io.shapereader.Reader import Reader

reader（读取器）就是用来读取你想要读取的shp文件。如何使用呢，下面通过两个例子来解释。

首先，引用我们需要的库包：

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，专门的气象地图软件上查看，具体如何操作呢？下面以恩施州分县地图来说明。