WZBSocialScienceCenter / geovoronoi

a package to create and plot Voronoi regions within geographic boundaries
https://datascience.blog.wzb.eu/2018/03/16/creating-and-plotting-voronoi-regions-for-geographic-data-with-geovoronoi/
Apache License 2.0
134 stars 27 forks source link

AttributeError: 'MultiPolygon' object has no attribute 'exterior' #7

Closed AtelierLibre closed 4 years ago

AtelierLibre commented 4 years ago

Hi Markus,

Thanks for sharing geovoronoi. I'm using it in a workflow to split polygons following k-means clustering similar to this post.

Occasionally however I run into this error:

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-3-f68e3fdefed9> in <module>
----> 1 geovoronoi.voronoi_regions_from_coords(centroids, polygon)

~/GitHub/geovoronoi/geovoronoi/_voronoi.py in voronoi_regions_from_coords(coords, geo_shape, shapes_from_diff_with_min_area, accept_n_coord_duplicates, return_unassigned_points)
     60 
     61     logger.info('generating Voronoi polygon lines')
---> 62     poly_lines = polygon_lines_from_voronoi(vor, geo_shape)
     63 
     64     logger.info('generating Voronoi polygon shapes')

~/GitHub/geovoronoi/geovoronoi/_voronoi.py in polygon_lines_from_voronoi(vor, geo_shape, return_only_poly_lines)
    147 
    148     # now add the lines that make up `far_points_hull` to the final `poly_lines` list
--> 149     far_points_hull_coords = far_points_hull.exterior.coords
    150     for i, pt in list(enumerate(far_points_hull_coords))[1:]:
    151         poly_lines.append(LineString((far_points_hull_coords[i-1], pt)))

AttributeError: 'MultiPolygon' object has no attribute 'exterior'

I can consistently reproduce it with quite a simple example:

centroids = np.array([[537300, 213400], [538700, 213700], [536100, 213400]])
polygon = Polygon([[540000, 214100], [535500, 213700], [535500, 213000], [539000, 213200]])
geovoronoi.voronoi_regions_from_coords(centroids, polygon)

I appreciate you saying you don't have a lot of time to work on the code but if you could give me any pointers that would be much appreciated. It only occurs on about 2 out of 180 shapes in the same area so I'd love to get over this last hurdle.

Many thanks,

Nick Screenshot from 2020-02-21 10-53-39

internaut commented 4 years ago

Thanks, this is an interesting case.

I investigated a bit and the problem is that the voronoi regions that are created do not even lie within the shape you specified (polygon). Of course this shouldn't happen. The problem lies with how the "loose ridges" of the voronoi regions are handled. These are infinite lines and in polygon_lines_from_voronoi() I try to restrict them to a certain area by calculating finite "far points". Turns out for some configurations (like yours) these far points were not "far away enough". I changed the way they're calculated and this should fix the problem.

I pushed to changes to branch issue_7. Could you please check them out and see if it works for you? If so, I'll create a new version and publish on PyPI.

AtelierLibre commented 4 years ago

Thanks for looking at this, yes, I will check out the new branch, I have a couple of other examples I can try it on.

Just in case its useful for you to look at the code, I think the same issue is present in the scipy voronoi implementation (at least in the visualisation using voronoi_plot_2d) whereas PySAL's implementation

libpysal.cg.voronoi_frames(points,` radius=None, `clip='extent')

handled it as I would have expected.

internaut commented 4 years ago

Thanks for the pointer to PySAL. Looks like I implemented almost the same as they did (they also use a default radius multiplier of 2 for calculating the far points), so it should work, but please try it out first.

AtelierLibre commented 4 years ago

That's great, the new version does seem to have fixed that case however I've been running it against the other cases that I have and I'm getting a new error.

centroids = np.array([[496712, 232672], [497987, 235942], [496425, 230252], [497482, 234933],
                   [499331, 238351], [496081, 231033], [497090, 233846], [496755, 231645],
                   [498604, 237018]])
polygon = Polygon([[495555,230875], [496938,235438], [499405,239403], [499676,239474],
                 [499733,237877], [498863,237792], [499120,237335], [498321,235010],
                 [497295,233185], [497237,231359], [496696,229620], [495982,230047],
                 [496154,230347], [496154,230347], [495555,230875]])

poly_shapes, pts, poly_to_pt_assignments = geovoronoi.voronoi_regions_from_coords(centroids, polygon)

fig, ax = subplot_for_map()
plot_voronoi_polys_with_points_in_area(ax, polygon, poly_shapes, centroids, poly_to_pt_assignments)
plt.show()

The error that I'm getting now is

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
~/anaconda3/envs/geovoronoi/lib/python3.7/site-packages/matplotlib/axes/_axes.py in _parse_scatter_color_args(c, edgecolors, kwargs, xshape, yshape, get_next_color_func)
   4289                     valid_shape = False
-> 4290                     raise ValueError
   4291             except ValueError:

ValueError: 

During handling of the above exception, another exception occurred:

ValueError                                Traceback (most recent call last)
<ipython-input-2-8e674a60519e> in <module>
     10 
     11 fig, ax = subplot_for_map()
---> 12 plot_voronoi_polys_with_points_in_area(ax, polygon, poly_shapes, centroids, poly_to_pt_assignments)
     13 plt.show()

~/GitHub/geovoronoi/geovoronoi/plotting.py in plot_voronoi_polys_with_points_in_area(ax, area_shape, poly_shapes, points, poly_to_pt_assignments, area_color, area_edgecolor, voronoi_and_points_cmap, voronoi_color, voronoi_edgecolor, points_color, points_markersize, points_marker, voronoi_labels, voronoi_label_fontsize, voronoi_label_color, point_labels, point_label_fontsize, point_label_color, plot_area_opts, plot_voronoi_opts, plot_points_opts)
    152     plot_points(ax, points, points_markersize, points_marker, color=points_color,
    153                 labels=point_labels, label_fontsize=point_label_fontsize, label_color=point_label_color,
--> 154                 **plot_points_opts)
    155 
    156 

~/GitHub/geovoronoi/geovoronoi/plotting.py in plot_points(ax, points, markersize, marker, color, labels, label_fontsize, label_color, label_draw_duplicates, **kwargs)
     95         coords = points
     96 
---> 97     ax.scatter(coords[:, 0], coords[:, 1], s=markersize, marker=marker, color=color, **kwargs)
     98 
     99     if labels:

~/anaconda3/envs/geovoronoi/lib/python3.7/site-packages/matplotlib/__init__.py in inner(ax, data, *args, **kwargs)
   1597     def inner(ax, *args, data=None, **kwargs):
   1598         if data is None:
-> 1599             return func(ax, *map(sanitize_sequence, args), **kwargs)
   1600 
   1601         bound = new_sig.bind(ax, *args, **kwargs)

~/anaconda3/envs/geovoronoi/lib/python3.7/site-packages/matplotlib/axes/_axes.py in scatter(self, x, y, s, c, marker, cmap, norm, vmin, vmax, alpha, linewidths, verts, edgecolors, plotnonfinite, **kwargs)
   4451             self._parse_scatter_color_args(
   4452                 c, edgecolors, kwargs, xshape, yshape,
-> 4453                 get_next_color_func=self._get_patches_for_fill.get_next_color)
   4454 
   4455         if plotnonfinite and colors is None:

~/anaconda3/envs/geovoronoi/lib/python3.7/site-packages/matplotlib/axes/_axes.py in _parse_scatter_color_args(c, edgecolors, kwargs, xshape, yshape, get_next_color_func)
   4295                         "acceptable for use with 'x' with size {xs}, "
   4296                         "'y' with size {ys}."
-> 4297                             .format(nc=n_elem, xs=xsize, ys=ysize)
   4298                     )
   4299                 else:

ValueError: 'c' argument has 7 elements, which is not acceptable for use with 'x' with size 9, 'y' with size 9.

I was getting a slightly different error message before with what I thought was the same example. Unfortunately I can't reproduce it now but I'm including it here incase it's useful.

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-258-556d4f47d3dc> in <module>
----> 1 poly_shapes, pts, poly_to_pt_assignments = geovoronoi.voronoi_regions_from_coords(cluster_centroids_, poly_)
      2 
      3 fig, ax = subplot_for_map()
      4 plot_voronoi_polys_with_points_in_area(ax, geo_shape, poly_shapes, coords, poly_to_pt_assignments)
      5 plt.show()

~/GitHub/geovoronoi/geovoronoi/_voronoi.py in voronoi_regions_from_coords(coords, geo_shape, shapes_from_diff_with_min_area, accept_n_coord_duplicates, return_unassigned_points, farpoints_max_extend_factor)
     77                                                                                accept_n_coord_duplicates=accept_n_coord_duplicates,
     78                                                                                return_unassigned_points=True,
---> 79                                                                                coords=coords)
     80 
     81     if return_unassigned_points:

~/GitHub/geovoronoi/geovoronoi/_voronoi.py in assign_points_to_voronoi_polygons(points, poly_shapes, accept_n_coord_duplicates, return_unassigned_points, coords)
    266             pt = points[i_pt]
    267 
--> 268             if pt.intersects(vor_poly):      # if this point is inside this Voronoi region, assign it to the region
    269                 assigned_pts.append(i_pt)
    270                 if n_assigned >= accept_n_coord_duplicates - n_assigned_dupl:

~/anaconda3/envs/geovoronoi/lib/python3.7/site-packages/shapely/geometry/base.py in intersects(self, other)
    709     def intersects(self, other):
    710         """Returns True if geometries intersect, else False"""
--> 711         return bool(self.impl['intersects'](self, other))
    712 
    713     def overlaps(self, other):

~/anaconda3/envs/geovoronoi/lib/python3.7/site-packages/shapely/predicates.py in __call__(self, this, other, *args)
     10 
     11     def __call__(self, this, other, *args):
---> 12         self._validate(this)
     13         self._validate(other, stop_prepared=True)
     14         try:

~/anaconda3/envs/geovoronoi/lib/python3.7/site-packages/shapely/topology.py in _validate(self, ob, stop_prepared)
     15 
     16     def _validate(self, ob, stop_prepared=False):
---> 17         if ob is None or ob._geom is None:
     18             raise ValueError("Null geometry supports no operations")
     19         if stop_prepared and not hasattr(ob, 'type'):

~/anaconda3/envs/geovoronoi/lib/python3.7/site-packages/shapely/geometry/proxy.py in _geom(self)
     31         if gtag != self._gtag or self._is_empty:
     32             self.empty()
---> 33             if len(self.context) > 0:
     34                 self.__geom__, n = self.factory(self.context)
     35         self._gtag = gtag

TypeError: object of type 'Point' has no len()

Happy to try it out again.

internaut commented 4 years ago

I pushed another update to issue_7, pls try it out.

It's basically the same problem as before (far points not placed "far away enough"); changed calc. again a little bit to place them farther away but I'm hesitating to use the same calc. as PySAL because I think their approach can quickly lead to a numerical overflow. Anyway, if the problem persists with another polygon/point configuration, try to increase farpoints_max_extend_factor when calling voronoi_regions_from_coords(). Default is 10.

Also, I added another validation so that an exception is raised when one point is assigned to two voronoi regions. This happened in your problem above and only later caused problems when trying to plot.

AtelierLibre commented 4 years ago

That's great, I have run through all of the examples that I have both for voronoi_regions_from_coords and for plot_voronoi_polys_with_points_in_area and they all work out. I have been trying them out in PySAL in parallel and there I found that I had to increase the voronoi_radius argument for some of the cases which I think is exactly the same 'not far away enough' issue.

Perhaps slightly off topic (and I'm not raising this as an issue) but do you know why the very simple case of a line bisecting two points is not implemented in any (?) of the libraries? As far as I understand, it is strictly a voronoi diagram. I hacked together a version myself but it would have been nice if it was built in.

Thanks again, I think this can be closed now.

internaut commented 4 years ago

Thanks for you contributions, I integrated the fix into the new version v0.2.0 that I just released.