本节提要：关于子图的一些问题、使用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就算设定子图大小一致，画出来也是不一样的，这可能涉及投影转换问题。比如：
代码语言：javascript复制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不同的第一个点，在第十一章中，我们设置了两个命令以改变边框：
代码语言：javascript复制ax.spines['right'].set_color('None')
ax.spines['bottom'].set_linewidth(5)原子图为Axes，前一个是设置框线颜色，后一个为设置框线粗细。但是当我们在GeoAxes中使用这个命令时，他是无效的。（然而并不会出现报错的情况），这是因为在GeoAxes中，设置了新的命令来管理框线——ax.outline_patch。
这时你可以重新设置框线了：
代码语言：javascript复制ax1.outline_patch.set_linewidth(5)
另一种方式是关闭GeoAxes的框线，然后开启Axes的框线，对Axes框线命令进行操作：
代码语言：javascript复制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的命令时是完全的。
第三小节，介绍如何在子图间添加连接线。
首先，我们绘制出两个子图（添加地理信息和填色图命令已略去）：
代码语言：javascript复制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)
然后导入添加连线的模块，添加绘制连线的命令：
代码语言：javascript复制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，其他参数与线装参数类似：
现在修改某些参数以展示各关键字意义：
代码语言：javascript复制mark_inset(ax1,ax2,loc1=1,loc2=2,alpha=0.5,fc='yellow',ec='g',ls='-')
最后，我们还可以通过在上一节中提到的边框命令解除湖北省的框线：
代码语言：javascript复制ax1.background_patch.set_visible(False)
ax1.outline_patch.set_visible(False)
在这种情况下，上图产生的红框和放大子图的黑框大小是等比例放大缩小的。
二、使用path添加框线
在某些时候，需要在一幅地图中框选出比较重要的区域，很多同学使用plt.plot()命令绘制，是比较简便的。但是，也需要介绍比较复杂的命令path，这个命令会在某些时候与set boundary相连接，而plot命令是不能设置边界的。
代码语言：javascript复制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台风实例本土化
在官网上，有一个台风行径影响图，我在学习时将其中国化并实用化，添加了一些中文注释，并增加了一个本土化例子，希望能帮助到学习的同学：
代码语言：javascript复制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:pythonProvince_9Province_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:pythonProvince_9Province_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()
本土化的例子：
代码语言：javascript复制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等
