Update geoaxes.py pcolor to catch the ValueError #1476

[htonchia]

pcolor crash #1476: catch the ValueError and remove the artist

Rationale

In some cases, pcolor raises a ValueError and crashes when plotting data from one projection to an other see #1476

It comes from the fact that it cannot transform one of the paths to the axis projection.
Instead of crashing the plot, this modification catches the error in pcolor and removes the faulty artist from the axis to avoid further crashes in draw().

Only once this PR is accepted, it will be possible to add a pcolor crash test in the test suite.

A better result would be obtained by plotting each path (grid cell) in turn and only removing the cells that crashes pcolor. That would be very time consuming see #1476 for an example.

This is a modification already proposed in #1467 from which it can be removed to simplify it, this is a first baby step.

Implications

The only implication to my knowledge is that it doesn't crash.

[dopplershift]

So if I understand this correctly, currently we have cases where calls to pcolor/pcolormesh result in draw-time failures and this PR converts those to a warning?

I feel like the error is the proper behavior, otherwise, since users frequently do not read or pay attention to warnings, we're going to get bug reports about missing plots.

I feel like the right solution is:

1. Fix pcolor to work with these data
2. Do better checking of data when calling pcolor so that the error occurs up-front rather than at draw-time.

[htonchia]

@dopplershift, yes we have draw time failures with pcolor, for cells that cannot be plotted by pcolor.
As pcolormesh uses pcolor for wrapping cells it cannot plot, the issues happens also with pcolormesh.

As you suggest we have three options:

1. crashing entirely the plot as of now
2. not plotting pcolor but still plotting the other artist as suggested by the PR
3. 'fixing matlplolib and shapely' to remove the crash

As to do it upfront, the quickest way to find the cells it cannot plot is to try to draw them as the best way to know which cells are not compatible with the source and destination projections using shapely is to go through the transformation with shapely, which is done at draw time.

Below is a plot with cells that can be plotted with pcolor (lon starting at -175) and cells that cannot be plotted with colors (lon starting at -176).

The crashes happens in shapely: A LinearRing must have at least 3 coordinate tuples. One way to go further we would have to modify shapely and matplotib. I haven't found any obvious difference between a cell that crashes and one that doesn't.

The advantage of 2 is for pcolormesh plot. In that case only a couple of cells at the boundary are sent to pcolor for plotting to manage wrapping and they crash the whole plot with A LinearRing must have at least 3 coordinate tuples and the person has no way to guess or understand what happens or to plot one's data.

I feel that 2 is a progress as we go from a crash to a partial plot that would be totally acceptable for pcolormesh case.

There may be a 2.5 way to deal with it. After removing the artist that fails, it could be possible to try to plot each cell in turn and only remove those that crashes at draw time. I am afraid it may be quite time consuming but it may be worth testing.

[htonchia]

Following the above comment, the following code removes only the cells that raises the value error A LinearRing must have at least 3 coordinate tuples. It would be the 2.5 way to deal with the issue.
To make it work, it uses the geoaxes.py from the PR with a slight modification to keep the artist that crashes and raises the value error so it can be caught and process by the code in the example.
This is the modified geoaxes code:

        try:
            result = matplotlib.axes.Axes.pcolor(self, *args, **kwargs)

            # Update the datalim for this pcolor.
            limits = result.get_datalim(self.transData)
            self.update_datalim(limits)
            self.autoscale_view()
        except ValueError as err:
            warnings.warn("Unexpected ValueError in pcolor"
                          " the pcolor plot artist is removed."
                          " Some cells could not be transformed"
                          " in the axis projection.")
            for child in self.get_children():
                if child is result:
                    self.crashedpcolorchild = child
                    child.remove()

            raise ValueError(err)

This is the code of an example that removes only the artist of the cells that crashes shapely.

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

def comparecoll(child, mesh):
    equal = True
    for c in child.properties():
        if c in ['array', 'paths', 'norm', 'cmap' ]:
            if c == 'paths':
                for pc, pm in zip(child.properties()[c], mesh.properties()[c]):
                    if (pc.vertices != pm.vertices).any():
                        equal = False
                        break
                    if (pc.codes != pm.codes).any():
                        equal = False
                        break
            else:
                if not ((child.properties()[c] == mesh.properties()[c]) is False):
                    continue
                else:
                    print(child.properties()[c])
                    print(mesh.properties()[c])
                    equal = False
                    break
    return equal

proj = ccrs.LambertAzimuthalEqualArea

lats = np.linspace(10, 42, 3)
lonsok = np.linspace(-175, -90, 3)
lonsno = np.linspace(-176, -90, 3)
lonsok, latsok = np.meshgrid(lonsok, lats)
lonsno, latsno = np.meshgrid(lonsno, lats)

axproj = proj(central_latitude=-30)

xyeok = axproj.transform_points(ccrs.PlateCarree(), lonsok, latsok)
xyeno = axproj.transform_points(ccrs.PlateCarree(), lonsno, latsno)

xeok = xyeok[..., 0]
yeok = xyeok[..., 1]
xeno = xyeno[..., 0]
yeno = xyeno[..., 1]

np.set_printoptions(precision=0)
np.set_printoptions(precision=3)

z4 = latsok[:-1, :-1] * (1 + 0.01 * lonsok[:-1, :-1])
fig = plt.figure(figsize=(13.1, 6.4), dpi=100)
plt.tight_layout()
fig.suptitle('LambertAzimuthalEqualArea central_latitude=-30')

ax21 = plt.subplot(2, 2, 3, projection=ccrs.PlateCarree())
ax22 = plt.subplot(2, 2, 4, projection=ccrs.PlateCarree())

ax1 = plt.subplot(2, 4, 1, projection=axproj)
ax2 = plt.subplot(2, 4, 2, projection=axproj)

ax4 = plt.subplot(2, 4, 3, projection=axproj)
ax5 = plt.subplot(2, 4, 4, projection=axproj)


ax21.pcolor(lonsok, latsok, z4)
ax21.set_title('pcolor data no crash (-175)\n plot in data projection')
mesh22 = ax22.pcolor(lonsno, latsno, -z4, cmap='cool')
ax22.set_title('pcolor data crash (-176)\n plot in data projection')


ax1.pcolor(xeok, yeok, z4)
ax1.set_title('pcolor\ndata no crash\n coords in\naxis projection')
ax2.pcolor(lonsok, latsok, z4, transform=ccrs.PlateCarree())
ax2.set_title('pcolor\ndata no crash\n coords in\ndata projection')


ax4.pcolor(xeno, yeno, -z4, cmap='cool')
ax4.set_title('pcolor\ndata crash\n coords in axis\nprojection')

args = (lonsno, latsno, -z4)
kwargs = {
        'transform' : ccrs.PlateCarree(),
        'cmap' : 'cool'}

try:
    mesh = ax5.pcolor(*args, **kwargs)
except ValueError as err:
    print('Exception ValueError HT', err)
    # geoaxes.py has been modified to store the crashed artist
    # in crashedpcolorchild 
    mesh = ax5.crashedpcolorchild 
    for child in ax5.get_children():
        if isinstance(child, matplotlib.collections.PolyCollection):
            if comparecoll(child, mesh):
                child.remove()
                print("Crashing artist removed")
    # try to plot each cell independently 
    # if it crashes the cell artist is removed from the plot
    properties = mesh.properties()
    paths = properties['paths']
    kwargs['vmin'], kwargs['vmax'] = properties['clim']
    array = properties['array']
    
    for arr, path in  zip(array, paths):
        try:
            coll = ax5.pcolor(
                    np.reshape(path.vertices[[0,1,3,2]][...,0], (2,2)),
                    np.reshape(path.vertices[[0,1,3,2]][...,1],(2,2)),
                    [[arr]], **kwargs)
        except ValueError as err:
            print( path.vertices[0:4:2][...,0],
                path.vertices[0:4:2][...,1],
                arr)
            print('Exception ValueError 2 HT', err)
            crashed = ax5.crashedpcolorchild
            for child in ax5.get_children():
                if isinstance(child, matplotlib.collections.PolyCollection):
                    if comparecoll(child, mesh):
                        child.remove()
                        print(err, "Crashing artist removed")                        


for ax in [ax21,  ax1, ax2]:
    ax.set_global()
    ax.coastlines(linewidth=0.2)

    ax.scatter(
        np.mean(np.array([lonsok[:-1, :-1], lonsok[1:, 1:]]), axis=0),
        np.mean(np.array([latsok[:-1, :-1], latsok[1:, 1:]]), axis=0),
        50, z4, edgecolors='black',
        transform=ccrs.PlateCarree())

for ax in [ax22, ax4, ax5]:
    ax.set_global()
    ax.coastlines(linewidth=0.2)
    ax.scatter(
        np.mean(np.array([lonsno[:-1, :-1], lonsno[1:, 1:]]), axis=0),
        np.mean(np.array([latsno[:-1, :-1], latsno[1:, 1:]]), axis=0),
        50, -z4, edgecolors='black', cmap='cool',
        transform=ccrs.PlateCarree())


ax5.set_title('pcolor\ndata crash\n coords in data\nprojection\n raises a ValueError')

plt.savefig("img/small_crash_pcolor_caught_onlybadcell.png",
            dpi=100, pad_inches=0)

This is the console output:

UserWarning: Unexpected ValueError in pcolor the pcolor plot artist is removed. Some cells could not be transformed in the axis projection.
  warnings.warn("Unexpected ValueError in pcolor"
...cartopy/mpl/geoaxes.py:1775: UserWarning: Unexpected ValueError in pcolor the pcolor plot artist is removed. Some cells could not be transformed in the axis projection.
  warnings.warn("Unexpected ValueError in pcolor"
Exception ValueError HT A LinearRing must have at least 3 coordinate tuples
[-176. -133.] [26. 42.] 19.5
Exception ValueError 2 HT A LinearRing must have at least 3 coordinate tuples

And the plot:

The left plot doesn't crash but the right one does. As this new code manages the exceptions, only the cell that crashes (the pink one) is not plotted while the others are properly plotted.

This code cannot be transferred as is in the PR as there are many cases to manage, cmap, facecolors, edgecolors... so it would need some lines of code and some testing. It could be done in a further improvement.

[htonchia]

Analysis of the cause of the crash.

Summary
The crashes happens when a cell is transformed from the source projection to the axis projection.
The cell is transformed in a polygon or a list of polygons.
In some case, those polygons are intersected with the projection boundary.
It crashes when the intersection of the polygon and the projection boundary is a series of two points segments instead of the expected polygon.

Images
The cell in source projection (the pink one):

One of the polygons of the pink cell after projection to the plot projection. The cell is defined as a hole:

The polygon after intersection (and inversion) with the projection boundary using shapely:

Resulting data
It produces the expected result graphically, but the data are not in the expected format.
GEOMETRYCOLLECTION (LINESTRING (12755636.1863 0, 12685759.47662721 1330784.929796209), LINESTRING (12685759.47662721 1330784.929796209, 12476894.93146266 2646989.501416921), LINESTRING (12476894.93146266 2646989.501416921, 12131330.91447088 3934193.102141946),...

Those lineString contains only two points and cannot be converted to linearRing hence the crash.

Further analysis
The intersection of the shape with the boundary happens in crs.py:

cartopy/lib/cartopy/crs.py

Line 617 in 136a047

polygon = boundary_poly.intersection(polygon)

And further done the line, in shapely geometry multipolygon, it considers that if the object is not a polygon, then the first one is the outer shell and the others are the holes. So in our case, it considers that the first linestring is the outer shell and that's what crashes further done the line.

    for i, ob in enumerate(obs):
        if isinstance(ob, polygon.Polygon):
            shell = ob.exterior
            holes = ob.interiors
        else:
            shell = ob[0]
            holes = ob[1]
        geom, ndims = polygon.geos_polygon_from_py(shell, holes)

Possible fixes
I assume shapely cannot be modified.

1. Find a way to remove the cells that are not in the right format after the intersection.
2. Find a way to change the geometry collection into a polygon and go on with the same process.
3. Process the geometry collection in a new branch, adapted to its type.

In the mean time, I proposed to go through with the PR.

[htonchia]

Test of a fix following the above analysis:

2. Find a way to change the geometry collection into a polygon and go on with the same process.

Summary

In some case the projected cell is transformed in a shapely geometryCollection that contains segments (lineString with 2 points) that crash in other shapely functions.
So the idea has been to remove the lineString from the geometryCollection.
It appeared that small bits of the cell were also transformed in the full boundary polygon, thus obscuring the plot, so the large bits (large area wise) compared to the full boundary have been removed too.

Modification to crs.py

cartopy/lib/cartopy/crs.py

Line 617 in 136a047

polygon = boundary_poly.intersection(polygon)

                    # if polygon is a GeometryCollection with lineString
                    # of 2 points, the plot will crash wih a valueError
                    # A LinearRing must have at least 3 coordinate tuples
                    # It is due to the fact, that further in the process
                    # it expects polygons or an outer shell and a list of inner
                    # holes, but not the geometryCollection
                    # To avoid a crash, this removes all the linestring
                    # from the geometryCollection.
                    # It appears that, despite it doesn't crash anymore
                    # there are small bits of polygon that once inverted and
                    # intersected with the boundary_poly will cover the
                    # full boundary_poly.
                    # Those are artefacts and should be removed.

                    if isinstance(polygon, sgeom.collection.GeometryCollection):
                        if len(polygon) > 0:
                            oblist = []
                            for obj in polygon:
                                if not isinstance(obj,
                                                  sgeom.linestring.LineString):
                                    if polygon.area < boundary_poly.area * 0.9:
                                        oblist.append(obj)
                            if len(oblist) != 1:
                                # there still may be issues if len(oblist) > 1
                                # but has not been encountered
                                polygon = sgeom.collection.GeometryCollection(
                                    oblist)
                            else:
                                polygon = oblist[0]

Result

It doesn't crash the plot anymore due to lineString in a geometryCollection.
Two sets of data have been plotted, one in viridis that doesn't crash pcolor , and a second in cool, very close to the first one, that crashes pcolor without the fix (as seen earlier above).
Below is the result with the crs.py fix.

Implications

It does modify crs in a function that may be used elsewhere.
But one thing is sure. A geometryCollection with linesString of two points will crash the code.
A second thing is sure. If the intersection with the boundary_poly is a geometryCollection, this object doesn't conform to what is expected later in the code: an outer shell and a list of holes.
So I cannot imagine a case where this geometryCollection could work as expected, but it doesn't mean it doesn't exist.

Code to reproduce the plot

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

proj = ccrs.LambertAzimuthalEqualArea

lats = np.linspace(10, 42, 3)
lonsok = np.linspace(-175, -90, 4)
lonsno = np.linspace(-176, -90, 4)
lonsok, latsok = np.meshgrid(lonsok, lats)
lonsno, latsno = np.meshgrid(lonsno, lats)

axproj = proj(central_latitude=-30)

xyeok = axproj.transform_points(ccrs.PlateCarree(), lonsok, latsok)
xyeno = axproj.transform_points(ccrs.PlateCarree(), lonsno, latsno)

xeok = xyeok[..., 0]
yeok = xyeok[..., 1]
xeno = xyeno[..., 0]
yeno = xyeno[..., 1]

np.set_printoptions(precision=0)
np.set_printoptions(precision=3)

z4 = latsok[:-1, :-1] * (1 + 0.01 * lonsok[:-1, :-1])
fig = plt.figure(figsize=(13.1, 6.4), dpi=100)
plt.tight_layout()
fig.suptitle('LambertAzimuthalEqualArea central_latitude=-30')

ax21 = plt.subplot(2, 2, 3, projection=ccrs.PlateCarree())
ax22 = plt.subplot(2, 2, 4, projection=ccrs.PlateCarree())

ax1 = plt.subplot(2, 4, 1, projection=axproj)
ax2 = plt.subplot(2, 4, 2, projection=axproj)

ax4 = plt.subplot(2, 4, 3, projection=axproj)
ax5 = plt.subplot(2, 4, 4, projection=axproj)


ax21.pcolor(lonsok, latsok, z4)
ax21.set_title('pcolor data no crash (-175)\n plot in data projection')
ax22.pcolor(lonsno, latsno, -z4, cmap='cool')
ax22.set_title('pcolor data for crash (-176)\n plot in data projection')


ax1.pcolor(xeok, yeok, z4)
ax1.set_title('pcolor data no crash\n coords in axis projection')
ax2.pcolor(lonsok, latsok, z4, transform=ccrs.PlateCarree())
ax2.set_title('pcolor data no crash\n coords in data projection')


ax4.pcolor(xeno, yeno, -z4, cmap='cool')
ax4.set_title('pcolor data for crash\n coords in axis projection')


ax5.pcolor(lonsno, latsno, -z4,
           transform=ccrs.PlateCarree(),
           cmap='cool')



for ax in [ax21, ax1, ax2]:
    ax.set_global()
    ax.coastlines(linewidth=0.2)

    ax.scatter(
        np.mean(np.array([lonsok[:-1, :-1], lonsok[1:, 1:]]), axis=0),
        np.mean(np.array([latsok[:-1, :-1], latsok[1:, 1:]]), axis=0),
        50, z4, edgecolors='black',
        transform=ccrs.PlateCarree())

for ax in [ax22, ax4, ax5]:
    ax.set_global()
    ax.coastlines(linewidth=0.2)
    ax.scatter(
        np.mean(np.array([lonsno[:-1, :-1], lonsno[1:, 1:]]), axis=0),
        np.mean(np.array([latsno[:-1, :-1], latsno[1:, 1:]]), axis=0),
        50, -z4, edgecolors='black', cmap='cool',
        transform=ccrs.PlateCarree())


ax5.set_title('pcolor data crash\n coords in data projection\n would raise a ValueError')

plt.savefig("img/small_crash_pcolor_caught_onlybadcell.png",
            dpi=100, pad_inches=0)

[htonchia]

Analysis of why crs.py _rings_to_multi_polygon issue of big cells

Summary

After fixing the valueError 'A LinearRing must have at least 3 coordinate tuples', it appeared that the result was not the expected one.
The analysis showed an issue in the handling of polygon holes outside of any shell.
To produce a polygon from a hole we need to invert the hole and intersect the result with the projection boundary. When there are several holes, those need to be joined before the inversion but crs joins them after, that's why it doesn't produce the expected result.

@dopplershift
I suggest to drop the current PR on geoaxes and replace it with a new PR on crs that fixes:

- the valueError exception 'A LinearRing must have at least 3 coordinate tuples'
- the inversion of the holes

Example on why the holes need to be joined before inversion

In the example below a polygon abcda is defined by three holes and a projection boundary.
ABCD is the box including the projection bounding box.
abcd is the polygon defined by three holes.
holes are holes if there are clock wise.
shells are ccw.

D ______________E_____C
 |              |     |
 |              |     |
 |              |     |
 |              |     |
d|______________|c    |
 |              |     |
a|______________|b    |
 |              |     |
 |              |     |
 |______________|_____|  
A               F     B

The holes defining abcd are:
dDEcd
FECBF
AabFA

The result expected is the polygon abcda.

The image below shows:

• the polygons and points
• the result with the current crs.py _rings_to_multi_polygon process
• the result with the union of the holes before the inversion

The expected result is the right plot, with the union of the holes before the inversion.

Code to reproduce the example

"""
this is to show the issue of holes inversion and union
in crs lines 590 and further
in def _rings_to_multi_polygon(self, rings, is_ccw):

In this function there is a case where several polygons
define the holes while there is no outer shell.
The polygon is thus defined as outside the holes

To define the polygon one has to make the union of holes
and invert them in the projection  bounding box
but instead of doing this in crs, they invert the holes
separately and plot each of them.
that's doesn't work when there is more than one hole.

Here is an example.
ABCD is box including the projection bounding box
abcd is the polygon defined by three holes
holes are holes if there are clock wise
shells are ccw.

D ______________E_____C
 |              |     |
 |              |     |
 |              |     |
 |              |     |
d|______________|c    |
 |              |     |
a|______________|b    |
 |              |     |
 |              |     |
 |______________|_____|
A               F     B

The holes defining abcd are:
    dDEcd
    FECBF
    AabFA

the bounding box is ABCD

In our case, we need also to intersect the result in
the boundary polygon of the projection, let's say a circle.

to plot abcd we need to:
    union the holes
    invert the union in the bounding box
    intersect the result with the boundary polygon
    plot the result


But crs does the following:
    invert each hole in the bounding box
    intersect each inverted hole with the boundary polygon
    plot each result

which isn't the right way to do it as we will see below
"""
import math
import numpy as np
import shapely.geometry as sgeom
import matplotlib.pyplot as plt

def define_pols():
    A = (0, 0)
    B = (10, 0)
    C = (10, 10)
    D = (0, 10)
    E = (7, 10)
    F = (7, 0)
    a = (0, 3)
    b = (7, 3)
    c = (7, 5)
    d = (0, 5)

    texts = ['A', 'B', 'C', 'D', 'E', 'F', 'a', 'b', 'c', 'd']
    poss = [A, B, C, D, E, F, a, b, c, d]

    hole1 = sgeom.Polygon([d, D, E, c, d])
    hole2 = sgeom.Polygon([F, E, C, B, F])
    hole3 = sgeom.Polygon([A, a, b, F, A])

    pbox = sgeom.Polygon([A, B, C, D])

    pboundary_poly = sgeom.Polygon([(5+5*math.cos(teta), 5+5*math.sin(teta))
        for teta in np.linspace(0, 2 * math.pi, 100)])

    axe = plt.subplot(1, 3, 1)
    _ = [axe.fill(*pol.exterior.xy) for pol in
         [hole1, hole2, hole3]]

    _ = [axe.plot(*pol.exterior.xy) for pol in
         [hole1, hole2, hole3, pbox, pboundary_poly]]
    axe.set_title('The polygons')

    for text, pos in zip(texts, poss):
        plt.text(*pos, text, fontsize=12)
    return [hole1, hole2, hole3], pbox, pboundary_poly

def crs_process(holes, box, boundary_poly):
    result = []
    for hole in holes:
        holeinv = box.difference(hole)
        pol = boundary_poly.intersection(holeinv)
        result.append(pol)
    axe = plt.subplot(1, 3, 2)
    axe.set_xlim((-1, 11))
    axe.set_ylim((-1, 11))
    for pol in result:
        plt.fill(*pol.exterior.xy)
    axe.set_title('crs process')

def good_process(holes, box, boundary_poly):
    poly = holes[0]
    for hole in holes:
        poly = poly.union(hole)
    holeinv = box.difference(poly)
    result = boundary_poly.intersection(holeinv)
    axe = plt.subplot(1, 3, 3)
    axe.set_xlim((-1, 11))
    axe.set_ylim((-1, 11))
    plt.fill(*result.exterior.xy)
    axe.set_title('The right process')



plt.figure('crs_good', figsize=(12, 6))
holes, box, boundary_poly = define_pols()
crs_process(holes, box, boundary_poly)
good_process(holes, box, boundary_poly)

Fixing crs.py to union the holes before inversion

The code

cartopy/lib/cartopy/crs.py

Line 590 in 136a047

if interior_rings:

should be modified according to the following, on top of the modification to handle the valueError exception:

        if interior_rings:
            boundary_poly = self.domain
            x3, y3, x4, y4 = boundary_poly.bounds
            bx = (x4 - x3) * 0.1
            by = (y4 - y3) * 0.1
            x3 -= bx
            y3 -= by
            x4 += bx
            y4 += by
            holes = sgeom.Polygon()
            # union the interior_rings before inverting them
            for ring in interior_rings:
                # Use shapely buffer in an attempt to fix invalid geometries
                polygon = sgeom.Polygon(ring).buffer(0)
                if not polygon.is_empty and polygon.is_valid:
                        holes = holes.union(polygon)
            if not holes.is_empty:       
                x1, y1, x2, y2 = holes.bounds
                bx = (x2 - x1) * 0.1
                by = (y2 - y1) * 0.1
                x1 -= bx
                y1 -= by
                x2 += bx
                y2 += by
                box = sgeom.box(min(x1, x3), min(y1, y3),
                                max(x2, x4), max(y2, y4))
    
                # Invert the union of the interior_rings
                polygon = box.difference(holes)
    
                # Intersect the inverted union of the interior_rings 
                # with the boundary
                polygon = boundary_poly.intersection(polygon)
 

I suggest to drop the current geoaxes PR for a crs PR.

Implications

This fixes a geometry union/difference error.
I cannot see any implication of doing it the right way.

[QuLogic]

It looks like you've figured out what to do here. Can you open a new PR with your proposed change for the hole calculations? (Please also give your commits and PRs more descriptive messages.)

[htonchia]

@QuLogic. I will open a new PR to clean all this up.