Python气象绘图教程(十三)

您所在的位置:网站首页 枪花乐队经典专辑 Python气象绘图教程(十三)

Python气象绘图教程(十三)

2023-11-01 03:50| 来源: 网络整理| 查看: 265

本节提要:关于子图的一些问题、使用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等



【本文地址】


今日新闻


推荐新闻


CopyRight 2018-2019 办公设备维修网 版权所有 豫ICP备15022753号-3