【ProPlot库】ProPlot3兰伯特投影-可添加刻度(三)

2022-10-09 09:53:52 浏览数 (1)

虽然 cartopy 下的 Plate Carrée 投影使用方便,但中纬度下使用 Lambert 投影能更好的呈现真实的地图。用一个正圆锥切于或割于球面,将地球面投影到圆锥面上,然后沿一母线展开成平面。下图是使用proplot绘制的最终效果:

在 proplot 中,可在以下链接找到相关投影名称表,其中兰伯特投影简称 'lcc'。

代码语言:javascript复制
https://proplot.readthedocs.io/en/latest/api/proplot.constructor.Proj.html#proj-table

下面直观感受一下两种投影的不同:

代码语言:javascript复制
Import proplot as plot
fig,axs = plot.subplots(ncols=2,width=5,height=3,projection=['cyl','lcc'])

fig.format(reso='hi',coast=True,metalinewidth=2,coastlinewidth=0.5)
axs[0].format(title='Equidistant Cylindrical (Plate Carrée)')
axs[1].format(title='Lambert Conformal Conic')

也可以使用 plot.Proj 来具体定义投影的参数和中央经线位置。同时可在format 里设置 labels=True (这是直接利用cartopy的 gridlines 方法)快速绘制坐标。

代码语言:javascript复制
import proplot as plot
proj = plot.Proj('lcc', lon_0=105)
fig,axs = plot.subplots(ncols=2,width=5,height=3,projection=['cyl',proj])

fig.format(reso='hi',coast=True,metalinewidth=2,coastlinewidth=0.5,
           lonlim=(80,130),latlim=(15,55),labels=True)

axs[0].format(title='Equidistant Cylindrical (Plate Carrée)')
axs[1].format(title='Lambert Conformal Conic')

也可以直接使用 cartopy 提供的 cartopy.mpl.geoaxes.GeoAxes.gridlines 方法添加 grid:

代码语言:javascript复制
import proplot as plot
proj = plot.Proj('lcc', lon_0=105)
fig,axs = plot.subplots(ncols=1,width=5,height=3,projection=proj)

axs.format(reso='hi',coast=True,metalinewidth=2,coastlinewidth=0.5,
           lonlim=(80,130),latlim=(15,55),title='Lambert Conformal Conic',titlesize=15,titlepad=10)

gl=axs[0].gridlines(
    xlocs=np.arange(-180, 180   1, 10), ylocs=np.arange(-90, 90   1, 10),
    draw_labels=True, x_inline=False, y_inline=False,
    linewidth=0.5, linestyle='--', color='gray')

gl.top_labels = False
gl.right_labels = False
gl.rotate_labels =0

如果只想用 proplot 的 format 函数,也可以达到相同的效果

代码语言:javascript复制
import proplot as plot
import numpy as np
proj = plot.Proj('lcc', lon_0=105)
fig,axs = plot.subplots(ncols=1,width=5,height=3,projection=proj)

axs.format(reso='hi',coast=True,metalinewidth=2,coastlinewidth=0.5,
           lonlim=(80,130),latlim=(15,55),title='Lambert Conformal Conic',titlesize=15,titlepad=10,grid=True,
           latlocator=np.arange(-180, 180   1, 10),labels=True,lonlabels=True,gridlabelsize=10,gridlabelpad=10,
           gridlinestyle='--', gridcolor='grey',gridlinewidth=0.5)
#https://proplot.readthedocs.io/en/latest/api/proplot.axes.GeoAxes.html?highlight=geo#proplot.axes.GeoAxes

在上述例子中,gridlabelsize=10 , gridlabelpad=10 分别控制经纬度标签的大小和离axis的距离;gridlinestyle='--', gridcolor='grey' , gridlinewidth=0.5 来改变 gridline 的风格。追求快速出图的朋友可以采用这种方案绘制底图。

添加自定义刻度

虽然 cartopy 提供了 gridlines 的方法,但经度的标签出现在了 y 轴上(80°E),并且和matplotlib 的 matplotlib.axes.Axes.set_xticks 等并不相融,难以给出 tickmarker ;同时,ncl 可以为顶部添加 ticks label。网上的解决方案并不多,自己尝试了很多办法。下面提供一种可以给 cartopy 兰伯特投影添加 major、minor ticks 的方法(略微有些复杂)。

该方法部分参照:

代码语言:javascript复制
https://zhajiman.github.io/post/cartopy_lambert/
https://gist.github.com/ajdawson/dd536f786741e987ae4e

代码如下:

代码语言:javascript复制
def add_alt(self,**kwargs):


    locator = self._make_inset_locator([0, 0, 1, 1], self.transAxes)
    ax = plot.CartesianAxes(self.figure, locator(self, None).bounds, **kwargs)
    ax.set_axes_locator(locator)

    cx=self.add_child_axes(ax)  # to facilitate tight layout
    cx.grid(False)
    limsx=self.get_xlim()
    limsy=self.get_ylim()

    cx.set_xlim(limsx)
    cx.set_ylim(limsy)
    return cx

def find_side(ls, side):
    """
    Given a shapely LineString which is assumed to be rectangular, return the
    line corresponding to a given side of the rectangle.

    """
    minx, miny, maxx, maxy = ls.bounds
    points = {'left': [(minx, miny), (minx, maxy)],
              'right': [(maxx, miny), (maxx, maxy)],
              'bottom': [(minx, miny), (maxx, miny)],
              'top': [(minx, maxy), (maxx, maxy)],}
    return sgeom.LineString(points[side])


def lambert_xticks(ax,ax2, ticks,po='bottom',formatter = LongitudeFormatter()):
"""Draw ticks on the bottom x-axis of a Lambert Conformal projection."""
""" ax用于计算交点,ax2用于呈现ticks """
    te = lambda xy: xy[0]
    
    lc = lambda t, n, b: np.vstack((np.zeros(n)   t, np.linspace(b[2] 0.1*(b[2]-b[3]), b[3]-0.1*(b[2]-b[3]), n))).T
    xticks, xticklabels = _lambert_ticks(ax, ticks, po, lc, te)
    #ax.xaxis.tick_bottom()
    ax2.set_xticks(xticks)
    #print(xticks)
    #formatter = LongitudeFormatter()
    ax2.grid(False)
    ax2.set_xticklabels([formatter(xtick) for xtick in xticklabels])


def lambert_yticks(ax,ax2, ticks,po='left',formatter=LatitudeFormatter()):
    """Draw ricks on the left y-axis of a Lamber Conformal projection."""
    te = lambda xy: xy[1]
    lc = lambda t, n, b: np.vstack((np.linspace(b[0], b[1], n), np.zeros(n)   t)).T
    yticks, yticklabels = _lambert_ticks(ax, ticks, po, lc, te)
    #ax.yaxis.tick_left()
    ax2.set_yticks(yticks)
    #formatter = LatitudeFormatter()
    ax2.set_yticklabels([formatter(ytick) for ytick in yticklabels])
    ax2.grid(False)

def _lambert_ticks(ax, ticks, tick_location, line_constructor, tick_extractor):
    """Get the tick locations and labels for an axis of a Lambert Conformal projection."""
    outline_patch = sgeom.LineString(ax.spines['geo'].get_path().vertices.tolist())
    axis = find_side(outline_patch, tick_location)
    n_steps = 30
    extent = ax.get_extent(ccrs.PlateCarree())
    _ticks = []
    
    for t in ticks:


        xy = line_constructor(t, n_steps, extent)

        proj_xyz = ax.projection.transform_points(ccrs.Geodetic(), xy[:, 0], xy[:, 1])

        xyt = proj_xyz[..., :2]
        ls = sgeom.LineString(xyt.tolist())

        locs = axis.intersection(ls)
        if not locs:
            tick = [None]
        else:
            tick = tick_extractor(locs.xy)
        #print(tick[0])
        _ticks.append(tick[0])
    # Remove ticks that aren't visible:    
    ticklabels = copy(ticks)
    while True:
        try:
            index = _ticks.index(None)
        except ValueError:
            break
        _ticks.pop(index)
        ticklabels.pop(index)

    return _ticks, ticklabels




def lambert_minorxticks(ax,ax2, ticks,po='bottom',formatter = LongitudeFormatter()):
    """Draw ticks on the bottom x-axis of a Lambert Conformal projection."""
    te = lambda xy: xy[0]
    lc = lambda t, n, b: np.vstack((np.zeros(n)   t, np.linspace(b[2] 0.1*(b[2]-b[3]), b[3]-0.1*(b[2]-b[3]), n))).T 
    
    xticks, xticklabels = _lambert_ticks(ax, ticks, po, lc, te)

    ax2.grid(False)
    ax2.xaxis.set_minor_locator(mticker.FixedLocator(xticks))

    
def lambert_minoryticks(ax,ax2, ticks,po='left',formatter=LatitudeFormatter()):
    """Draw ricks on the left y-axis of a Lamber Conformal projection."""
    te = lambda xy: xy[1]
    lc = lambda t, n, b: np.vstack((np.linspace(b[0], b[1], n), np.zeros(n)   t)).T
    yticks, yticklabels = _lambert_ticks(ax, ticks, po, lc, te)
    ax2.yaxis.set_minor_locator(mticker.FixedLocator(yticks))
    ax2.grid(False)


#以下是画图的代码:

proj = plot.Proj('lcc', lon_0=105)
proj2 = plot.Proj('lcc', lon_0=115)

fig, axs = plot.subplots(ncols=2,refwidth=5,refheight=5.5,projection=proj,outerpad=4,share=False)


xticks = np.arange(65, 185 1, 20).tolist()
yticks = np.arange(0,  50 1,  10).tolist()
xmiticks = np.arange(65, 185 1, 4 ).tolist()
ymiticks = np.arange(0,  50 1,  2.5).tolist()
xticks2 = np.arange(65, 185, 10)
yticks2 = np.arange(0, 50, 10)

axs.format(grid=True,gridalpha=0.5,gridcolor='grey',lonlim=(80,130),latlim=(15,55),reso='hi',labels=False,coast=True,metalinewidth=2,
         lonlocator=xticks,latlocator=yticks
         ,gridlinestyle='--',gridlinewidth = 2)

fig.canvas.draw() #在proplot中,如果没有插入shape地图,则要加上这句,原因未知(matplotlib不用)

ax=axs[0]

ax.format(title='without top and right ticks',titlesize=20,titlepad=10)


lambert_xticks(ax,ax, xticks)
lambert_yticks(ax,ax,yticks)

lambert_minorxticks(ax,ax,xmiticks)
lambert_minoryticks(ax,ax,ymiticks)

ax.tick_params(which='major', width=1.5, length=8,bottom=True,left=True,right=False,top=False , )
ax.tick_params(labelsize=15, pad=5)
ax.tick_params(which='minor', width=0.8,length=4, color='k',bottom=True,left=True,right=False,top=False ,)

#添加南海子图
ins = ax.inset_axes([0.825,0.001,0.2,0.2],projection =proj2) #proplot中inset_axes可以自己选择子图的projection,非常便捷


ins.format(grid=True,coast=True,lonlim=(105,125),latlim=(0,25),metalinewidth=1.5,reso='hi',coastlinewidth=0.2,
          lonlocator=xticks2, latlocator=yticks2,gridlinestyle='--',gridlinewidth = 1)

ax=axs[1]


lambert_xticks(ax,ax, xticks)
lambert_yticks(ax,ax,yticks)

lambert_minorxticks(ax,ax,xmiticks)
lambert_minoryticks(ax,ax,ymiticks)

#添加secondary axis (CartesianAxes projection) 在matplotlib中可直接用ax.secondary_xaxis 和 ax.secondary_yaxis
#也可以直接ax2 = ax.inset_axes([0.,0.00,1,1],projection =proj)
ax2=add_alt(ax)

lambert_xticks(ax,ax2, xticks,'top') #由ax判断交点,在ax2上绘制
lambert_yticks(ax,ax2, yticks,'right')

lambert_minorxticks(ax,ax2,xmiticks,'top')
lambert_minoryticks(ax,ax2,ymiticks,'right')


ax.tick_params(which='major', width=1.5, length=8,bottom=True,left=True,right=False,top=False , )
ax.tick_params(labelsize=15, pad=5)
ax.tick_params(which='minor', width=0.8,length=4, color='k',bottom=True,left=True,right=False,top=False ,)



ax2.tick_params(which='major', width=1.5, length=8,bottom=False,left=False,right=True,top=True,color='k' ,labelleft=False,labelbottom=False
                ,labelright=True, labeltop=True)
ax2.tick_params(labelsize=15, pad=5)
ax2.tick_params(which='minor',width=0.8, length=4,bottom=False,left=False,right=True,top=True,color='k')


#添加南海子图
ins = ax.inset_axes([0.825,0.001,0.2,0.2],projection =proj2) 

ins.format(grid=True,coast=True,lonlim=(105,125),latlim=(0,25),metalinewidth=1.5,reso='hi',coastlinewidth=0.2,
          lonlocator=xticks2, latlocator=yticks2,gridlinestyle='--',gridlinewidth = 1)

利用 Shapley 判断网格线与一个轴线相交的位置,进而找到刻度的位置。传入set_xticksset_yticks 中,再利用ax.set_xticklabels 添加刻度。在proplot中,目前我只找到了用ax.tick_params添加geoaxes 的tickmarker 的例子,在format 里给地图设置刻度标签的方法可能还未实现。

目前网上较多的案例是添加左侧和底部的majortick,而由于 top 和right 轴的坐标范围发生变化,于是模仿Secondary Axis 方法,叠加新的“图层”来为顶部和右侧添加处于不同位置的 ticks,add_alt(self,**kwargs) 是专门给 proplot 改编的,如果想用matplotlib 绘制,便可以直接采用axes.Axes.secondary_xaxis 和axes.Axes.secondary_yaxis ,为geoAxes 添加一个笛卡尔坐标系的轴,再使用 lambert_xticks 函数添加标签。

在实际画图中也可以只选择采用左图的方案。另外,proplot 的inset_axes 方法可以添加任意投影(proj参数)的小子图,比matplotlib 的inset_axes 方便很多。

再添加如下代码,完成等值线和 colorbar 的绘制;当然,也可也不画顶部和右侧的 tick,或是 ticklabel 和 tickmarker。

代码语言:javascript复制
m=ax.contourf(T[0],cmap=cmaps.BlueWhiteOrangeRed,lw=0.2,levels=12,vmax=316,vmin=264)

cb=ax.colorbar(m,loc='b',format='%.0f',pad=0.0, shrink=1 ,drawedges=True,width=0.2,linewidth=1,
              ticklabelsize=12,label='')

m2=ins.contourf(T[0],cmap=cmaps.BlueWhiteOrangeRed,lw=0.2,levels=12,vmax=316,vmin=260)

ax.format(grid=False,ltitle='2m Temperature',titlesize=15,titlepad=15,rtitle='(K)',suptitle='Your Second Plot' ,suptitlesize=20,suptitlepad=15)

0 人点赞