Points in which polygon

Question

I have several (order of M) points [(x0,y0),...,(xn,yn)]. I also have limited number of hexagons. I want to find each point falls in which hexagon. Using shapely this can be done for 1 point at a time so the loop below does the job. But is there any other way to do it faster?

from shapely.geometry import Point
from shapely.geometry.polygon import Polygon
zones = np.zeros(n)-1
for j,p in enumerate(points):
    point = Point(p)
    for i in range(len(poly_hex)): #poly_hex = list of hexagonal polygones
        polygon = Polygon(poly_hex[i])
        if polygon.contains(point):
            zones[j] = int(i)
            break

Is this something that will be done more than once (and often)? — martineau, Nov 22 '19 at 18:28
You can run multiple processes in parallel using the ```multiprocessing``` library, as checking if a point is inside a polygon is an independent task. — Victor Deleau, Nov 22 '19 at 19:31
Are they regular hexagons? Are they all of the same size? Are they in a hexagonal grid? Are they allowed to intersect? — kaya3, Nov 22 '19 at 20:57
Possible duplicate of [Looking for a fast way to find the polygon a point belongs to using Shapely](https://stackoverflow.com/questions/20297977/looking-for-a-fast-way-to-find-the-polygon-a-point-belongs-to-using-shapely) — Georgy, Nov 23 '19 at 10:40

JohanC · Answer 1 · 2019-11-22T20:46:21.793

Some ideas to speed things up:

Don't convert your hexagons to polygon = Polygon(poly_hex[i]) at each step of the inner loop. Calculate them once before the loop and store them in a list.
As the surface of a hexagon is close to a circle, also make a list of just the centre and the radius of the hexagons. Inside the loop, first test the distance between the point p and the centre of each hexagon, comparing to the radius. Only in case p is less than the radius to the centre, do an explicit test to see whether it is inside the hexagon.
If this still isn't fast enough, and depending on how small, and how far apart the hexagons are relative to your points, you could speed things up further with a quadtree or a grid. But this is less trivial. If most of the points are outside all hexagons, you could first test whether p is inside the convex hull of the hexagons.

Some remarks:

I don't understand your zone[j] = int(i). You convert i to an integer, while it already is an integer. Also, the way you store this, means a point can only be inside maximum one hexagon. Therefore, you can stop the loop as soon as you encountered such a hexagon. (There is also the rare case where a point is just on the edge between two hexagon. This might assign it to the first hexagon encountered. Or, due to rounding errors, even might have it outside both hexagons.)
If you want to save the zones as integers, create it as: zones = np.zeros(n, dtype=int)-1.
If your hexagons form a regular hexagonal grid, you can find the exact hexagon simply from the x and y coordinates without the need for in-polygon testing.

Thanks for the ideas. For the remarks, without int(i) it was saving the zone numbers as float (1.0, 2.0, ...), which was not what I wanted. That's correct, the hexagons don't layover so each point is at max in one hexagon. I am actually breaking the inner loop when the first hexagon is found. — Nima Chitsazan, Nov 22 '19 at 20:31
The zones are still saved as floats, because you have a general numpy array (which gets type np.float64). Try `zones = np.zeros(n, dtype=int)-1` to get integers. — JohanC, Nov 22 '19 at 20:42

Points in which polygon

1 Answers1