Cartopy: Contour plot breaks with certain contour levels

[chengzhuzhang]

We encountered this seemingly random issue (only occurs with certain datasets) with Cartopy. Below is an example script. With less contour levels, script works:

levels =[-35, -30, -25, -20, -15, -10, 0, 10, 15, 20, 25, 30, 35]

However, with denser contour levels as below, it returns the following error : AttributeError: 'PreparedGeometry' object has no attribute 'is_valid'.

levels =[-35, -30, -25, -20, -15, -10, -5, -2, 2, 5, 10, 15, 20, 25, 30, 35]

No handlers could be found for logger "shapely.geos"
Traceback (most recent call last):
  File "cartopy_issue.py", line 30, in <module>
    extend='both',
  File "/Users/zhang40/anaconda2/envs/test_env/lib/python2.7/site-packages/cartopy/mpl/geoaxes.py", line 1339, in contourf
    if col.get_paths()])
  File "/Users/zhang40/anaconda2/envs/test_env/lib/python2.7/site-packages/matplotlib/collections.py", line 214, in get_datalim
    paths = [transform.transform_path_non_affine(p) for p in paths]
  File "/Users/zhang40/anaconda2/envs/test_env/lib/python2.7/site-packages/matplotlib/transforms.py", line 2386, in transform_path_non_affine
    return self._a.transform_path_non_affine(path)
  File "/Users/zhang40/anaconda2/envs/test_env/lib/python2.7/site-packages/cartopy/mpl/geoaxes.py", line 187, in transform_path_non_affine
    geom, self.source_projection)
  File "/Users/zhang40/anaconda2/envs/test_env/lib/python2.7/site-packages/cartopy/crs.py", line 175, in project_geometry
    return getattr(self, method_name)(geometry, src_crs)
  File "/Users/zhang40/anaconda2/envs/test_env/lib/python2.7/site-packages/cartopy/crs.py", line 330, in _project_polygon
    return self._rings_to_multi_polygon(rings, is_ccw)
  File "/Users/zhang40/anaconda2/envs/test_env/lib/python2.7/site-packages/cartopy/crs.py", line 541, in _rings_to_multi_polygon
    if prep_polygon.contains(interior_ring):
  File "/Users/zhang40/anaconda2/envs/test_env/lib/python2.7/site-packages/shapely/impl.py", line 37, in wrapper
    return func(*args, **kwargs)
  File "/Users/zhang40/anaconda2/envs/test_env/lib/python2.7/site-packages/shapely/prepared.py", line 45, in contains
    return bool(self.impl['prepared_contains'](self, other))
  File "/Users/zhang40/anaconda2/envs/test_env/lib/python2.7/site-packages/shapely/predicates.py", line 18, in __call__
    self._check_topology(err, this, other)
  File "/Users/zhang40/anaconda2/envs/test_env/lib/python2.7/site-packages/shapely/topology.py", line 34, in _check_topology
    if not geom.is_valid:
AttributeError: 'PreparedGeometry' object has no attribute 'is_valid'

Below is the code we used to generate this error. Do note that you'll need the netcdf4 library.

from netCDF4 import Dataset
import numpy as np
import matplotlib.pyplot as plt
from cartopy.util import add_cyclic_point
import cartopy.crs as ccrs

# Read test data
nc_fid = Dataset('test_data.nc','r')
lat = nc_fid.variables['lat'][:]
lon = nc_fid.variables['lon'][:]
var = nc_fid.variables['LWCF(test - reference)'][:]
cyclic_data, cyclic_lons = add_cyclic_point(var, coord=lon)

# Plot it
fig = plt.figure()
ax = fig.add_axes()

# These contour levels work fine...
levels =[-35, -30, -25, -20, -15, -10, 0, 10, 15, 20, 25, 30, 35]
# ...but slightly tweaking them results in a failure:
#levels =[-35, -30, -25, -20, -15, -10, -5, -2, 2, 5, 10, 15, 20, 25, 30, 35]

proj = ccrs.PlateCarree(central_longitude=180)
ax = plt.axes(projection= proj)
p1 = ax.contourf(cyclic_lons, lat, cyclic_data,
    transform=ccrs.PlateCarree(),
    norm=None,
    levels=levels,
    cmap='bwr',
    extend='both',
    )
ax.coastlines(lw=0.3)
ax.set_xticks([-180, -120, -60, 0, 60, 120, 180])
ax.set_yticks([-90, -60, -30, 0, 30, 60, 90])

# Color bar
cbax = fig.add_axes()
cbar = fig.colorbar(p1, cax=cbax)

plt.savefig('test.png', dpi=150)

We provided a zip file with both the script and data. We hope that this can be reproduced and fixed by the Cartopy team.
cartopy_issue.v2.zip

More detailed discussion can be found from our git repository:
https://github.com/ACME-Climate/acme_diags/issues/98

Let us (me, @zshaheen, or @golaz) know if we can assist with anything.

[pelson]

Thanks for supplying the data - certainly makes it easer reproducing these kinds of issues.

We've made a lot of improvements to the geometry transform code in cartopy since the last release. It would be good to test this against master or v0.16 (when it is released).

[akrherz]

Sorry for the misfire on the previous comment, indeed, I reproduce this on current conda-forge and cartopy 0.16

[BenMGeo]

I have the same issue as described here. I'm calculating percentiles and a specific level for a specific data set cannot be plotted. The contour-levels are [0.00, 0.04, 0.07, 0.11, 0.15, 0.18, 0.22, 0.26, 0.30, 0.33, 0.37]. If I add ,e.g., 0.00001 to those levels, as in [1e-05, 0.040010000000000004, 0.07001, 0.11001, 0.15001, 0.18001, 0.22001, 0.26001, 0.30001, 0.33001, 0.37001], I get the plots done. (I added a try-except environment to my code, but this is an ugly workaround.)

Looks like a pretty random bug, to me.

Edit: 0.16 conda -forge cartopy

[akrherz]

I had another look at this and as I debug, I see this logged
TopologyException: side location conflict at -180 0.11677411238420946

if I make this subtle code change to your example, the script works
```diff --git a/cartopy/issue946.py b/cartopy/issue946.py
index e036040..1929b0c 100644
--- a/cartopy/issue946.py
+++ b/cartopy/issue946.py
@@ -20,7 +20,7 @@ ax = fig.add_axes()
# ...but slightly tweaking them results in a failure:
levels =[-35, -30, -25, -20, -15, -10, -5, -2, 2, 5, 10, 15, 20, 25, 30, 35]

-proj = ccrs.PlateCarree(central_longitude=180)
+proj = ccrs.PlateCarree(central_longitude=179.99)
ax = plt.axes(projection= proj)
p1 = ax.contourf(cyclic_lons, lat, cyclic_data,
transform=ccrs.PlateCarree(),
```
shrug, I continue to attempt to grok issues like these, but half of the time I wonder if we are just seeing flaky issues due to stderr exception logging that is not properly handled.

[BenMGeo]

Update: I realized that if I set the contours to levels that produce a plot object ( GeoAxes: ), I still cannot print the object. I think recalculations to the preferred projection are taking place that revise contours and cannot be handled.

[jonnyhtw]

Hi there

I am also seeing this error (FYI @pelson, @pdearnshaw it's the Iris-based climate validation note code); seemingly at random!

Is there any news on a fix?

happy to help with testing etc.

Cheers

Jonny

[zshaheen]

I tested the above cartopy_issue.v2.zip with cartopy 0.16.0, build py27h81b52dc_2 from the conda-forge channel and I'm still having an issue.

[jonnyhtw]

same...

thanks for the info

jonny