cartopy.mpl.feature_artist.FeatureArtist.draw does not preserve ordering for facecolors, etc

[jtbraun]

Example code: ax.add_feature(cfeature.ShapelyFeature(shapefile.geometries(), crs, facecolors=['red','green','blue'])

Intent: I want to assign a different color to each geometry in shapefile.geometries(). I can loop through shapefile.records().attr[] to create a list of facecolors (red/green/blue in the example above).

However, the cartopy.mpl.feature_artist.FeatureArtist's draw() routine culls the list of geometries via geoms = self._feature.intersecting_geometries(extent), and then creates a set of N paths for each geometry.

For this to work, draw() will need to select the correct kwargs['facecolors'] for each geometry, and create `final_kwargs['facecolors'], duplicating that color for each path resulting from the geometry.

These should be handled similarly: edgecolors, linewidths, linestyles and antialiaseds

[pelson]

Thanks @jtbraun. I confess I've wanted this myself, but have always told myself that the extra complexity wasn't worth the extra functionality. As you can see in the example at http://scitools.org.uk/cartopy/docs/latest/examples/hurricane_katrina.html colouring geometries means you (currently) have to go down a slightly different code path than you would with features.

The truth is, even if we close the gap between functionality in the way you describe, there will always be a further case which warrants further functionality. For instance, it may be that I wanted to color geometries based on their area, rather than by specifying specific colors for each geometry.

Would it be helpful if there were an example in the gallery of using the shapereader to colour a country based on some other attribute of the geometry/shape?

[has2k1]

I run into this issue while trying to create a geom_carto for plotnine.

Maybe a solution could start with adding an option for the Feature.intersecting_geometries method to
return the indices of the selected geometries.

def intersecting_geometries(self, extent, include_indices=False):
    ...

This would be backwards compatible.

Then in FeatureArtist.draw the list-like kwargs (with the same number of elements as there are geometries) to Matplotlib can be adjusted.

[has2k1]

@pelson, is the suggestion in the previous comment sensible? I can work on a PR for it.

[pelson]

@has2k1 - definitely interested in getting your hands on the code!
Is solving this simply a matter of processing callbacks around https://github.com/SciTools/cartopy/blob/master/lib/cartopy/mpl/feature_artist.py#L172?

For example:

shapes = shapefile.geometries()

def colorization(shape):
    if shapes.record.name == "USA":
        color = "Red"
    else:
        color = "Blue"
    return color

ax.add_feature(cfeature.ShapelyFeature(shapes, crs),
               facecolors=colorization)

I'm glossing over the need to include CRS, projected and unprojected geometry, and the record information (if still available at that point...), but in principle I'd rather be handling that than assuming the incoming shapes all need to be in memory at the same time (maybe they already do though...).

[has2k1]

I had not considered using callbacks. Although they would provide an improvement over the workaround in the current hurricane Katrina example, they would have a few issues.

1. They would not work for my use-case.
2. They create a special case (different from the rest of Matplotlib) for setting parameters used by artists.
3. If we consider that the reason for this issue is that when the geometries are culled it affects the user's expectation of how things should work. Using callbacks would leak the solution into the user API, instead of addressing the expectation.

However, as there are different types downloadable features, with callbacks the user does not have to get the geometries into memory. Even so, I think this would only save the user 2-3 lines of code i.e. creating a list of all the geometries and looping through them to determine any of colours, linewidths, etc. Also, as you would do the calculations for geometries that will be culled anyway, I do not know how big of a deal that would be for most cases.

[pelson]

I've just pushed a PR (#1019) to allow this functionality.
The work in progress would mean the following code would produce zoom invariant coloring (made a little more complex than necessary by using matplotlib's cycler functionality):

"""
Consistent coloring of geometries
---------------------------------

When we zoom in to a map, cartopy only draws the geometries that are
necessary. This can be problematic if you wish to assign a constant colour to
a geometry based on its position in your list of geometries.

This example demonstrates using the cartopy styler functionality to get around
this issue.


"""
import matplotlib.patches as mpatches
import matplotlib.pyplot as plt
import shapely.geometry as sgeom

import cartopy.crs as ccrs
import cartopy.io.shapereader as shpreader



def main():
    fig = plt.figure()
    ax = fig.add_axes([0, 0, 1, 1], projection=ccrs.LambertConformal())

    ax.set_extent([-125, -66.5, 20, 50], ccrs.Geodetic())

    shapename = 'admin_1_states_provinces_lakes_shp'
    states_shp = shpreader.natural_earth(resolution='110m',
                                         category='cultural', name=shapename)

    ax.set_title('Zoom invariant coloring of geometries')

    geometry_styles = {}
    from cycler import cycler

    style_cycle = cycler(facecolor=['r', 'g', 'b'])
    from itertools import cycle
    style_cycle = cycle(style_cycle)

    def deterministic_colorization(geometry):
        style = geometry_styles.setdefault(id(geometry), next(style_cycle))
        return style

    ax.add_geometries(
        shpreader.Reader(states_shp).geometries(),
        ccrs.PlateCarree(),
        styler=deterministic_colorization)

    plt.show()