Cartopy:轴标签 - 解决方法

yng*_*waz 5 python matplotlib cartopy

我正在寻找一种解决方法,在Lambert投影中将x和y轴刻度和标签添加到Cartopy贴图中.

我提出的解决方案只是一个近似值,对于较大的地图会产生更糟糕的结果:它涉及使用transform_points方法将所需的刻度位置转换为地图投影.为此,我使用y轴(或x轴)的中位数经度(或纬度)以及所需的纬度(或经度)刻度位置来计算地图投影坐标.见下面的代码.

因此,我假设沿y轴的恒定经度(沿x轴的纬度),这是不正确的,因此导致偏差.(注意附加结果图中的差异:在set_extent中设置46°并产生刻度位置).

有没有更精确的解决方案?有什么提示我怎么能解决这个问题呢?

感谢任何想法!

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import numpy as np

def main():
    #my desired Lambert projection:
    myproj = ccrs.LambertConformal(central_longitude=13.3333, central_latitude=47.5,
                                   false_easting=400000, false_northing=400000,
                                   secant_latitudes=(46, 49))

    arat = 1.1 #just some factor for the aspect ratio
    fig_len = 12
    fig_hig = fig_len/arat
    fig = plt.figure(figsize=(fig_len,fig_hig), frameon=True)
    ax = fig.add_axes([0.08,0.05,0.8,0.94], projection = myproj)

    ax.set_extent([10,16,46,49])
    #This is what is not (yet) working in Cartopy due to Lambert projection:
    #ax.gridlines(draw_labels=True) #TypeError: Cannot label gridlines on a LambertConformal plot.  Only PlateCarree and Mercator plots are currently supported.
    x_lons = [12,13,14] #want these longitudes as tick positions
    y_lats = [46, 47, 48, 49] #want these latitudes as tick positions
    tick_fs = 16
    #my workaround functions:
    cartopy_xlabel(ax,x_lons,myproj,tick_fs)
    cartopy_ylabel(ax,y_lats,myproj,tick_fs)

    plt.show()
    plt.close()

def cartopy_xlabel(ax,x_lons,myproj,tick_fs):    
    #transform the corner points of my map to lat/lon
    xy_bounds = ax.get_extent()
    ll_lonlat = ccrs.Geodetic().transform_point(xy_bounds[0],xy_bounds[2], myproj)
    lr_lonlat = ccrs.Geodetic().transform_point(xy_bounds[1],xy_bounds[2], myproj)
    #take the median value as my fixed latitude for the x-axis
    l_lat_median = np.median([ll_lonlat[1],lr_lonlat[1]]) #use this lat for transform on lower x-axis
    x_lats_helper = np.ones_like(x_lons)*l_lat_median

    x_lons = np.asarray(x_lons)
    x_lats_helper = np.asarray(x_lats_helper)
    x_lons_xy = myproj.transform_points(ccrs.Geodetic(), x_lons,x_lats_helper)
    x_lons_xy = list(x_lons_xy[:,0]) #only lon pos in xy are of interest     
    x_lons = list(x_lons)

    x_lons_labels =[]
    for j in xrange(len(x_lons)):
        if x_lons[j]>0:
            ew=r'$^\circ$E'
        else:
            ew=r'$^\circ$W'
        x_lons_labels.append(str(x_lons[j])+ew)
    ax.set_xticks(x_lons_xy)
    ax.set_xticklabels(x_lons_labels,fontsize=tick_fs)

def cartopy_ylabel(ax,y_lats,myproj,tick_fs):        
    xy_bounds = ax.get_extent()
    ll_lonlat = ccrs.Geodetic().transform_point(xy_bounds[0],xy_bounds[2], myproj)
    ul_lonlat = ccrs.Geodetic().transform_point(xy_bounds[0],xy_bounds[3], myproj)
    l_lon_median = np.median([ll_lonlat[0],ul_lonlat[0]]) #use this lon for transform on left y-axis
    y_lons_helper = np.ones_like(y_lats)*l_lon_median

    y_lats = np.asarray(y_lats)    
    y_lats_xy = myproj.transform_points(ccrs.Geodetic(), y_lons_helper, y_lats)
    y_lats_xy = list(y_lats_xy[:,1]) #only lat pos in xy are of interest 

    y_lats = list(y_lats)

    y_lats_labels =[]
    for j in xrange(len(y_lats)):
        if y_lats[j]>0:
            ew=r'$^\circ$N'
        else:
            ew=r'$^\circ$S'
        y_lats_labels.append(str(y_lats[j])+ew)
    ax.set_yticks(y_lats_xy)
    ax.set_yticklabels(y_lats_labels,fontsize=tick_fs)

if __name__ == '__main__': main()
Run Code Online (Sandbox Code Playgroud)

在此输入图像描述

ajd*_*son 5

我的(相当粗略的)解决方法在此笔记本中进行了详细说明:http : //nbviewer.ipython.org/gist/ajdawson/dd536f786741e987ae4e

笔记本电脑需要的Cartopy> = 0.12。

我所要做的就是找到适当的网格线与地图边界的交点。我假设地图边界将始终为矩形,并且只能标记底部和左侧。希望这可能会有所帮助。