wzbsocialsciencecenter / geovoronoi Goto Github PK

Hi,
I feel that a method returning a dictionary containing the point x and y data along with the respective voronoi polygon would be a nice addition, to easily generate a GeoDataFrame. Or is there a method, that i oversaw? Without further inspection of the source code my quick and dirty postprocessing approach looks like this:

# points_arr: numpy array of some points

gpd_input_dict= {
    'x': [],
    'y': [],
    'geometry': []
}

for pt, poly in region_pts.items():
    gpd_input_dict['x'] += [points_arr[poly[0]][0]]
    gpd_input_dict['y'] += [points_arr[poly[0]][1]]
    gpd_input_dict['geometry'] += [region_polys[pt]]


gdf_voronoi = gpd.GeoDataFrame(gpd_input_dict)

However, these information must be somewhere in the voronoi_regions_from_coords method as the mapping of the dictionaries is derived from it.

It might as well be numpy arrays or lists, just having something that returns data which can be directly sent to geopandas would be nice. Do you think it's worth adding some method like this?

Best regards and have a nice weekend

Hi!
I found your module while working with geolocation data. My problem is, that final map is not as precise as I would expect. So I'm wandering if using your module could help me.

My workflow is following:

1. I have some geolocation data (lat,lon) and I'm putting it inside scipy voronoi. I'm calculating vertices.
2. I then calculate distance between vertices and points using geopy.distance.vincenty(coords_1, coords_2).km
3. I'm selecting biggest circle without points based on above calculations and I'm plotting them using folium.

Here's my result.

This circle should be biggest possible circle in Poland without marked points. I believe it should touch three points. There are to possible sources of this error:

1. I made a mistake somewhere.
2. One of modules calculates something not precisely enough. I know that scipy.spatial is designed to work with 2d-flat data, and folium and geopy should work nicely with non flat earth data.

So here is my question: If I would replace default scipy voronoi with your geovoronoi, would it calculate voronoi cells correctly (on earth surface) based on latitude and longitude? Did you used any geographical distance formula for voronoi cells?

If you are curious,
my code is here: https://nbviewer.jupyter.org/gist/bbcf796d75dfe7a8d06040fcd73d35eb

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

In the function, "voronoi_regions_from_coords" on line 72 when the function "assign_points_to_voronoi_polygons" is called, the argument "return_unassigned_points" is set to true instead of the variable "return_unassigned_points".

After far too much effort and being stupid, I finally have a Voronoi map. Nice! I do have a question though: is it possible to have this map drawn in folium, say, with a street map drawn on it? Then I'd be able (I hope) to replace the data point labels with folium popups or tooltips. And show the map in its geographic location.

I already have a nice folium map which shows the region and its data points. So what I want to with it is put a Voronoi map over it (clearly then the colors would have to semi-transparent). Is this possible?

Thanks!

Hi,

I have a set of points over the Netherlands and also the polygon of the country, both in .shp format. Originally, both the points and the polygon were encoded in EPSG:4326 (WGS84), but they were also reprojected to EPSG:3395 (Web Mercator). For the record, the polygon is actually a MultiPolygon type.

So I ran the following code:

from shapely.ops import cascaded_union
from geovoronoi.plotting import subplot_for_map, plot_voronoi_polys_with_points_in_area
from geovoronoi import voronoi_regions_from_coords, points_to_coords

# Read shapefile 
gdf = gpd.read_file(path_loc) 
print(gdf.head())
print("Dataframe CRS:", gdf.crs)

boundary = gpd.read_file(path_nl)
gdf_proj = gdf.to_crs(boundary.crs)
print("Boundary CRS: ", boundary.crs)

boundary_shape = cascaded_union(boundary.geometry)
coords = points_to_coords(gdf_proj.geometry)

# Calculate Voronoi Regions
poly_shapes, pts, poly_to_pt_assignments = voronoi_regions_from_coords(coords, boundary_shape)

But unfortunately the program ends with the following error:

  station_id  longitude  ...        Y_webm                        geometry
0        279   6.572446  ...  6.902455e+06  POINT (731641.307 6902454.683)
1        225   4.554683  ...  6.850046e+06  POINT (507024.938 6850045.637)
2        235   4.780770  ...  6.935182e+06  POINT (532192.902 6935182.126)
3        240   4.789374  ...  6.823666e+06  POINT (533150.640 6823666.308)
4        242   4.920471  ...  6.993052e+06  POINT (547744.338 6993052.456)

[5 rows x 9 columns]

Dataframe CRS: {'init': 'epsg:3395'}
Boundary CRS:  {'init': 'epsg:3395'}

Traceback (most recent call last):
  File "[MYPATH]/eda/00_make_voronoi_diagram.py", line 41, in <module>
    poly_shapes, pts, poly_to_pt_assignments = voronoi_regions_from_coords(coords, boundary_shape)
  File "[MYPATH]/miniconda3/envs/defpy37/lib/python3.7/site-packages/geovoronoi/_voronoi.py", line 202, in voronoi_regions_from_coords
    geom_polys, geom_pts = region_polygons_from_voronoi(vor, geom, return_point_assignments=True, **kwargs)
  File "[MYPATH]/miniconda3/envs/defpy37/lib/python3.7/site-packages/geovoronoi/_voronoi.py", line 364, in region_polygons_from_voronoi
    raise RuntimeError('ridge line must intersect with surrounding geometry from `geom`')
RuntimeError: ridge line must intersect with surrounding geometry from `geom`

I also tried the same procedure, but using WGS84 first, then EPSG:28992 (official CRS of the Netherlands) for both datasets, but getting the same error. Since I am not sure of the source of the problem, I thought I could open this issue to see if I can get some more info on the error and perhaps some hints to fix this problem.

Thanks a lot for your help! :-)

I liked to calculate the areas of the voronoi regions,
but in my example they was much smaller than they should be (0.11715301m² for whole of Berlin).

When I tried do fix the bug, I took a look at random_points_and_area.py
and printet the sum of the areat there

area_sum=0
for a in poly_areas:
   print("area is: %s"%a)
    area_sum=area_sum+a
print("area_sum is: %s"%area_sum)
print("area_sum from shape is: %s"%calculate_polygon_areas([area_shape], m2_to_km2=True))

where I got:

area is: 38047.414248663496

area is: 81558.23012468345

area is: 45652.50539830615

area is: 37314.4644019487

area is: 85299.03205647763

area is: 97225.72636558591

area is: 60036.51052754671

area is: 124761.16667811942

area is: 128687.40068640596

area is: 28103.058823670108

area is: 68784.31109230645

area is: 36575.84996001104

area is: 16695.936399656446

area is: 15257.965532032309

area_sum is: 863999.5722954138

area_sum from shape is: [863999.57229541]

But according to Wikipedia it should be: 505.970 km²

I noticed that part after computing the polygons the code assigns all points to these polygons. This results in an O(N^2) implementation which is untenable (memory and speed) when there are a lot of points (I have 200000). Assigning most points to cells by computing the matrix of distances is unnecessary. Most points (90-95%) are unambiguously assigned by Voronoi from scipy.spatial. It would be fantastic if the library (which is definitely filling a much needed niche) could be sped up by using these unambiguous assignments to cells that don't have edges that extend to infinity.

Thanks!

Hi!
I am trying to run the example that you placed on the first page. But a got an error that is the same for other examples that I have found all over the internet. It will happen when I call plot_voronoi_polys_with_points_in_area(ax, area_shape, region_polys, coords, region_pts) The error is as follows:

{ "name": "IndexError", "message": "too many indices for array: array is 0-dimensional, but 2 were indexed", }
For reproducing this error you should just run your example.

I thought I'd finally found the holy grail - a voronoi tool in python. Unfortunately, the results come out all janky. Below is the data and code to recreate this unlovable monster. You can compare it against the QGIS voronoi tool.

points.zip

import geovoronoi as gv
from shapely.geometry import box
from pathlib import Path
import pandas as pd, numpy as np
import geopandas as gpd,
g = 'geometry'

pts = gpd.read_file('points.gpkg')
ptz = pts[g].to_list()
crds = gv.points_to_coords(ptz)

bbox = coordz.dissolve().bounds
bbox = box(*bbox.iloc[0].to_list())
bbox
polys,_ = gv.voronoi_regions_from_coords(crds,bbox,per_geom=False)
vor = gpd.GeoSeries(polys,crs=pts.crs)
vor
gdf= gpd.GeoDataFrame(geometry=vor,crs=vor.crs)
gdf
gdf.to_file('uglyVoronoi.gpkg')

    I used the method: Voronoi_regions_from_coords to form Voronoi polygons in geographical regions, and then calculated the area of each Voronoi polygon, but the result was different from ArcGIS calculation。My code is as follows：

polys, pts = voronoi_regions_from_coords(v_coords, boundary_shape)
boundary_area = Decimal(boundary_shape.area)

for index,key in enumerate(polys):
    area = Decimal(polys[key].area)
    weigt = area / boundary_area


 Looking forward to your reply！   

I have been using geovoronoi for a while.
Recently I had to make voronoi polygons for 30K or more points and geovoronoi got way too slow for my purposes.
I tried PyGEOS, which has a voronoi_polygons function that is actually faster than scipy's, but there is no way that I found to refer back to the original points without doing a spatial join.
So I decided to try my hand at using scipy's scipy.spatial.Voronoi myself (like geovoronoi).
I quickly ran into the same place where geovoronoi gets all its runtime and complexity - dealing with infinite ridges/polygons.
Now, I had a big breakthrough.
I just add four more points to the voronoi far away - this way the infinte polygons are all oustide the bounding geometry.
Then every polygon I care about is finite and I just ignore the last 4.
I get the polygons now in 300ms vs 5 seconds using geovoronoi.

I'll clean up my code and put it here, but the idea is simple.

I receive the following warning often when using geovoronoi.voronoi_regions_from_coords():

"The 'cascaded_union()' function is deprecated. Use 'unary_union()' instead."

Looks like there might be just one place in the actual code that uses cascaded_union():

...and a couple of places in examples and tests...

I'll see about submitting a PR to update.

I am running python3.7 with everything installed and when running the examples I get the following error:

Traceback (most recent call last):
File "duplicate_points.py", line 111, in
points_markersize=np.array(count_per_pt)*10)
File "/Users/dianacrisan/Library/Python/3.7/lib/python/site-packages/geovoronoi/plotting.py", line 137, in plot_voronoi_polys_with_points_in_area
_plot_polygon_collection_with_color(ax, [area_shape], color=area_color, edgecolor=area_edgecolor, **plot_area_opts)
File "/Users/dianacrisan/Library/Python/3.7/lib/python/site-packages/geovoronoi/plotting.py", line 204, in _plot_polygon_collection_with_color
geoms, _ = _flatten_multi_geoms(geoms)
File "/Users/dianacrisan/Library/Python/3.7/lib/python/site-packages/geopandas/plotting.py", line 30, in _flatten_multi_geoms
if not geoms.geom_type.str.startswith('Multi').any():
AttributeError: 'list' object has no attribute 'geom_type'

Hi,

We've been trying to switch to Shapely's implementation as per your suggestion but are having trouble. Our previous example looked something like this:

# Old approach
# Now get coords in the correct format
coords = points_to_coords(s_gpd.geometry)
region_polys, region_pts = voronoi_regions_from_coords(coords, bounding_shape)

# Reconstruct a dataframe with points and polys
s_gpd["site_voronoi"] = pd.Series()

# Loop through the region polys and assign them to the correct index of points
for key in region_polys:
    idx = region_pts[key][0]
    s_gpd.loc[s_gpd.index[idx], "site_voronoi"] = region_polys[key]

and we think the new approach would just be:

# New approach with Shapely
s_gpd["new_site_voronoi"] = voronoi_polygons(s_gpd.geometry, extend_to=bounding_shape)

But this seems to give a series of multipoint geometry with the bounding shape completely filled. Do you have any suggestions as to why that might be?

I'm trying to generate geovoronoi regions using a GeoDataFrame. The input points come from the dataframe's geometry. Once I have the resulting polygons, I need a mapping of the rows in the dataframe to the regions generated by the library. I don't find any way to do this using the outputs from voronoi_regions_from_coords. Here follows an example program using cities in South Africa. The idea would be to get the right city name back for each Voronoi region. It's not correct at the moment. What am I missing?

import logging

import matplotlib.pyplot as plt
import geopandas as gpd

from geovoronoi import coords_to_points, voronoi_regions_from_coords
from geovoronoi.plotting import subplot_for_map, plot_voronoi_polys_with_points_in_area

logging.basicConfig(level=logging.INFO)
geovoronoi_log = logging.getLogger('geovoronoi')
geovoronoi_log.setLevel(logging.INFO)
geovoronoi_log.propagate = True

mercator = 3395
wsg84 = 4326
hartebeesthoek94 = 4148 # Frequently used for southern Africa
crs = hartebeesthoek94

country = 'South Africa'

print(f'loading country {country} from naturalearth_lowres')
world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))

area = world[world.name == country]
assert(len(area) == 1)

print('CRS:', area.crs)

area = area.to_crs(epsg=crs)
area_shape = area.iloc[0].geometry

print('loading local cities from naturalearth_cities' )
cities = gpd.read_file(gpd.datasets.get_path('naturalearth_cities'))
cities = cities.to_crs(epsg=crs)
local_cities = cities.loc[cities['geometry'].within(area_shape)]
assert(len(local_cities) >= 1) # This doesn't work for all countries!

# Calculate the `per_geom=False` voronoi diagram for the country, using the GeoDataFrame's `geometry` column
# to provide the coordinates.
voronoi_regions, voronoi_points = voronoi_regions_from_coords(local_cities.geometry, area_shape, per_geom=False)

fig, ax = subplot_for_map()

plot_voronoi_polys_with_points_in_area(ax, area_shape, voronoi_regions, local_cities.geometry, voronoi_points)

ax.set_title(f'Cities and their Voronoi regions in {country} (with per_geom=False)')

plt.tight_layout()
plt.savefig('national_city_regions.png')
plt.show()

city = 0
assert(city < len(local_cities)-1) # No point looking for cities with higher indexes

print('Goal is to map the original GeoPandas df to the voronoi_regions output')
print('\nvoronoi_points showing dict index:')
pprint(voronoi_points)
print('\nvoronoi_regions showing dict index:')
pprint(voronoi_regions)

print('\nlocal_cities showing index after the filter out of the global cities originally loaded:')
print(local_cities)

print(f'\nlocal_cities.iloc[city]["name"]: {local_cities.iloc[city]["name"]}')

print(f'\nvoronoi_points[city]: {voronoi_points[city]}')
print(f'\nvoronoi_points[city][0]: {voronoi_points[city][0]}')

print(f'\nvoronoi_regions[voronoi_points[city][0]]:')
voronoi_regions[voronoi_points[city][0]]

The region polygon in voronoi_regions[voronoi_points][0]] is Cape Town, not Bloemfontein.

Hello,

I would like to use the voronoi tesselation to aggregate points. I have got a data set with routing points of the baltic and north sea, so i'm using a shape file for the bounding box and read the points which are stored in a CSV file as WGS84 points. After that I try to generate the voronoi tesselation but it breaks with the error: generated polygon is not valid or is empty after intersection with the surrounding geometry 'geom'. How can I debug this error or fix it?
This is the code:

import pandas as pd
import geopandas as gp
import geovoronoi as gv
import geovoronoi.plotting as gvp

boundary = gp.read_file("northbalticsea.shx").geometry[0]
df = pd.read_csv('data.csv.gz', parse_dates=[0, ])

df = gp.GeoDataFrame(df, crs="EPSG:4326", geometry=gp.points_from_xy(df.LONGITUDE, df.LATITUDE))
region_polys, region_pts = gv.voronoi_regions_from_coords(df["geometry"], boundary)

fig, ax = gvp.subplot_for_map()
gvp.plot_voronoi_polys_with_points_in_area(ax, boundary, region_polys, df["geometry"], region_pts)
plt.show()

In my case my data.csv.gz contains routes in the sea, so it could be possible, that some regions contains no data points

Hi

I'm a naive geo user rather than a geo-expert and whilst trying to use outputs from the (really useful - thanks:-) geovoronoi package to plot choropleth maps using the folium package, I found the geo-json created from a geopandas representation of the data was invalid.

from geovoronoi import voronoi_regions_from_coords, points_to_coords

# convert the pandas Series of Point objects to NumPy array of coordinates
coords = points_to_coords(gdf.Coordinates)

# Calculate the Voronoi regions for the identifed points and then 
# clip within the geographic area shape
poly_shapes, pts, poly_to_pt_assignments = voronoi_regions_from_coords(coords, boundary_shape)

gdf['voronoi_boundary'] = poly_shapes
gdf = gdf.set_geometry('voronoi_boundary')

However, I could fix it with the following:

import json
from geojson_rewind import rewind

gdf = gdf.set_geometry('voronoi_boundary')
rewoundGeoJSON = rewind( json.dumps(gdf['voronoi_boundary'].__geo_interface__) )

The rewoundGeoJSON was then well formed and I could use it as the boundaryfile to generate a folium.Choropleth() map.

In the workflow, I'm not sure where the co-ordinate winding is most appropriately handled, or what the side effects of various windings are?