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

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

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

https://proplot.readthedocs.io/en/latest/api/proplot.constructor.Proj.html#proj-table

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

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 方法)快速绘制坐标。

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:

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 函数,也可以达到相同的效果

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 的方法(略微有些复杂)。

该方法部分参照:

https://zhajiman.github.io/post/cartopy_lambert/
https://gist.github.com/ajdawson/dd536f786741e987ae4e

代码如下:

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。

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)

原文地址:https://cloud.tencent.com/developer/article/2135668

版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 dio@foxmail.com 举报,一经查实,本站将立刻删除。

相关推荐


学习编程是顺着互联网的发展潮流,是一件好事。新手如何学习编程?其实不难,不过在学习编程之前你得先了解你的目的是什么?这个很重要,因为目的决定你的发展方向、决定你的发展速度。
IT行业是什么工作做什么?IT行业的工作有:产品策划类、页面设计类、前端与移动、开发与测试、营销推广类、数据运营类、运营维护类、游戏相关类等,根据不同的分类下面有细分了不同的岗位。
女生学Java好就业吗?女生适合学Java编程吗?目前有不少女生学习Java开发,但要结合自身的情况,先了解自己适不适合去学习Java,不要盲目的选择不适合自己的Java培训班进行学习。只要肯下功夫钻研,多看、多想、多练
Can’t connect to local MySQL server through socket \'/var/lib/mysql/mysql.sock问题 1.进入mysql路径
oracle基本命令 一、登录操作 1.管理员登录 # 管理员登录 sqlplus / as sysdba 2.普通用户登录
一、背景 因为项目中需要通北京网络,所以需要连vpn,但是服务器有时候会断掉,所以写个shell脚本每五分钟去判断是否连接,于是就有下面的shell脚本。
BETWEEN 操作符选取介于两个值之间的数据范围内的值。这些值可以是数值、文本或者日期。
假如你已经使用过苹果开发者中心上架app,你肯定知道在苹果开发者中心的web界面,无法直接提交ipa文件,而是需要使用第三方工具,将ipa文件上传到构建版本,开...
下面的 SQL 语句指定了两个别名,一个是 name 列的别名,一个是 country 列的别名。**提示:**如果列名称包含空格,要求使用双引号或方括号:
在使用H5混合开发的app打包后,需要将ipa文件上传到appstore进行发布,就需要去苹果开发者中心进行发布。​
+----+--------------+---------------------------+-------+---------+
数组的声明并不是声明一个个单独的变量,比如 number0、number1、...、number99,而是声明一个数组变量,比如 numbers,然后使用 nu...
第一步:到appuploader官网下载辅助工具和iCloud驱动,使用前面创建的AppID登录。
如需删除表中的列,请使用下面的语法(请注意,某些数据库系统不允许这种在数据库表中删除列的方式):
前不久在制作win11pe,制作了一版,1.26GB,太大了,不满意,想再裁剪下,发现这次dism mount正常,commit或discard巨慢,以前都很快...
赛门铁克各个版本概览:https://knowledge.broadcom.com/external/article?legacyId=tech163829
实测Python 3.6.6用pip 21.3.1,再高就报错了,Python 3.10.7用pip 22.3.1是可以的
Broadcom Corporation (博通公司,股票代号AVGO)是全球领先的有线和无线通信半导体公司。其产品实现向家庭、 办公室和移动环境以及在这些环境...
发现个问题,server2016上安装了c4d这些版本,低版本的正常显示窗格,但红色圈出的高版本c4d打开后不显示窗格,
TAT:https://cloud.tencent.com/document/product/1340