GeoAxes.set_extent on non-cylindrical projections

[QuLogic]

I think I am not the only one who finds it a bit confusing that while GeoAxes.set_extent accepts a crs, it only really has an effect on the "corners". That is, it essentially re-projects the corners and then takes the largest x/y boundary that encases them. It is probably fine in something like Mercator or PlateCarree where latitude and longitude are at right-angles.

But if you pick something conical like Albers or Lambert Conformal, then lines of constant latitude or longitude are not parallel to the figure edges:

There is a relevant example for how to set a more complex boundary, but there are various subtleties. For example, you cannot just create a rectangle because of #363. I can add this snippet to go along with the star example, but I'm wondering if this should be available as a builtin method?

    npoints_side = npoints - 1
    verts = np.empty((npoints_side * 4 + 1, 2))
    # left-bottom to right-bottom
    verts[:npoints, 0] = np.linspace(min_lon, max_lon, npoints)
    verts[:npoints, 1] = min_lat
    # right-bottom to right-top
    verts[npoints_side:npoints_side + npoints, 0] = max_lon
    verts[npoints_side:npoints_side + npoints, 1] = np.linspace(min_lat, max_lat, npoints)
    # right-top to left-top
    verts[2 * npoints_side:2 * npoints_side + npoints, 0] = np.linspace(max_lon, min_lon, npoints)
    verts[2 * npoints_side:2 * npoints_side + npoints, 1] = max_lat
    # left-top to left-bottom
    verts[3 * npoints_side:, 0] = min_lon
    verts[3 * npoints_side:, 1] = np.linspace(max_lat, min_lat, npoints)

    codes = mpath.Path.LINETO * np.ones(verts.shape[0])
    codes[0] = mpath.Path.MOVETO
    codes[-1] = mpath.Path.CLOSEPOLY

    rect = mpath.Path(verts, codes)
    ax.set_extent([min_lon, max_lon, min_lat, max_lat]) 
    ax.set_boundary(rect, transform=proj.as_geodetic())

[pelson]

I agree the behaviour may not be immediately obvious. One advantage of keeping the axes rectangular is that we can make use of the pan/zoom functionality given to us by matplotlib. If we took the extent provided, and used that as the axes extent limits explicitly, we are precluding that functionality.

As for the example you gave, I agree - we can make this a huge amount easier with very little effort. For the record, your example can be simplified if we just make use of matplotlib's Path machinery (much of which was added to for cartopy). Your example drops out to be:

import cartopy.crs as ccrs
import cartopy.feature as feat
import matplotlib.pyplot as plt
import matplotlib.path as mpath


ax = plt.axes(projection=ccrs.LambertConformal())
ax.gridlines()
ax.add_feature(feat.NaturalEarthFeature('physical', 'land', '50m',
                                        facecolor=feat.COLORS['land'],
                                        edgecolor='black',
                                        linewidth=1.2))
xlim = [-120, -60]
ylim = [60, 80]

rect = mpath.Path([[xlim[0], ylim[0]],
                   [xlim[1], ylim[0]],
                   [xlim[1], ylim[1]],
                   [xlim[0], ylim[1]],
                   [xlim[0], ylim[0]],
                   ]).interpolated(20)

proj_to_data = ccrs.PlateCarree()._as_mpl_transform(ax) - ax.transData
rect_in_target = proj_to_data.transform_path(rect)

ax.set_boundary(rect_in_target)

# Notice the ugly hack to stop any further clipping - this is
# the same problem as #363.
ax.set_extent([xlim[0], xlim[1], ylim[0] - 2, ylim[1]])
plt.show()

[pelson]

Proposal:

• Extend both set_extent and set_boundary to accept shapely geometries
• Extend set_boundary to accept extents as per set_extent.
• Extend set_boundary to allow passing a CRS to transform

Because of argument renaming, this would be a breaking change.

[kdpenner]

Since cartopy has changed significantly since this issue was opened, I want to clarify a few things.

set_extent does interpolate between the corners:

https://github.com/SciTools/cartopy/blob/2148e13846944475456390d738b9cbb4391c3add/lib/cartopy/mpl/geoaxes.py#L857-L858

which leads to:

https://github.com/SciTools/cartopy/blob/2148e13846944475456390d738b9cbb4391c3add/lib/cartopy/crs.py#L206-L207

which leads to:

https://github.com/SciTools/cartopy/blob/d12c86ca228365c61d5fb8b7121b49201f19fea8/lib/cartopy/trace.pyx#L627

The problem behind #1367 is that the interpolator doesn't use enough points:

The issue behind #1362 and an issue I came across with transverse Mercator is that interpolation of the bounding geometry when a Plate Carree CRS is passed does not lead to an extent in projected space that most people expect.

Note that [-180, 180, -90, 90] does not work, perhaps because of a numerical issue. pyproj projects (-180, -90) to (-4.8e-26, -1e7).

When I pass a region in Plate Carre corresponding to the globe, I expect a plot in projected space with as much as possible of the globe.

The issue behind #1617 is indeed that no interpolation is done.

[pgf]

Continuing here #1804. Although this behavior is not a set_extent bug, I would like to point out a possible solution. As pointed out by @pelson the default rectangular axes has the advantage of the pan/zoom functionality. Anyway, if one just wants to produce some publication quality map where the pan/zoom functionality is not relevant, a workaround is to avoid set_extent and directly use set_xlim/set_ylim on the boundary extremes:

Code to reproduce

import matplotlib.pyplot               as plt
import matplotlib.path                 as mpath
import cartopy.crs                     as ccrs
import cartopy.mpl.ticker              as ctk


lon1, lon2, lat1, lat2 = [-20, 20, 50, 80]

rect = mpath.Path([[lon1, lat1], [lon2, lat1],
    [lon2, lat2], [lon1, lat2], [lon1, lat1]]).interpolated(50)

name='NearsidePerspective'
proj=ccrs.NearsidePerspective(central_longitude=(lon1+lon2)*0.5,
    central_latitude=(lat1+lat2)*0.5)

fig, (ax1, ax2) = plt.subplots(1,2, subplot_kw={'projection':proj,})

proj_to_data = ccrs.PlateCarree()._as_mpl_transform(ax1) - ax1.transData
rect_in_target = proj_to_data.transform_path(rect)
ax1.set_boundary(rect_in_target)
ax1.set_extent([lon1, lon2, lat1, lat2], crs=ccrs.PlateCarree())
ax1.coastlines()

gl=ax1.gridlines(draw_labels=True, x_inline=False, y_inline=False, linestyle='dashed')
gl.top_labels=False
gl.right_labels=False
gl.rotate_labels=False
gl.xlocator=ctk.LongitudeLocator(4)
gl.ylocator=ctk.LatitudeLocator(6)
gl.xformatter=ctk.LongitudeFormatter(zero_direction_label=True)
gl.yformatter=ctk.LatitudeFormatter()

ax1.set_title(name+'\nset_extent()')

proj_to_data = ccrs.PlateCarree()._as_mpl_transform(ax2) - ax2.transData
rect_in_target = proj_to_data.transform_path(rect)
ax2.set_boundary(rect_in_target)
ax2.set_xlim(rect_in_target.vertices[:,0].min(), rect_in_target.vertices[:,0].max())
ax2.set_ylim(rect_in_target.vertices[:,1].min(), rect_in_target.vertices[:,1].max())
ax2.coastlines()

gl=ax2.gridlines(draw_labels=True, x_inline=False, y_inline=False, linestyle='dashed')
gl.top_labels=False
gl.right_labels=False
gl.rotate_labels=False
gl.xlocator=ctk.LongitudeLocator(4)
gl.ylocator=ctk.LatitudeLocator(6)
gl.xformatter=ctk.LongitudeFormatter(zero_direction_label=True)
gl.yformatter=ctk.LatitudeFormatter()

ax2.set_title(name+'\nset_xlim()/set_ylim()')

fig.tight_layout()
plt.show()

[kdpenner]

@pgf your situation is similar to the one I wrote about for #1367 above. The interpolator in trace.pyx doesn't sample enough points to determine the extent in projected space from the bounds of the box encompassing all the points.

[pgf]

@pgf your situation is similar to the one I wrote about for #1367 above. The interpolator in trace.pyx doesn't sample enough points to determine the extent in projected space from the bounds of the box encompassing all the points.

Concerning your example #1367 , again this kind of behaviour is better addressed using set_xlim/set_ylim insted of set_extent.

Code to reproduce

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

def fix_extent(ax, extent):
    mlon = np.mean(extent[:2])
    mlat = np.mean(extent[2:])
    xtrm_data    = np.array([[extent[0], mlat], [mlon, extent[2]], [extent[1], mlat], [mlon, extent[3]]])
    proj_to_data = ccrs.PlateCarree()._as_mpl_transform(ax) - ax.transData
    xtrm         = proj_to_data.transform(xtrm_data)
    ax.set_xlim(xtrm[:,0].min(), xtrm[:,0].max())
    ax.set_ylim(xtrm[:,1].min(), xtrm[:,1].max())

extent=(-180, 180, -60, 60)

ax = plt.axes(projection=ccrs.Robinson())
ax.coastlines()
ax.gridlines()
fix_extent(ax, extent)

plt.show()

[kdpenner]

again this kind of behaviour is better addressed using set_xlim/set_ylim insted of set_extent

I disagree---your work and #1367 reveal a bug in set_extent that may be solved by increasing the number of points sampled by trace.pyx. What you've presented is a workaround that should be unnecessary.

[pgf]

I disagree---your work and #1367 reveal a bug in set_extent that may be solved by increasing the number of points sampled by trace.pyx. What you've presented is a workaround that should be unnecessary.

After spending more time delving into this set_extent behaviour I agree that my workaround should be unnecessary. Why don't let set_extent use more points instead of just relying on the four corners? Following this approach I experimented a bit with how set_extent defines domain_in_crs and found a possible fix replacing sgeom.polygon.LineString in: https://github.com/SciTools/cartopy/blob/82a6b66f609cb1102d62cc4de6c3cf89951a7b01/lib/cartopy/mpl/geoaxes.py#L836-L839

with:

x1, x2, y1, y2 = extents 
rect = mpath.Path([[x1, y1], [x2, y1],
    [x2, y2], [x1, y2], [x1, y1]], closed=True).interpolated(50)

domain_in_crs = cpatch.path_to_geos(rect)[0]

This change gives the same result as my previous workaround but using the usual set_extent syntax. I also tested it on a rectangular region using several projections (see #1804) and it works. Could this be a valid fix? Can you see any possible drawback?

For the sake of completeness I would like to add one more information. The request of a closed Path generates a Polygon (domain_in_crs.geom_type == 'Polygon'). This is handled in crs.py by project_geometry through _project_polygon . With no closure request path_to_geos generates a MultiLineString, which is handled through _project_multiline. Oddly, in the latter case the fix doesn't work and gives a weird result.

[dopplershift]

@pgf Your solution using Path and calling interpolated() is pretty much functionally solving the problem identified by @kdpenner: that there are insufficient points being interpolated. Your solution just completely goes around it by using Matplotlib's Path.

I agree with @kdpenner that the real solution would probably be to make _project_line_string or interpolator use more points. I don't know this code well, @kdpenner is there an easy way to make this a tunable knob, at least internally? To me, that's the best fix.

Failing that, I'm completely fine with the approach to use Matplotlib's Path to get the interpolated polygon if it fixes the issue and doesn't break anything--I can't think of any reason why it would break, all it really is doing is cheating and using Matplotlib's support for interpolating a Path to do what CartoPy should already be doing.

[kdpenner]

@dopplershift with the caveat that I know no cython/c++ and haven't used GEOS directly: the important parameter may be threshold in trace.pyx. Related to this TODO?

https://github.com/SciTools/cartopy/blob/d12c86ca228365c61d5fb8b7121b49201f19fea8/lib/cartopy/trace.pyx#L446-L447

threshold is defined here:

https://github.com/SciTools/cartopy/blob/d12c86ca228365c61d5fb8b7121b49201f19fea8/lib/cartopy/trace.pyx#L614

so a first pass at this would be to alter the projection's threshold property. I may be able to do this on the timescale of a month or so.

[dopplershift]

@kdpenner Well, an easy way to test might be #1815 .

[kdpenner]

excellent timing.

FWIW #1739 should fix the point I made earlier in this thread about #1362.

[pgf]

@kdpenner @dopplershift I've already experimented with the threshold value with no success. I've got some improvement reducing the value at: https://github.com/SciTools/cartopy/blob/d12c86ca228365c61d5fb8b7121b49201f19fea8/lib/cartopy/trace.pyx#L449 Reducing from 1.0e-6 to 1.0e-8 results in some more points, but not enough to fix this issue: Original value, 1.0e-6: 1.0e-8: Smaller values do not add new points unfortunatelly.

[kdpenner]

I didn't clone #1815 but I was able to produce the following plot by altering the following line, i.e. by changing only threshold:

https://github.com/SciTools/cartopy/blob/7c75f5ab6a397bca3197774b063820b6aabe87f9/lib/cartopy/crs.py#L2028

[blazing216]

@dopplershift with the caveat that I know no cython/c++ and haven't used GEOS directly: the important parameter may be threshold in trace.pyx. Related to this TODO?

https://github.com/SciTools/cartopy/blob/d12c86ca228365c61d5fb8b7121b49201f19fea8/lib/cartopy/trace.pyx#L446-L447

threshold is defined here:

https://github.com/SciTools/cartopy/blob/d12c86ca228365c61d5fb8b7121b49201f19fea8/lib/cartopy/trace.pyx#L614

so a first pass at this would be to alter the projection's threshold property. I may be able to do this on the timescale of a month or so.

I have been looking into threshold during the last couple of days. As far as I understand, threshold basically determines how far (defined as distance in projected coordinates; see straightAndDomain in trace.pyx) each projected line segment can deviate from the perfectly interpolated curve.

_project_segment in trace.pyx projects each line segment (between two vertices; a two-point line) to a multiple-point line. When threshold is significantly small, the projected multiple-point line follows the 'perfect interpolation'. When threshold is too large, the projected multiple-point line will have only two end points, the same as no interpolation at all.

So reducing threshold can solve the seemingly problem of set_extent. Currently, the default threshold for each projection seems too high. An adaptive threshold might be the ultimate solution. #8

[kdpenner]

@blazing216 yup, that's my understanding as well. Can you tell if threshold has the same meaning for a LinearRing?

In this thread we are talking about projecting a LineString, but in #1739 I replaced the LineString here:

https://github.com/SciTools/cartopy/blob/f019487b3cb3717bf5d56a8b23097708ec3a1b8b/lib/cartopy/mpl/geoaxes.py#L837-L839

with a Polygon to solve a different issue with set_extent. What's actually projected for a Polygon are the interior and exterior LinearRings.

[pgf]

@kdpenner Just to clarify, for the Robinson projection I've changed: https://github.com/SciTools/cartopy/blob/7c75f5ab6a397bca3197774b063820b6aabe87f9/lib/cartopy/crs.py#L1874-L1877 to:

        self._threshold = 1e4

    @property
    def threshold(self):
        return self._threshold

and then modified #1367 with:

print(ax.projection._threshold)
ax.projection._threshold=1.e-2
print(ax.projection._threshold)
ax.set_extent(extent, ccrs.PlateCarree())

Trying several values didn't fix this particular case.

[kdpenner]

@pgf thanks for changing threshold on the Robinson projection. In my comment I wrote too hastily in attributing all of the issue to not enough interpolates. My version of your problematic plot has better longitude bounds than does yours:

My development branch fixes the issue:

and I haven't changed threshold there. The problem was probably fixed by changing the extent geometry to a Polygon from a LineString.

[pgf]

and I haven't changed threshold there. The problem was probably fixed by changing the extent geometry to a Polygon from a LineString.

@kdpenner This makes sense. Indeed sounds like my previous experience:

For the sake of completeness I would like to add one more information. The request of a closed Path generates a Polygon (domain_in_crs.geom_type == 'Polygon'). This is handled in crs.py by project_geometry through _project_polygon . With no closure request path_to_geos generates a MultiLineString, which is handled through _project_multiline. Oddly, in the latter case the fix doesn't work and gives a weird result.

Maybe the issues is related to differences between _project_polygon and _project_multiline?

[blazing216]

@blazing216 yup, that's my understanding as well. Can you tell if threshold has the same meaning for a LinearRing?

@kdpenner I think so. Under the hood, both LinearRing and LineString use project_linear. The difference is that _project_line_string is a thin wrapper of project_linear, while _project_linear_ring first projects LinearRing to a MultipleLineString in the projected coordinates, then forms LinearRing again if possible.

project_linear seems the root of a lot of problems, including #1797, which drives me to look into trace.pyx.

See the below the beginning of _project_linear_ring which calls project_linear

def _project_linear_ring(self, linear_ring, src_crs):
    """
    Project the given LinearRing from the src_crs into this CRS and
    returns a list of LinearRings and a single MultiLineString.

    """
    debug = False
    # 1) Resolve the initial lines into projected segments
    # 1abc
    # def23ghi
    # jkl41
    multi_line_string = cartopy.trace.project_linear(linear_ring,
                                                     src_crs, self)

[kdpenner]

The difference is that _project_line_string is a thin wrapper of project_linear, while _project_linear_ring first projects LinearRing to a MultipleLineString in the projected coordinates, then forms LinearRing again if possible.

OK, to clarify, _project_linear_ring passes a LinearRing to trace.pyx. The conversion to a MultiLineString happens in trace:

https://github.com/SciTools/cartopy/blob/7c75f5ab6a397bca3197774b063820b6aabe87f9/lib/cartopy/trace.pyx#L641-L645

The reason that an extent defined as a Polygon differs from the same extent defined as a LineString or a LinearRing is _attach_lines_to_boundary in _project_polygon:

https://github.com/SciTools/cartopy/blob/4ca93c1eefa86cb85177467407852b4636bb3e70/lib/cartopy/crs.py#L351-L352

project_linear seems the root of a lot of problems, including #1797, which drives me to look into trace.pyx.

Oops, I looked into #1797 and found a different way to fix the problem. I will follow up in that thread.

[blazing216]

I found set_extent already has everything we need to set an 'intuitive' extent. https://github.com/SciTools/cartopy/blob/7c75f5ab6a397bca3197774b063820b6aabe87f9/lib/cartopy/mpl/geoaxes.py#L857-L863

I added the following lines before calling set_xlim and set_ylim. The new set_extent can set a rectilinear boundary if calling with crs=ccrs.PlateCarree() and a rectangle boundary if calling with crs=None (default).

The advantage of the new set_extent is that we can avoid the hack to manually adjust the xlim and ylim to include the entire boundary. The bounds of the boundary is computed by x1, y1, x2, y2 = projected.bounds.

It also does not break the pan/zoom functionality. But we cannot see the content beyond the boundary.

I think we have to think about this: what do we want set_extent to do, or what do USERS expect it to do?

if crs is not None:
            path, = cpatch.geos_to_path(projected)
            self.patch.set_boundary(path, self.transData)
            self.spines['geo'].set_boundary(path, self.transData)

Code to produce the example:

import cartopy.crs as ccrs
import cartopy.feature as feat
import matplotlib.pyplot as plt
import matplotlib.path as mpath
import shapely.geometry as sgeom
import cartopy.mpl.patch as cpatch
import matplotlib.patches as patches

high_res_proj = ccrs.LambertConformal()
high_res_proj.threshold = 1e3

xlim = [-120, -60]
ylim = [60, 80]

# with smaller threshold, it is no need to interpolate
rect = mpath.Path([[xlim[0], ylim[0]],
                   [xlim[1], ylim[0]],
                   [xlim[1], ylim[1]],
                   [xlim[0], ylim[1]],
                   [xlim[0], ylim[0]],
                   ])#.interpolated(20)

plt.figure(figsize=(6*2,2*2))                                    

ax = plt.subplot(131, projection=high_res_proj)
ax.gridlines()
ax.add_feature(feat.NaturalEarthFeature('physical', 'land', '50m',
                                        facecolor=feat.COLORS['land'],
                                        edgecolor='black',
                                        linewidth=1.2))

proj_to_data = ccrs.PlateCarree()._as_mpl_transform(ax) - ax.transData
rect_in_target = proj_to_data.transform_path(rect)
ax.add_patch(patches.PathPatch(rect_in_target, edgecolor='r', 
    facecolor='none', lw=4, zorder=2))

ax.set_boundary(rect_in_target)

# Notice the ugly hack to stop any further clipping - this is
# the same problem as #363.
ax.set_extent([xlim[0], xlim[1], ylim[0] - 2, ylim[1]])

plt.title('Path')

ax = plt.subplot(132, projection=high_res_proj)
ax.gridlines()
ax.add_feature(feat.NaturalEarthFeature('physical', 'land', '50m',
                                        facecolor=feat.COLORS['land'],
                                        edgecolor='black',
                                        linewidth=1.2))

ax.add_patch(patches.PathPatch(rect_in_target, edgecolor='r', 
    facecolor='none', lw=4, zorder=2))
ax.set_extent([xlim[0], xlim[1], ylim[0], ylim[1]])

plt.title('New set_extent(crs=None)')


ax = plt.subplot(133, projection=high_res_proj)
ax.gridlines()
ax.add_feature(feat.NaturalEarthFeature('physical', 'land', '50m',
                                        facecolor=feat.COLORS['land'],
                                        edgecolor='black',
                                        linewidth=1.2))
ax.add_patch(patches.PathPatch(rect_in_target, edgecolor='r', 
    facecolor='none', lw=4, zorder=2))

ax.set_extent([xlim[0], xlim[1], ylim[0], ylim[1]],
    crs=ccrs.PlateCarree())

plt.title('New set_extent(crs=PlateCarree())')

plt.tight_layout()
plt.savefig('boundary_using_set_extent.png', dpi=150)

[kdpenner]

if crs is not None:
            path, = cpatch.geos_to_path(projected)
            self.patch.set_boundary(path, self.transData)
            self.spines['geo'].set_boundary(path, self.transData)

@blazing216 where do you add the code above? It breaks a script of mine if I add it here:

https://github.com/SciTools/cartopy/blob/7c75f5ab6a397bca3197774b063820b6aabe87f9/lib/cartopy/mpl/geoaxes.py#L871

[blazing216]

Yes, I added to where you pointed. I am just sharing an idea that you can add set_boundary in set_extent to implement the curvilinear boundary. Probably it is not a universal solution.

[KeithWM]

Old issue, I know, but I think I'm running into it too when trying to create a print of a subset of some datapoints. I want to be able to list the datapoints that should be included in a zoomed-in map and to create a map that is big enough to show those points (but not more than that). I'm using the Peirce Quincuncial projection (because why not...), but then it is a real pain to set the boundaries in terms of min-max lon-lat.

Does

projected = self.projection.project_geometry(domain_in_crs, crs)
[...]
x1, y1, x2, y2 = projected.bounds

work if the domain_in_crs is something else than the sgeom.polygon.LineString representation of lon-lat square that is currently implemented? I really like the idea of being able to call ax.set_extent with the first argument being a collection of points that should be included.