Search code examples
pythonmatplotlibcartopy

Ignore projection limits when setting the extent


Is there a way to set the extent of the figure beyond the projection limits?

For example, when using the "Rijksdriehoek" projection (EPSG 28992), the limits from Cartopy (proj4?) are wrong, way too narrow.

That projection is designed cover all of the Netherlands, but the imposed limits even cause part of the country to be cut off. Whereas I would rather set the extent slightly wider then the official boundaries to provide some extra context.

Rd example

Unfortunately, the set_extent method gives an error:

ValueError: Failed to determine the required bounds in projection 
coordinates. Check that the values provided are within the valid range 
(x_limits=[646.3608848793374, 284347.25011780026], 
y_limits=[308289.55751689477, 637111.0245778429]).

The set_xlim/set_ylim methods don't seem to do anything, which would work for a normal matplotlib axes.

Example code:

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

projection = ccrs.epsg(28992)

fig, ax = plt.subplots(figsize=(5,10), subplot_kw=dict(projection=projection))

ax.coastlines(resolution='10m')
ax.add_feature(cfeature.NaturalEarthFeature('cultural', 'admin_0_boundary_lines_land', '10m', facecolor='none', edgecolor='k'))

The extent of the figure is automatically set to the limits of the projection:

print(projection.bounds)
print(ax.get_extent())

(646.3608848793374, 284347.25011780026, 308289.55751689477, 637111.0245778429)
(646.3608848793374, 284347.25011780026, 308289.55751689477, 637111.0245778429)

According to the documentation about the projection, the actual limits should be: (-700 300000 289000 629000). But even those seem a bit strict for visualization purposes.

See for example the "Scope of validity section" at:

https://translate.google.com/translate?hl=en&sl=nl&u=https://nl.wikipedia.org/wiki/Rijksdriehoeksco%25C3%25B6rdinaten


Solution

  • The answer by @pp-mo is excellent. However, here is an alternative solution. The working code is:

    import matplotlib.pyplot as plt
    import cartopy.crs as ccrs
    import cartopy.feature as cfeature
    
    # subclassing for a modified projection
    class nd_prj(ccrs.Stereographic):
        """
        Hacked projection for ccrs.epsg(28992) to get extended plotting limits
        """
        def __init__(self):
            globe = ccrs.Globe(ellipse=u'bessel')
            super(nd_prj, self).__init__(central_latitude=52.15616055555555, \
                         central_longitude=5.38763888888889, \
                         #true_scale_latitude=52.0, \
                         scale_factor=0.9999079, \
                         false_easting=155000, false_northing=463000, globe=globe)
    
        @property
        def x_limits(self):
            return (500, 300000)   # define the values you need
    
        @property
        def y_limits(self):
            return (300000, 650000) # define the values you need
    
    projection = nd_prj()  # make use of the projection
    fig, ax = plt.subplots(figsize=(5,10), subplot_kw=dict(projection=projection))
    
    ax.coastlines(resolution='10m')
    ax.add_feature(cfeature.NaturalEarthFeature('cultural', 'admin_0_boundary_lines_land', \
                                                '10m', facecolor='none', edgecolor='k'))
    plt.show()
    

    The resulting plot:

    enter image description here

    Hope this is useful.