9

I am not sure if this is a shapely bug or a problem with my polygons.

I have a bunch of polygons that look like this: enter image description here They are level-12 watersheds provided by HydroBasins, and are presumably non-overlapping.

I do some analysis to select a group of these polygons based on their basin ID, then join that group together to make a single, large polygon: enter image description here

I have done this in four other areas without problem. However, for this polygon, I find that it actually has tiny holes, but in unusual locations: enter image description here

Here's a zoom-view: enter image description here

In this view, the yellow lines represent the border of the large polygon: enter image description here

I should note that when displaying in QGIS, no matter how far I zoom in, the smaller polygons do not have gaps between them.

What I find unusual is that these tiny gaps seem to occur at roughly constant lines of latitude. The region is Arctic (i.e. high-latitudes), so I am wondering if they appear due to the WGS84 distortions at higher latitudes.

I am using geopandas to manage and select the smaller polygons, and shapely to join them together:

polys = [] for i in idcs: polys.append(gdf.iloc[i].geometry) polyout = shapely.ops.unary_union(polys)

where "gdf" is a pandas geodataframe and "idcs" are the indices of smaller polygons I want to join together.

How can I tell if this is a shapely bug or an issue with my smaller polygons, and either way, what is the best approach to remove these? I could buffer each of the individual polygons by a very small amount before union-ing, but I would like to avoid that if possible.

Improve this question
3
  • Are the "holes" present at original borders of the polygons? Because in that case it might be that the polygons don't fit perfectly one next to another. Commented Mar 27, 2018 at 18:27
  • @joris Nope, they're not there as far as I can tell via zooming in as far as possible in QGIS. I added this to the main body, thanks. Commented Mar 27, 2018 at 18:33
  • But they are still at the location of the original borders. So even if it is not visisble, it probably is some floating point difference between the original borders. As you say, you can solve this with a tiny buffer, or you could also construct a new Polygon based on the exterior of the union like here: gis.stackexchange.com/a/262049/9828 Commented Mar 27, 2018 at 19:00

2 Answers 2

Reset to default
16

This is typically a result of the the borders not fitting perfectly one next to another (and this is very easy to get with floating point coordinates).

As an example, I use the world dataset available in geopandas, and take the unary union of the Africa continent:

import geopandas gdf = geopandas.read_file(geopandas.datasets.get_path('naturalearth_lowres')) gdf_africa = gdf[gdf['continent'] == 'Africa'] africa = gdf_africa.geometry.unary_union # only take main continent, disregard Madagascar africa = list(africa)[1]

You can see on the following plot that also this union is not perfect:

enter image description here

This is a typical problem to encounter when doing unary unions, so now how to solve this:

  1. One possible solution is to take the exterior of this polygon (so disregarding the interior holes) (taken from this answer), and create a new polygon of it:

    from shapely.geometry import Polygon africa2 = Polygon(africa.exterior)

    which gives this result:

    enter image description here

    which is already much better. There is still a small distortion at the left, because such a 'hole' touched the overall border.

  2. Another solution is to apply a very small buffer to the polygon, and undo it again (taken from this answer)

    from shapely.geometry import JOIN_STYLE eps = 0.00001 africa3 = africa.buffer(eps, 1, join_style=JOIN_STYLE.mitre).buffer(-eps, 1, join_style=JOIN_STYLE.mitre)

    which gives:

    enter image description here

    In this case, this solution looks better from the figure (no remaining distortion on the left), but this method can slightly change the coordinates (although this is often negligible).

Which of the methods works best for you, probably depends on the actual case.

Improve this answer
3
  • Excellent, thanks! I wasn't aware of the exterior method. Commented Mar 27, 2018 at 19:36
  • I would also like to state that I find it interesting that in my case, the holes seem to follow lines of latitude. Commented Mar 27, 2018 at 19:58
  • Not directly an idea about that one. Commented Mar 27, 2018 at 20:40
1

For readers in 2025+, shapely + geopandas now provide a grid_size parameter which fixes exactly this.

union(a, b, grid_size=None, **kwargs)

Merge geometries into one.

If grid_size is nonzero, input coordinates will be snapped to a precision grid of that size and resulting coordinates will be snapped to that same grid. If 0, this operation will use double precision coordinates. If None, the highest precision of the inputs will be used, which may be previously set using set_precision. Note: returned geometry does not have precision set unless specified previously by set_precision.

Cf: https://shapely.readthedocs.io/en/latest/reference/shapely.union.html

Improve this answer

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.