Skip to content

MNT: update contour for mpl v3.8 #2213

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
greglucas merged 1 commit into from
Jul 17, 2023
Merged

MNT: update contour for mpl v3.8 #2213

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
7 changes: 3 additions & 4 deletions lib/cartopy/crs.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1247,7 +1247,8 @@ def quick_vertices_transform(self, vertices, src_crs):
is required (see :meth:`cartopy.crs.Projection.project_geometry`).

"""
return_value = None
if vertices.size == 0:
return vertices

if self == src_crs:
x = vertices[:, 0]
Expand All @@ -1258,9 +1259,7 @@ def quick_vertices_transform(self, vertices, src_crs):
y_limits = (self.y_limits[0] - epsilon, self.y_limits[1] + epsilon)
if (x.min() >= x_limits[0] and x.max() <= x_limits[1] and
y.min() >= y_limits[0] and y.max() <= y_limits[1]):
return_value = vertices

return return_value
return vertices


class _RectangularProjection(Projection, metaclass=ABCMeta):
Expand Down
126 changes: 76 additions & 50 deletions lib/cartopy/mpl/contour.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4,10 +4,11 @@
# See COPYING and COPYING.LESSER in the root of the repository for full
# licensing details.


import matplotlib as mpl
from matplotlib.contour import QuadContourSet
import matplotlib.path as mpath
import numpy as np
import packaging


class GeoContourSet(QuadContourSet):
Expand All @@ -20,66 +21,91 @@ class GeoContourSet(QuadContourSet):
# fiddling with instance.__class__.

def clabel(self, *args, **kwargs):
# nb: contour labelling does not work very well for filled
# contours - it is recommended to only label line contours.
# This is especially true when inline=True.

# This wrapper exist because mpl does not properly transform
# paths. Instead it simply assumes one path represents one polygon
# (not necessarily the case), and it assumes that
# transform(path.verts) is equivalent to transform_path(path).
# Unfortunately there is no way to easily correct this error,
# so we are forced to pre-transform the ContourSet's paths from
# the source coordinate system to the axes' projection.
# The existing mpl code then has a much simpler job of handling
# pre-projected paths (which can now effectively be transformed
# naively).

for col in self.collections:
# Snaffle the collection's path list. We will change the
# list in-place (as the contour label code does in mpl).
paths = col.get_paths()
if packaging.version.parse(mpl.__version__).release[:2] < (3, 8):
# nb: contour labelling does not work very well for filled
# contours - it is recommended to only label line contours.
# This is especially true when inline=True.

# This wrapper exist because mpl does not properly transform
# paths. Instead it simply assumes one path represents one polygon
# (not necessarily the case), and it assumes that
# transform(path.verts) is equivalent to transform_path(path).
# Unfortunately there is no way to easily correct this error,
# so we are forced to pre-transform the ContourSet's paths from
# the source coordinate system to the axes' projection.
# The existing mpl code then has a much simpler job of handling
# pre-projected paths (which can now effectively be transformed
# naively).

for col in self.collections:
# Snaffle the collection's path list. We will change the
# list in-place (as the contour label code does in mpl).
paths = col.get_paths()

# Define the transform that will take us from collection
# coordinates through to axes projection coordinates.
data_t = self.axes.transData
col_to_data = col.get_transform() - data_t

# Now that we have the transform, project all of this
# collection's paths.
new_paths = [col_to_data.transform_path(path)
for path in paths]
new_paths = [path for path in new_paths
if path.vertices.size >= 1]

# The collection will now be referenced in axes projection
# coordinates.
col.set_transform(data_t)

# Clear the now incorrectly referenced paths.
del paths[:]

for path in new_paths:
if path.vertices.size == 0:
# Don't persist empty paths. Let's get rid of them.
continue

# Split the path if it has multiple MOVETO statements.
codes = np.array(
path.codes if path.codes is not None else [0])
moveto = codes == mpath.Path.MOVETO
if moveto.sum() <= 1:
# This is only one path, so add it to the collection.
paths.append(path)
else:
# The first MOVETO doesn't need cutting-out.
moveto[0] = False
split_locs = np.flatnonzero(moveto)

split_verts = np.split(path.vertices, split_locs)
split_codes = np.split(path.codes, split_locs)

for verts, codes in zip(split_verts, split_codes):
# Add this path to the collection's list of paths.
paths.append(mpath.Path(verts, codes))

else:
# Where contour paths exist at the edge of the globe, sometimes a
# complete path in data space will become multiple paths when
# transformed into axes or screen space. Matplotlib's contour
# labelling does not account for this so we need to give it the
# pre-transformed paths to work with.

# Define the transform that will take us from collection
# coordinates through to axes projection coordinates.
data_t = self.axes.transData
col_to_data = col.get_transform() - data_t
col_to_data = self.get_transform() - data_t

# Now that we have the transform, project all of this
# collection's paths.
paths = self.get_paths()
new_paths = [col_to_data.transform_path(path) for path in paths]
new_paths = [path for path in new_paths if path.vertices.size >= 1]
paths[:] = new_paths
Copy link
Member Author

@rcomer rcomer Jul 14, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

set_paths is not implemented for the ContourSet although it is for other collections. I wonder if we should raise a feature request, although this use-case is pretty niche. I guess we would only have a problem if get_paths started returning a copy.

Copy link
Contributor

@greglucas greglucas Jul 14, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yeah, I think this could make sense. It does seem a bit suspect to rely on changing the returned object here like this.

Copy link
Member Author

@rcomer rcomer Jul 18, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Feature request now at matplotlib/matplotlib#26340.


# The collection will now be referenced in axes projection
# coordinates.
col.set_transform(data_t)

# Clear the now incorrectly referenced paths.
del paths[:]

for path in new_paths:
if path.vertices.size == 0:
# Don't persist empty paths. Let's get rid of them.
continue

# Split the path if it has multiple MOVETO statements.
codes = np.array(
path.codes if path.codes is not None else [0])
moveto = codes == mpath.Path.MOVETO
if moveto.sum() <= 1:
# This is only one path, so add it to the collection.
paths.append(path)
else:
# The first MOVETO doesn't need cutting-out.
moveto[0] = False
split_locs = np.flatnonzero(moveto)

split_verts = np.split(path.vertices, split_locs)
split_codes = np.split(path.codes, split_locs)

for verts, codes in zip(split_verts, split_codes):
# Add this path to the collection's list of paths.
paths.append(mpath.Path(verts, codes))
self.set_transform(data_t)

# Now that we have prepared the collection paths, call on
# through to the underlying implementation.
Expand Down
33 changes: 20 additions & 13 deletions lib/cartopy/mpl/geoaxes.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -45,7 +45,8 @@
from cartopy.mpl.slippy_image_artist import SlippyImageArtist


assert packaging.version.parse(mpl.__version__).release[:2] >= (3, 4), \
_MPL_VERSION = packaging.version.parse(mpl.__version__)
assert _MPL_VERSION.release >= (3, 4), \
'Cartopy is only supported with Matplotlib 3.4 or greater.'

# A nested mapping from path, source CRS, and target projection to the
Expand Down Expand Up @@ -1602,12 +1603,15 @@ def contour(self, *args, **kwargs):
result = super().contour(*args, **kwargs)

# We need to compute the dataLim correctly for contours.
bboxes = [col.get_datalim(self.transData)
for col in result.collections
if col.get_paths()]
if bboxes:
extent = mtransforms.Bbox.union(bboxes)
self.update_datalim(extent.get_points())
if _MPL_VERSION.release[:2] < (3, 8):
bboxes = [col.get_datalim(self.transData)
for col in result.collections
if col.get_paths()]
if bboxes:
extent = mtransforms.Bbox.union(bboxes)
self.update_datalim(extent.get_points())
else:
self.update_datalim(result.get_datalim(self.transData))
Copy link
Contributor

@greglucas greglucas Jul 14, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Do we need to update the datalim here or is that already taken care of because it is a full collection that is added to the axes now? The other collection methods like scatter and pcolormesh don't update the datalims from their return objects. Maybe this has to do with the contour wrapping stuff and not getting the bounds until actually calling update_datalim()? This is mostly just my curiosity here.

Copy link
Member Author

@rcomer rcomer Jul 14, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I saw several test failures when I tried not having it. But now I look more closely I am questioning whether they matter. For example

E           AssertionError: 
E           Arrays are not almost equal to 6 decimals
E           
E           Mismatched elements: 1 / 4 (25%)
E           Max absolute difference: 1.
E           Max relative difference: 0.00555556
E            x: array([-179.,  180.,  -25.,   25.])
E            y: array([-180,  180,  -25,   25])

which comes from the second assertion here. So it fails because the non-transformed-first version is now the same as the transformed-first version. Without much context, more consistency seems like an improvement?

cartopy/lib/cartopy/tests/mpl/test_contour.py

Lines 137 to 146 in 41880a8

# When calculating the contour in projection-space the extent
# will now be the extent of the transformed points (-179, 180, -25, 25)
test_func(xx, yy, z, transform=ccrs.PlateCarree(),
transform_first=True)
assert_array_almost_equal(ax.get_extent(), (-179, 180, -25, 25))
# The extent without the transform_first should be all the way out to -180
test_func(xx, yy, z, transform=ccrs.PlateCarree(),
transform_first=False)
assert_array_almost_equal(ax.get_extent(), (-180, 180, -25, 25))

There are also some image test failures which, by eye, look like they might be because of a similar extent change to the above.

We would need to add in something like self.ignore_existing_data_limits = False to cover my special dodgy latitude case from #2129.

Copy link
Contributor

@greglucas greglucas Jul 14, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ahh, OK. I think we definitely need this then, thanks for checking.

The two cases should definitely be different and that is preferable in these cases because their extents truly are different. The gallery kind of shows that if we force the extent to be global on the transformed-first one there is a gap and the extent is only -179/179 and doesn't touch the extreme edges. But, if we let the path logic do its magic it should be out to -180 and fill that gap.

https://scitools.org.uk/cartopy/docs/latest/gallery/scalar_data/contour_transforms.html#sphx-glr-gallery-scalar-data-contour-transforms-py


self.autoscale_view()

Expand Down Expand Up @@ -1650,12 +1654,15 @@ def contourf(self, *args, **kwargs):
result = super().contourf(*args, **kwargs)

# We need to compute the dataLim correctly for contours.
bboxes = [col.get_datalim(self.transData)
for col in result.collections
if col.get_paths()]
if bboxes:
extent = mtransforms.Bbox.union(bboxes)
self.update_datalim(extent.get_points())
if _MPL_VERSION.release[:2] < (3, 8):
bboxes = [col.get_datalim(self.transData)
for col in result.collections
if col.get_paths()]
if bboxes:
extent = mtransforms.Bbox.union(bboxes)
self.update_datalim(extent.get_points())
else:
self.update_datalim(result.get_datalim(self.transData))

self.autoscale_view()

Expand Down
5 changes: 4 additions & 1 deletion lib/cartopy/tests/mpl/test_examples.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,7 @@
import pytest

import cartopy.crs as ccrs
from cartopy.tests.mpl import MPL_VERSION


@pytest.mark.natural_earth
Expand All @@ -31,7 +32,9 @@ def test_global_map():


@pytest.mark.natural_earth
@pytest.mark.mpl_image_compare(filename='contour_label.png', tolerance=0.5)
@pytest.mark.mpl_image_compare(
filename='contour_label.png',
tolerance=3.9 if MPL_VERSION.release[:2] >= (3, 8) else 0.5)
def test_contour_label():
from cartopy.tests.mpl.test_caching import sample_data
fig = plt.figure()
Expand Down