问题

当绘图时,往往并不需要绘制整块区域,而是想聚焦于局部地区,此时我们需要绘制扇形图。
在cartopy中,只提供普通正方形的框架,如果我们需要其他,边界,需要自己去绘制,最常见的是使用set_boundary绘制边界,set_extend规定边界,如:

latmin = 10
latmax = 60
lonmin = 70
lonmax = 140
lats = np.linspace(latmax, latmin, latmax - latmin + 1)
lons = np.linspace(lonmin, lonmax, lonmax - lonmin + 1)
vertices = [(lon, latmin) for lon in range(lonmin, lonmax + 1, 1)] + \
    [(lon, latmax) for lon in range(lonmax, lonmin - 1, -1)]
boundary = mpath.Path(vertices)
ax.set_boundary(boundary, transform=ccrs.PlateCarree())
ax.set_extent(extent)

需要注意的是,在使用set_boundary裁剪边界后,需要使用set_extend规定边界,否则会出现如下情况:
在这里插入图片描述

但是,这种方法有着自己的bug与缺陷:比如

  1. 需要自己添加网格和经纬度标签
  2. set_boundary和 set_extend一起使用时,绘制到边界消失,如:cartopy二次曲线外观保持
  3. 添加坐标时可能会有A LinearRing must have at least 3 coordinate tuples的报错。

解决方式

实际上,这个问题原因还是由于投影转换的问题,在set_extend时,绘制的上下边界仍然是方形、未被正确投影的边界,与我们的set_boundary存在冲突,最根本的原因还是在于cartopy对于投影计算的一些缺陷。
这类问题的解决在:GeoAxes.set_extent on non-cylindrical projections得到了充分的讨论与解决,主要是使用了这几行代码:

lon1, lon2, lat1, lat2 = [-20, 20, 50, 80]
rect = mpath.Path([[lon1, lat1], [lon2, lat1],
    [lon2, lat2], [lon1, lat2], [lon1, lat1]]).interpolated(50)
proj_to_data = ccrs.PlateCarree()._as_mpl_transform(ax) - ax.transData
rect_in_target = proj_to_data.transform_path(rect)

将我们绘制的边界进行了转换。
我自己绘制到图形示例如下:

import matplotlib.pyplot               as plt
import matplotlib.path                 as mpath
import cartopy.crs                     as ccrs
import cartopy.mpl.ticker              as ctk
import cartopy.feature as cfeature
import cmaps


cmap = cmaps.BlueDarkRed18
leftlon, rightlon, lowerlat, upperlat = [0,120,65,85]

rect = mpath.Path([[leftlon, lowerlat], [rightlon, lowerlat],
    [rightlon, upperlat], [leftlon, upperlat], [leftlon, lowerlat]]).interpolated(50)

proj=ccrs.NearsidePerspective(central_longitude=(leftlon+rightlon)*0.5,
    central_latitude=(lowerlat+upperlat)*0.5)
fig= plt.figure(figsize=(8,8),dpi=300)#设置画布大小
ax= fig.add_axes([0.2,0.3,0.5,0.5],projection =proj)
#ax.set_extent(extent, ccrs.PlateCarree())#这样截出来是方形的
ax.add_feature(cfeature.COASTLINE.with_scale('110m'))

proj_to_data = ccrs.PlateCarree()._as_mpl_transform(ax) - ax.transData
rect_in_target = proj_to_data.transform_path(rect)
ax.set_boundary(rect_in_target)
ax.set_xlim(rect_in_target.vertices[:,0].min(), rect_in_target.vertices[:,0].max())
ax.set_ylim(rect_in_target.vertices[:,1].min(), rect_in_target.vertices[:,1].max())

gl=ax.gridlines(draw_labels=True, x_inline=False, y_inline=False, linestyle='dashed')
gl.top_labels=False
gl.right_labels=False
gl.rotate_labels=False
gl.xlocator=ctk.LongitudeLocator(4)
gl.ylocator=ctk.LatitudeLocator(6)
gl.xformatter=ctk.LongitudeFormatter(zero_direction_label=True)
gl.yformatter=ctk.LatitudeFormatter()
ax.set_title('AMJ-pc1 & SON-SIC',loc='center',fontsize =12)

c1 = ax.contourf(lon1,lat1, s, zorder=0,levels =np.arange(-0.09,0.12,0.03),
                 extend = 'both',cmap=cmap,transform=ccrs.PlateCarree())
c1b =ax.contourf(lon1,lat1, p,[0,0.1 ,1], zorder=1,hatches=['...', None],
                 colors="none",transform=ccrs.PlateCarree())

position=fig.add_axes([0.2, 0.2, 0.5, 0.02])
fig.colorbar(c1,cax=position,orientation='horizontal')#,format='%.2f',)
plt.show()

在这里插入图片描述
另外我们需要指出的是:**该方法不适用于极地投影,即NorthPolarStereo,由于NorthPolarStereo本身投影特性只需一个参数,本身并不适合。
**

极地局部绘制(不推荐)

我们绘制极地投影时,同样也是使用set_boundary绘制圆形边界,那么当我们想要绘制扇形时,可以通过只绘制一部分的圆形,通过调整绘制圆的参数来局部绘制。
详细例子可见:极地局部图绘制
在我进行实验后发现它存在一些缺陷:
1、难以准确截出自己需要的区域
2、图形位置不好确定。
3、同样不好添加网格标签。
该方法本质是:绘制整个圆形边界,而只显示部分区域,在绘图时会给人不过完整的感觉。
直接绘制如下图,需要自行裁剪或调整绘制:
在这里插入图片描述
请根据喜好自行选择,

Logo

开放原子开发者工作坊旨在鼓励更多人参与开源活动,与志同道合的开发者们相互交流开发经验、分享开发心得、获取前沿技术趋势。工作坊有多种形式的开发者活动,如meetup、训练营等,主打技术交流,干货满满,真诚地邀请各位开发者共同参与!

更多推荐