Картопия: метка оси - временное решение

Я ищу обходной путь для добавления отметок и меток осей x и y на карту Cartopy в проекции Ламберта.

Решение, которое я придумал, представляет собой всего лишь приближение, которое даст худшие результаты для больших карт: оно включает в себя преобразование желаемых местоположений отметок в проекцию карты с использованием метода transform_points. Для этого я использую среднюю долготу (или широту) моей оси y (или оси x) вместе с желаемыми положениями отметок широты (или долготы) для вычисления координат проекции карты. См. Код ниже.

Таким образом, я предполагаю постоянные долготы по оси y (широты по оси x), что неверно и, следовательно, приводит к отклонениям. (Обратите внимание на разницу в прилагаемом результирующем рисунке: значение 46 ° в set_extent и результирующая позиция галочки).

Есть ли более точные решения? Любые подсказки, как я мог бы подойти к этой проблеме иначе?

Спасибо за любые идеи!

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

введите описание изображения здесь


person yngwaz    schedule 15.01.2015    source источник
comment
Интересно, что Cartopy пока поддерживает только две проекции.   -  person Chang    schedule 01.08.2019


Ответы (4)


Мой (довольно грубый) обходной путь подробно описан в этой записной книжке: http://nbviewer.ipython.org/gist/ajdawson/dd536f786741e987ae4e

Ноутбук требует картографии> = 0.12.

Все, что я сделал, это нашел пересечение соответствующей линии сетки с границей карты. Я предполагал, что граница карты всегда будет прямоугольной, и я могу пометить только нижнюю и левую стороны. Надеюсь, это может быть полезно для дальнейшего развития.

person ajdawson    schedule 07.04.2015
comment
Спасибо за такой подход! Линии сетки выглядят несколько искаженными, но функция установки галочки работает очень хорошо! - person yngwaz; 04.05.2015

Я сам не пробовал, но заметил в salem пакет документов появилась возможность обрабатывать линии сетки других проекций с помощью собственной утилиты построения графиков, которая не меняет matplotlib проекцию осей.

person ryanjdillon    schedule 03.09.2017

Маркировка линий сетки теперь поддерживается в любой проекции картографии, начиная с версии 0.18.0. https://twitter.com/QuLogic/status/1257148289838911488

person Ryan    schedule 19.06.2020

К сожалению, все еще с версией 0.18 разметка координатной сетки во всех проекциях у меня была невозможна. Мне пришлось изменить решение, предложенное в этом репозитории на github, чтобы оно работало для меня. Вот мой рабочий обходной путь:

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
person Vinod Kumar    schedule 16.04.2021