Add support for all projections to lon/lat labelling

[stefraynaud]

Rationale

Longitude / latitude labelling is restricted to Mercator and PlateCarree projections.
This restriction can be removed by computing the position of labels as the intersection of grid lines with the map and plot boundaries. This makes it compatible with all projections.

In the current version, only labels that are the border of the plot are displayed, not those at the border of the map. This limitation can be overcome by better checking if a label is plotted over another object, like a label or the map.

The code below produces the attached figures.

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


plain_crs = ccrs.PlateCarree()
mercator_crs = ccrs.Mercator()
rotated_crs = ccrs.RotatedPole(pole_longitude=120.0, pole_latitude=45.0)
igh_crs = ccrs.InterruptedGoodeHomolosine()
polar_crs = ccrs.SouthPolarStereo()
uk_extent = [-10, 3, 49, 59]
north_atlantic_extent = [-90, 0, -10, 70]
polar_extent = [-180, 180, -90, -60]

to_test = [(plain_crs, uk_extent),
           (mercator_crs, uk_extent),
           (rotated_crs, uk_extent),
           (rotated_crs, north_atlantic_extent),
           (igh_crs, north_atlantic_extent),
           (polar_crs, polar_extent)]

for projection, extent in to_test:

    fig = plt.figure(figsize=(6, 6))
    ax = fig.add_subplot(1, 1, 1, projection=projection)
    ax.set_extent(extent, crs=plain_crs)
    ax.coastlines(resolution='50m', color='k')
    ax.gridlines(color='lightgrey', linestyle='-', draw_labels=True)
    plt.savefig(__file__[:-2] + projection.__class__.__name__.lower() + '.png')

plt.show()

Implications

Longitude and latitude labels can be plotted for all projections.
Another implication is that x/y labels are no longer associated with respectively longitudes/latitudes. This must be clarified in the internal API at least.

Note : This PR must be merged with PR #1089

[QuLogic]

Cool! I always thought something along these lines would be the way to go, but never got around to attempting it.

[kaedonkers]

This is a really great feature that's been needed for a while, thank you @stefraynaud!
I do not have enough experience with this part of the cartopy code to provide a full code review (@pelson, could you help with this?).
I ran the code you provided above to reproduce the plots and found:

1. I was not able to create some plots with the expected labels e.g. the Mercator plot of UK had no labels. Is this because this PR is waiting on Add inline grid meridian/parallel labeling #1089?
2. I got a verbose output of all the labels which were produced. Is this intended?

[stefraynaud]

I'm about to amend the PR with commits in which the distinction x/y and parralels/meridians is made. So it's gonna break a bit the API.

These commits will also include a commit to ticker to add support for minutes and seconds to the formatters, and to add longitude and latitude locators that also work at very small scale.

The changes will be added to this PR to make it cleaner and because they are used by the gridliner.

[stefraynaud]

Her is the kind of tick labels the last commit allows.

Degrees, minutes, and seconds can be auto hidden when not required.

[stefraynaud]

I'm working on a more universal way to plot labels, not only on plot axes limits, but also everywhere grid lines intersects the map boundary, whatever its shape.
I also added a clip keyword to the GeoAxes.set_boundary(), so it will be possible to easily make plots like this one : https://oceancurrents.rsmas.miami.edu/atlantic/img_mgsva/north-atlantic-drift-YYY.gif

So I'll amend this PR in the future.

[corinnebosley]

Hi @stefraynaud, I see that you're working on a new PR to achieve this goal a little bit better, but I think this PR along with #1089 are incredibly valuable.

Can you please let me know what's going on with this one now? I am happy to review it and to help progress it (because I really want this in cartopy, finally, and you've done a lot of work which I don't want to lose), but if it's dead and there is a replacement in progress then I will wait for that instead.

Also, I noticed that this PR must be merged with #1089; is that because of something functional or something to do with the docs? What I want to know is; if I borrow this branch to test it locally, will it work?

[stefraynaud]

Hi @corinnebosley

I'm still working actively on this PR and about to finish a version.
As you can seen on the figure below, I added a clipping fonctionality, and for this reason I generalized the way labels are plotted so that they can be plotted in most situations : the position and rotation of labels is determined by the intersection of grid lines with map boundary.

However, an important aspect of the new algorithm is that the visibility of all labels must updated during draw events to make sure labels are not overlapping themselves and not overlaping the map patch, contrary to what happened when saving the above figure for northern labels.
It works except when saving a figure, because the background path of each text label must be computed once the figure has its final shape and before the Axes.draw if finished. I'm currently having a look at Axes.apply_aspect

I still didn't tried to integrate #1089 because still trying to solve the above issue, but it will be very valuable to have the inline labels capability added too.

[stefraynaud]

The following code produces the following figures:

from __future__ import (absolute_import, division, print_function)

import matplotlib
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature

matplotlib.rc('font', size=6)

pcar = ccrs.PlateCarree()

# Global projections

ncol = 4
sy = 3.5
sx = 3.5
global_projs = [pcar,
              ccrs.Mercator(),
              ccrs.AlbersEqualArea(),
              ccrs.AzimuthalEquidistant(central_longitude=90),
#              ccrs.EquidistantConic(),
              ccrs.LambertConformal(),
              ccrs.LambertCylindrical(),
              ccrs.Miller(),
              ccrs.Mollweide(),
              ccrs.Orthographic(),
              ccrs.Robinson(),
              ccrs.Sinusoidal(),
              ccrs.Stereographic(),
              ccrs.TransverseMercator(),
              ccrs.InterruptedGoodeHomolosine(),
              ccrs.OSGB(),
              ccrs.EuroPP(),
              ccrs.Geostationary(),
              ccrs.NearsidePerspective(central_latitude=50.72,
                                       central_longitude=-3.53,
                                       satellite_height=10000000.0),
#              ccrs.EckertI(),
#              ccrs.EckertIII(),
#              ccrs.EqualEarth(),
              ccrs.Gnomonic(),
              ccrs.LambertAzimuthalEqualArea(),
              ccrs.NorthPolarStereo(),
              ccrs.OSNI(),
              ]

nrow = (len(global_projs)-1)//ncol + 1
fig = plt.figure(figsize=(sx*ncol, nrow * sy))
print('GLOBAL')
for ip, proj in enumerate(global_projs):
    print(' '+proj.__class__.__name__)
    ax = plt.subplot(*(nrow, ncol, ip+1), projection=proj)
    ax.gridlines(pcar, draw_labels=True)
    for feature in [cfeature.OCEAN, cfeature.LAND, cfeature.COASTLINE]:
        ax.add_feature(feature)
    ax.set_title(proj.__class__.__name__ + ' global')
fig.savefig(__file__[:-2]+'global.png')  # to make tight_layout working
fig.tight_layout(h_pad=1.5, w_pad=1.5, pad=0.1)
fig.savefig(__file__[:-2]+'global.png')


# Special cases

ncol = 2
sy = 3
sx = 4
uk_extent = [-10, 3, 49, 59]
north_atlantic_extent = [-90, 0, -10, 70]
special_projs = [(pcar, uk_extent, False),
                 (ccrs.RotatedPole(pole_longitude=120.0, pole_latitude=45.0),
                  north_atlantic_extent, False),
                 (ccrs.RotatedPole(pole_longitude=120.0, pole_latitude=45.0),
                  north_atlantic_extent, True)]

nrow = (len(special_projs)-1)//ncol + 1
fig = plt.figure(figsize=(sx*ncol, nrow * sy))
for ip, (proj, extents, clip) in enumerate(special_projs):
    ax = plt.subplot(*(nrow, ncol, ip+1), projection=proj)
    ax.set_extent(extents, crs=pcar, clip=clip)
    ax.gridlines(pcar, draw_labels=True)
    for feature in [cfeature.OCEAN, cfeature.LAND, cfeature.COASTLINE]:
        ax.add_feature(feature)
    ax.set_title(proj.__class__.__name__ + ' with extent [clip={}]'.format(clip))
fig.savefig(__file__[:-2]+'special.png')  # to make tight_layout working
fig.tight_layout(h_pad=1.5, w_pad=1.5, pad=0.1)
fig.savefig(__file__[:-2]+'special.png')

With this update, the position of labels is extrapolated from the intersection of a gridline with the map boundary, wherever it occurs. The visibility of each labels is checked so that the do not overlap themselve and do not overlap the map boundary. When an overlapping occurs, the other label directions are plotted within a given limit. Despite the algorithm can be improved, it should provide a decent basis.

Plots with global projections:

Special cases:

[corinnebosley]

Hi @stefraynaud, I've been taking a look at your work this morning. Those are some very pretty plots you've produced.

As you can see by the checks, I couldn't look at any test results because you need to resolve some merge conflicts before Appveyor can run them. I didn't find any test failures when I ran them on your branch locally though, so that bodes well.

There are a couple of sections in your code additions which I would like some clarification on, but I'll dig into that when we get some tests running.

Thank you for your continuing effort with this.

[pp-mo]

👍 just awesome !!

[crazyapril]

Any plans for a new release? This PR clears the ways for many developers to replace Basemap which will be EOL'ed just several months later. I think this PR alone can make a new release.

[QuLogic]

This doesn't appear to be used for anything.

[QuLogic]

There is no use of set_extent(..., clip=True), so how is this related to this PR? There are no examples or tests of it either, so I don't know what it's supposed to do either.

[stefraynaud]

You are right. It allows to adapt the map boundaries to specified lat/lon box, so you can make plots such as this one : #1117 (comment)

PR #1117 allows to add lon/lat labels even when map boundaries are not rectangular, like with clip=True if the projection is not rectangular.

Maybe this piece of code should be removed from master, and readded with another PR with a valid documentation.

[QuLogic]

But non-rectangular boundaries can be done with set_boundary, right? There are two examples for it.

[stefraynaud]

Setting clip=True makes two things

• compute de boundary according to the extent (here is the added value)
• apply the boundary set set_boundary.

[QuLogic]

Seems valuable, but I think it really should go through its own PR with documentation, examples, etc.

[QuLogic]

This change is incorrect; the CRS is not necessarily defined in terms of latitude and longitude.

[stefraynaud]

I didn't know about the curvilinear grid support of matplotlib tick labelling. Maybe it should be appropriately used for the cartopy labelling system. However, some adaptation would be needed to support labels placed around the edges of the map that do not coincide with the axes edges.

[zhonghua-zheng]

Hello, I installed the latest version
conda install -c conda-forge/label/dev cartopy

But I still got the error:

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-2-a98763482089> in <module>
     25     ax.set_extent(extent, crs=plain_crs)
     26     ax.coastlines(resolution='50m', color='k')
---> 27     ax.gridlines(color='lightgrey', linestyle='-', draw_labels=True)
     28 
     29 plt.show()

/opt/anaconda3/envs/only_plot/lib/python3.7/site-packages/cartopy/mpl/geoaxes.py in gridlines(self, crs, draw_labels, xlocs, ylocs, **kwargs)
   1222         gl = Gridliner(
   1223             self, crs=crs, draw_labels=draw_labels, xlocator=xlocs,
-> 1224             ylocator=ylocs, collection_kwargs=kwargs)
   1225         self._gridliners.append(gl)
   1226         return gl

/opt/anaconda3/envs/only_plot/lib/python3.7/site-packages/cartopy/mpl/gridliner.py in __init__(self, axes, crs, draw_labels, xlocator, ylocator, collection_kwargs)
    183         # public attributes are changed after instantiation.
    184         if draw_labels:
--> 185             self._assert_can_draw_ticks()
    186 
    187         #: The number of interpolation points which are used to draw the

/opt/anaconda3/envs/only_plot/lib/python3.7/site-packages/cartopy/mpl/gridliner.py in _assert_can_draw_ticks(self)
    397                             '{prj.__class__.__name__} plot.  Only PlateCarree'
    398                             ' and Mercator plots are currently '
--> 399                             'supported.'.format(prj=self.axes.projection))
    400         return True
    401 

TypeError: Cannot label gridlines on a RotatedPole plot.  Only PlateCarree and Mercator plots are currently supported.

I am wondering if I did something wrong?

[lukelbd]

Unfortunately it appears this feature has not yet been published in a release :/. See the version 0.18 milestone.

[QuLogic]

There is no 0.18.0b1 on https://anaconda.org/conda-forge/cartopy/labels, so you need to wait until they update.

[zhonghua-zheng]

Hi @lukelbd @QuLogic, thank you very much for your kind help!
I downloaded this repo, then run
python setup.py install
pip install Cartopy==0.18.0b1
Now it works.

[philippemiron]

I'm on v0.18.0 and getting:
set_extent() got an unexpected keyword argument 'clip'

is it suppose to be supported?

[QuLogic]

clip was not added, and is not an important part of this PR.

[philippemiron]

Do you know of any example how to manually add labels after using ax.set_boundary()? Thanks.

e.g.