Confusion with geopandas dataframe as Polygons input

[msarahan]

I have a GDF that I've created by overlaying a few other GDF's (state district data for various offices). I want to plot that on top of tile data.

Plotting the GDF with matplotlib via geopandas seems to work OK:

but creating a Polygons object seems to go wrong:

My ultimate goal is a colorized map with a hover tool: each district's color represents a unique combination of districts that cover that area, and the hover tool will tell you which districts cover a given area. In Travis County, there are 58 such combinations, so a legend is not really feasible. Color-coding would still be nice, though, to make the separations more clear.

I'm happy to supply the GDF - please suggest the best way to do so.

[philippjfr]

That looks pretty bad. One issue that may be causing this is that bokeh does not yet have support for polygon holes, so if there is any negative space inside the polygons you're going to get overlap. We are currently working on a project that explicitly requires this support, so I'm hoping to start working on that soon.

I'm happy to supply the GDF - please suggest the best way to do so.

Just in case it's not that, could you email it to me?

[rsignell-usgs]

https://github.com/bokeh/bokeh/pull/8340

[philippjfr]

Thanks @rsignell-usgs, that PR was actually sponsored by our ongoing project and I was going to start implementing this in holoviews/geoviews today.

[philippjfr]

@msarahan Finally starting to make some progress:

[philippjfr]

@msarahan I am encountering some weirdness which is causing various issues when projecting the data, e.g. this polygon seems very strange and when projected ends up hiding other polygons behind it:

[philippjfr]

The exact issue can be reproduced with:

wkt = 'POLYGON ((-97.6453334387044833 30.4773161392542349, -97.6453334387044691 30.4773161392542455, -97.6453334387044549 30.4773161392542491, -97.6448419999999970 30.4775050000000007, -97.6453334387044833 30.4773161392542349))'
poly = shapely.wkt.loads(wkt)
poly

ccrs.GOOGLE_MERCATOR.project_geometry(poly, ccrs.PlateCarree())

[philippjfr]

Dropping all polygons with an area < 1e-15 seems to fix things, suspect it's a floating point precision issue in the GEOS projection code.

[philippjfr]

Once merged this will work:

%%opts Polygons [width=800 height=600 tools=['hover', 'tap'] show_bounds=False] (cmap='gist_ncar' alpha=0.5)
intersection = intersection.reset_index(drop=True)
intersection_p = gv.Polygons(intersection, vdims=['index', 'Congress_CD115FP', 'State_house_NAME', 'State_senate_NAME']).options(alpha=0.5)
gv.tile_sources.Wikipedia * intersection_p

[msarahan]

Indeed, by plotting each geometry, it was clear by inspection that some were just tiny slivers where the boundaries didn't exactly match up. Floating point error is a great explanation. Hooray for crummy data.

[philippjfr]

Ah that makes sense, that makes two floating point related issues then. When computing the boundaries floating point differences result in tiny slivers being generated and when those slivers are projected onto the map they somehow explode and become larger than the entire globe.

[jbednar]

Should we have an optional utility that someone could call to match up slightly mismatched boundaries and eliminate tiny slivers, using a heuristic based on the total area enclosed (to determine how big a "sliver" is)?

[philippjfr]

Once my PR is merged we will automatically drop slivers with an area <1e-15, which seemingly cannot be accurately projected by GEOS.

[philippjfr]

Should we have an optional utility that someone could call to match up slightly mismatched boundaries

Also, this kind of thing is definitely out-of-scope for geoviews. It's something that should be handled at the shapely or even GEOS level.

[jbednar]

I've never gotten very far with opinions about what should be handled by some lower-level library I don't control. :-) In practice if it happens for users at a higher level, it has to be dealt with somehow, which usually means hacks at the very highest level (user scripts)...

[philippjfr]

In practice if it happens for users at a higher level, it has to be dealt with somehow, which usually means hacks at the very highest level (user scripts)...

These slivers are generated by applying intersect/overlap operations on shapely objects, which is not the responsibility of geoviews to deal with. As I said, we automatically eliminate these slivers, but fixing operations which were not performed using geoviews will never be our responsibility.