Python气象绘图教程（十三）—Cartopy_4

气象学家

发布于 2020-06-17 17:23:14

8.2K1

发布于 2020-06-17 17:23:14

本节提要：关于子图的一些问题、使用path添加示意框线、Cartopy台风实例本土化

一、关于子图的一些问题

在某些时候，我们需要展示某个地区在整个地图中的位置，常规的方法是绘制两幅地图，比如一张为全国地图，一张为湖北省地图。常见的subplot和subplot2grid函数一般来说绘制的地图大小是一样的，不容易展示比例大小，所以我们选择add_axes()命令来绘制两个大小不一样的子图。

ax1=fig.add_axes([x1,y1,m1,n1])

ax2=fig.add_axes([x2,y2,m2,n2]，projection=ccrs.PlateCarree())

上面这两个命令，即添加了两个子图。在前面已经讲解了x,y,m,n具体意义，此不赘述。唯有一点希望读者注意，在此时ax1与ax2已经不是一类子图了，因为ax2在使用了projection命令之后，已经转变为cartopy中的GeoAxes。而ax1还是matplotlib中的Axes，在大部分时候，两个容器的命令是一致的，也有一定的例外，这种例外会导致Axes的命令不能对GeoAxes生效，但也不会报错。具体如何解决，到时候总结完成再说。

还有一点，添加了投影的ax2与没投影的ax1就算设定子图大小一致，画出来也是不一样的，这可能涉及投影转换问题。比如：

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cf
fig=plt.figure(figsize=(3,3),dpi=300)
ax1=fig.add_axes([0,0,0.4,0.4])
ax2=fig.add_axes([0.5,0,0.4,0.4],projection=ccrs.PlateCarree())
ax2.coastlines()

设置两子图的长宽都为相对于画布的0.4个单位，但是经过投影转换后的ax2是无法填充满子图，而会产生空白，导致视觉上显得ax2比较小，其实两子图大小一致：

第二小节，就是GeoAxes与Axes不同的第一个点，在第十一章中，我们设置了两个命令以改变边框：

ax.spines['right'].set_color('None')
ax.spines['bottom'].set_linewidth(5)

原子图为Axes，前一个是设置框线颜色，后一个为设置框线粗细。但是当我们在GeoAxes中使用这个命令时，他是无效的。（然而并不会出现报错的情况），这是因为在GeoAxes中，设置了新的命令来管理框线——ax.outline_patch。

这时你可以重新设置框线了：

ax1.outline_patch.set_linewidth(5)

另一种方式是关闭GeoAxes的框线，然后开启Axes的框线，对Axes框线命令进行操作：

ax1.outline_patch.set_visible(False)#关闭GeoAxes框线
ax1.spines['left'].set_visible(True)#开启Axes框线
ax1.spines['bottom'].set_visible(True)#开启Axes框线
ax1.spines['right'].set_visible(True)#开启Axes框线
ax1.spines['top'].set_visible(True)#开启Axes框线
ax1.spines['bottom'].set_linewidth(5)#设置Axes框线宽度
ax1.spines['left'].set_linewidth(5)

两种方式均可在视觉上修改框线，但也可以看出一点小区别，比如在图的左下角，修改GeoAxes的ax.outline_patch的命令出现了一个小缺失，而使用Axes的命令时是完全的。

第三小节，介绍如何在子图间添加连接线。

首先，我们绘制出两个子图（添加地理信息和填色图命令已略去）：

ax1= fig.add_axes([0,0.1,0.8,0.8],projection=proj)
ax2= fig.add_axes([0.3,0,0.4,0.4],projection=proj)

然后导入添加连线的模块，添加绘制连线的命令：

from mpl_toolkits.axes_grid1.inset_locator import mark_inset
mark_inset(ax1,ax2,loc1=1,loc2=3,alpha=0.5,fc='None',ec='r',ls='-.')

其中，ax1,ax2代表需要产生连接的子图，loc1，loc2表示需要产生连线的结点，fc即facecolor，ec表示 edgecolor，其他参数与线装参数类似：

现在修改某些参数以展示各关键字意义：

mark_inset(ax1,ax2,loc1=1,loc2=2,alpha=0.5,fc='yellow',ec='g',ls='-')

最后，我们还可以通过在上一节中提到的边框命令解除湖北省的框线：

ax1.background_patch.set_visible(False)
ax1.outline_patch.set_visible(False)

在这种情况下，上图产生的红框和放大子图的黑框大小是等比例放大缩小的。

二、使用path添加框线

在某些时候，需要在一幅地图中框选出比较重要的区域，很多同学使用plt.plot()命令绘制，是比较简便的。但是，也需要介绍比较复杂的命令path，这个命令会在某些时候与set boundary相连接，而plot命令是不能设置边界的。

vertices = []
codes = []

codes = [Path.MOVETO] + [Path.LINETO]*3 + [Path.CLOSEPOLY]
vertices = [(1, 1), (1, 50), (50, 50), (50, 1), (0, 0)]
vertices = np.array(vertices, float)
path = Path(vertices, codes)
pathpatch = PathPatch(path, facecolor='None', edgecolor='green')
fig=plt.figure(figsize=(4,4),dpi=500)
ax=fig.add_axes([0,0,1,1],projection=ccrs.PlateCarree())
ax.add_feature(cfeat.OCEAN)
ax.add_feature(cfeat.LAND)
ax.set_extent([0,100,0,50])
ax.add_patch(pathpatch)

codes代表每一步参考的命令，vertices代表画点坐标。具体参考官网手册。

三、Cartopy台风实例本土化

在官网上，有一个台风行径影响图，我在学习时将其中国化并实用化，添加了一些中文注释，并增加了一个本土化例子，希望能帮助到学习的同学：

import numpy as np
import cartopy.crs as ccrs
import cartopy.feature as cfeat
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import cartopy.io.shapereader as shpreader
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
import matplotlib.path as mpath
import matplotlib.patches as mpatches
import shapely.geometry as sgeom
plt.rcParams['font.sans-serif']=['SimHei']
###虚构台风路径数据###################################################
def sample_data(): 
    lons=[125,122,120.5,120.1,120,119,118,116,114.4]
    lats=[15,20.5,23,24.5,27,30.4,33,35.5,36]
    return lons,lats
###自定义台风符号#######################################################
def get_hurricane():
    u = np.array([  [2.444,7.553],
                    [0.513,7.046],
                    [-1.243,5.433],
                    [-2.353,2.975],
                    [-2.578,0.092],
                    [-2.075,-1.795],
                    [-0.336,-2.870],
                    [2.609,-2.016]  ])
    u[:,0] -= 0.098
    codes = [1] + [2]*(len(u)-2) + [2] 
    u = np.append(u, -u[::-1], axis=0)
    codes += codes
    return mpath.Path(3*u, codes, closed=False)
def main():
###准备之后要用的常规数据############################################
    extent = [110, 140, 10, 60]#定义绘图范围
    shp_path = r'D:\python\Province_9\Province_9.shp'
    proj = ccrs.PlateCarree() #缩写投影
    hurricane = get_hurricane()#获得台风符号
    lons,lats=sample_data()#获得台风经纬数据
    fig=plt.figure(figsize=(5,6),dpi=600)#添加子图
    ax = fig.add_axes([0, 0, 1, 1], projection=proj,
                      frameon=False)
    ax.set_extent(extent, proj)
    provinces_shp=r'D:\python\Province_9\Province_9.shp'
    ax.set_title('一次台风路径及影响省份')
    track = sgeom.LineString(zip(lons, lats))#将台风线条转化为地理几何线条集合，Zip是Python基础库命令
    track_buffer = track.buffer(2)#缓冲两度，经纬两度
    def colorize_state(geometry):
        facecolor = (0.9375, 0.9375, 0.859375)
        if geometry.intersects(track):#intersects命令判断我们的台风路径线条是否与这个省的封闭几何相交
            facecolor = 'red'
        elif geometry.intersects(track_buffer):
            facecolor = '#FF7E00'
        return {'facecolor': facecolor, 'edgecolor': 'black'}
    ax.add_geometries(shpreader.Reader(provinces_shp).geometries(),ccrs.PlateCarree(),lw=0.5,styler=colorize_state)
    ###这一步，使各省着色###
    ax.add_geometries([track_buffer], ccrs.PlateCarree(),
                      facecolor='#C8A2C8', alpha=0.5)
    ###添加浅红色外围###
    ax.add_geometries([track], ccrs.PlateCarree(),
                      facecolor='none', edgecolor='k')
    ###添加台风路径线###
    ax.scatter(lons,lats, s=150, marker=hurricane, 
            edgecolors="red", facecolors='none', linewidth=2.3,zorder=3)#添加台风符号
    ############添加其他信息##############
    ax.add_feature(cfeat.COASTLINE.with_scale('50m'),lw=0.7)
    ax.add_feature(cfeat.LAND.with_scale('50m'))
    ax.add_feature(cfeat.OCEAN.with_scale('50m'))
    gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True, linewidth=0.7, 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))
    #######添加图例#######
    direct_hit = mpatches.Rectangle((0, 0), 1, 1, facecolor="red")
    within_2_deg = mpatches.Rectangle((0, 0), 1, 1, facecolor="#FF7E00")
    labels = ['正面经过',
              '外围影响']
    ax.legend([direct_hit, within_2_deg], labels,
              loc='lower left', bbox_to_anchor=(0.6, 0.05), fancybox=True)
    fig.savefig('台风路径',bbox_inches='tight')
    plt.show()
if __name__ == '__main__':
    main()

本土化的例子：

import numpy as np
import cartopy.crs as ccrs
import cartopy.feature as cfeat
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import cartopy.io.shapereader as shpreader
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
import matplotlib.path as mpath
import matplotlib.patches as mpatches
import shapely.geometry as sgeom
plt.rcParams['font.sans-serif']=['SimHei']
###虚构中尺度路径数据###################################################
def sample_data(): 
    lons=[108.5,108.8,109,109.5,109.8,109.9]
    lats=[31.5,31.5,31.5,31.5,31,31]
    return lons,lats
def main():
    extent=[107.5,111.5,28.5,32]#定义绘图范围
    shp_path = r'E:\shp\行政边界.shp'
    proj = ccrs.PlateCarree() #创建坐标系
    lons,lats=sample_data()
    fig=plt.figure(figsize=(5,5),dpi=600)
    ax = fig.add_axes([0, 0, 1, 1], projection=proj,
                      frameon=False)
    ax.set_extent(extent, proj)
    cities_shp=r'E:\shp\行政边界.shp'
    ax.set_title('一次中尺度对流影响的县市',fontsize=20)
    track = sgeom.LineString(zip(lons, lats))#将台风线条转化为地理信息
    track_buffer = track.buffer(1)#缓冲1度
    track_buffer2=track.buffer(2)#缓冲2度
    def colorize_state(geometry):
        facecolor = (0.9375, 0.9375, 0.859375)
        if geometry.intersects(track):
            facecolor = 'red'
        elif geometry.intersects(track_buffer):
            facecolor = '#FF7E00'
        elif geometry.intersects(track_buffer2):
            facecolor='green'
        return {'facecolor': facecolor, 'edgecolor': 'black'}
    ax.add_geometries(shpreader.Reader(cities_shp).geometries(),ccrs.PlateCarree(),lw=0.5,styler=colorize_state)
    ax.add_geometries([track_buffer2], ccrs.PlateCarree(),
                      facecolor='cyan', alpha=0.2)
    ax.add_geometries([track_buffer], ccrs.PlateCarree(),
                      facecolor='#C8A2C8', alpha=0.5)
    ax.add_geometries([track], ccrs.PlateCarree(),
                      facecolor='none', edgecolor='k')
    ax.scatter(lons,lats, s=150, marker='.', 
            edgecolors="red", facecolors='none', linewidth=2.3,zorder=3)
    ############添加其他信息##############
    ax.add_feature(cfeat.RIVERS.with_scale('10m'))
    gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True, linewidth=0.7, 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]+2, 1))
    gl.ylocator = mticker.FixedLocator(np.arange(extent[2], extent[3]+2, 1))
    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=65 , color = "k" , zorder = 3)
        ax.text(value[0]-0.09 , value[1]+0.05 , key , fontsize = 9 , color = "k")
    ###########################################
    direct_hit = mpatches.Rectangle((0, 0), 1, 1, facecolor="red")
    within_1_deg = mpatches.Rectangle((0, 0), 1, 1, facecolor="#FF7E00")
    within_2_deg = mpatches.Rectangle((0, 0), 1, 1, facecolor="green")
    labels = ['强对流过境',
              '风雨影响','外围缓和区']
    ax.legend([direct_hit, within_1_deg,within_2_deg], labels,
              loc='lower left', bbox_to_anchor=(0.7, 0.04), fancybox=True)
    fig.savefig('一次中尺度对流影响的县市',bbox_inches='tight')
    plt.show()
if __name__ == '__main__':
    main()

在程序中，有两个命令比较关键，一个是buffer。在第二张图中，设置了在黑线相交的地市填充为红色，在一个经纬度影响范围内填为橙色，在两个经纬度范围内填成绿色。另一个即track = sgeom.LineString(zip(lons, lats))一句，使得台风路径变为几何信息，这根黑线和地图线是一类，而不是plt.plot()这样的线条。

下期预告：Legend与Colorbar等

本文参与腾讯云自媒体分享计划，分享自微信公众号。

原始发表：2020-06-11，如有侵权请联系 cloudcommunity@tencent.com 删除

add