6

我正在寻找一种解决方法,将 x 和 y 轴刻度和标签添加到 Lambert 投影中的 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()

在此处输入图像描述

4

4 回答 4

5

本笔记本详细介绍了我的(相当粗略的)解决方法:http: //nbviewer.ipython.org/gist/ajdawson/dd536f786741e987ae4e

笔记本需要 cartopy >= 0.12。

我所做的就是找到适当的网格线与地图边界的交点。我假设地图边界总是矩形的,我只能标记底部和左侧。希望这可能是有用的东西。

于 2015-04-07T09:55:07.527 回答
2

我自己没有尝试过,但我注意到在salem包文档中可以使用他们自己开发的绘图实用程序处理其他投影的网格线,这不会改变matplotlib“轴”投影。

于 2017-09-03T18:32:00.567 回答
0

自 cartopy v0.18.0 起,任何 cartopy 投影现在都支持标注网格线。https://twitter.com/QuLogic/status/1257148289838911488

于 2020-06-19T16:55:06.127 回答
0

不幸的是,仍然使用 0.18 版,我无法在所有投影中标记轴网格。我不得不修改这个 github repo上提出的解决方案来为我工作。这是我的解决方法:

def gridlines_with_labels(ax, top=True, bottom=True, left=True,
                      right=True, **kwargs):
"""
Like :meth:`cartopy.mpl.geoaxes.GeoAxes.gridlines`, but will draw
gridline labels for arbitrary projections.
Parameters
----------
ax : :class:`cartopy.mpl.geoaxes.GeoAxes`
    The :class:`GeoAxes` object to which to add the gridlines.
top, bottom, left, right : bool, optional
    Whether or not to add gridline labels at the corresponding side
    of the plot (default: all True).
kwargs : dict, optional
    Extra keyword arguments to be passed to :meth:`ax.gridlines`.
Returns
-------
:class:`cartopy.mpl.gridliner.Gridliner`
    The :class:`Gridliner` object resulting from ``ax.gridlines()``.
Example
-------
>>> import matplotlib.pyplot as plt
>>> import cartopy.crs as ccrs
>>> plt.figure(figsize=(10, 10))
>>> ax = plt.axes(projection=ccrs.Orthographic(-5, 53))
>>> ax.set_extent([-10.0, 0.0, 50.0, 56.0], crs=ccrs.PlateCarree())
>>> ax.coastlines('10m')
>>> gridlines_with_labels(ax)
>>> plt.show()
"""

# Add gridlines
gridliner = ax.gridlines(**kwargs)

ax.tick_params(length=0)

# Get projected extent
xmin, xmax, ymin, ymax = ax.get_extent()

# Determine tick positions
sides = {}
N = 500
if bottom:
    sides['bottom'] = np.stack([np.linspace(xmin, xmax, N),
                                np.ones(N) * ymin])
if top:
    sides['top'] = np.stack([np.linspace(xmin, xmax, N),
                            np.ones(N) * ymax])
if left:
    sides['left'] = np.stack([np.ones(N) * xmin,
                              np.linspace(ymin, ymax, N)])
if right:
    sides['right'] = np.stack([np.ones(N) * xmax,
                               np.linspace(ymin, ymax, N)])

# Get latitude and longitude coordinates of axes boundary at each side
# in discrete steps
gridline_coords = {}
for side, values in sides.items():
    gridline_coords[side] = ccrs.PlateCarree().transform_points(
        ax.projection, values[0], values[1])

lon_lim, lat_lim = gridliner._axes_domain()
ticklocs = {
    'x': gridliner.xlocator.tick_values(lon_lim[0], lon_lim[1]),
    'y': gridliner.ylocator.tick_values(lat_lim[0], lat_lim[1])
}

# Compute the positions on the outer boundary where
coords = {}
for name, g in gridline_coords.items():
    if name in ('bottom', 'top'):
        compare, axis = 'x', 0
    else:
        compare, axis = 'y', 1
    coords[name] = np.array([
        sides[name][:, np.argmin(np.abs(
            gridline_coords[name][:, axis] - c))]
        for c in ticklocs[compare]
    ])

# Create overlay axes for top and right tick labels
ax_topright = ax.figure.add_axes(ax.get_position(), frameon=False)
ax_topright.tick_params(
    left=False, labelleft=False,
    right=right, labelright=right,
    bottom=False, labelbottom=False,
    top=top, labeltop=top,
    length=0
)
ax_topright.set_xlim(ax.get_xlim())
ax_topright.set_ylim(ax.get_ylim())

for side, tick_coords in coords.items():
    if side in ('bottom', 'top'):
        axis, idx = 'x', 0
    else:
        axis, idx = 'y', 1

    _ax = ax if side in ('bottom', 'left') else ax_topright

    ticks = tick_coords[:, idx]

    valid = np.logical_and(
        ticklocs[axis] >= gridline_coords[side][0, idx],
        ticklocs[axis] <= gridline_coords[side][-1, idx])

    if side in ('bottom', 'top'):
        _ax.set_xticks(ticks[valid])
        _ax.set_xticklabels([LONGITUDE_FORMATTER.format_data(t)
                             for t in ticklocs[axis][valid]])
    else:
        _ax.set_yticks(ticks[valid])
        _ax.set_yticklabels([LATITUDE_FORMATTER.format_data(t)
                             for t in np.asarray(ticklocs[axis])[valid]])

return gridliner
于 2021-04-16T20:25:35.843 回答