python中reader_Python气象绘图教程（八）——Cartopy_2

Cartopy进阶——自由的接口

一、复习回顾

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

二、数据读取接口

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

from cartopy.io.shapereader.Reader import Reader

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

例1

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

import numpy as npimport cartopy.crs as ccrsimport cartopy.feature as cfeatfrom cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTERfrom cartopy.io.shapereader import Readerimport matplotlib.pyplot as pltimport 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)

这样即可灵活实现地图的绘制，满足我们日常的需求。

三、下期预告

下期推送将讲解cartopy和填色图。

往期精彩：

Python气象绘图教程(一)

Python气象绘图教程(二)

Python气象绘图教程(三)

Python气象绘图教程(四)

Python气象绘图教程(五)

Python气象绘图教程(六)

Python气象绘图教程(七)——Cartopy

Python气象绘图教程特刊(一)

欢迎关注云台书使获取更多资讯