Cartopy: axis label - workaround

Question

I am looking for a workaround to add x and y axis ticks and labels to a Cartopy map in Lambert projection.

The solution I have come up with is just an approximation which will yield worse results for larger maps: It involves transforming desired tick locations to map projection using the transform_points method. For this I use the median longitude (or latitude) of my y axis (or x axis) together with the desired latitudes (or longitudes) tick positions to compute map projection coordinates. See code below.

Thus, I am assuming constant longitudes along the y-axis (latitudes along the x-axis), which is not correct and hence leads to deviations. (Note the difference in the attached resulting figure: 46° set in set_extent and resulting tick position).

Are there any solutions out there which are more accurate? Any hints how I could approach this problem otherwise?

Thank's for any ideas!

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()

Answer 1

My (quite crude) work-around to this is detailed in this notebook: http://nbviewer.ipython.org/gist/ajdawson/dd536f786741e987ae4e

The notebook requires cartopy >= 0.12.

All I've done is find the intersection of the appropriate gridline with the map boundary. I've assumed the map boundary will always be rectangular, and I can only label the bottom and left sides. Hopefully this might be useful as something to build on.