CODEWITHSUNDEEP
pythonnumpyshapely
Alex
Alex

Reputation: 619

How can I use Shapely to detect all points that are closer than N meters?

Say I've got a list of points (the coordinates are in meters):

points = [(1, 1), (2, 2), (-1, 3), ..., (1000, 1000)]

And I want to use Shapely library to return all the points (from points) that are within N meters radius to some specified origin e.g., (3, 3).

Update: if I execute:

import numpy as np
points = [(1, 1), (2, 2), (3, 4), (100, 100), (5, 5), (1, 2)]
points_array = np.array(points)
print(points_array.shape)  # prints [n,2] where n is the number of points
max_distance = 2
origin = np.array([3, 3])
distance_array = np.sqrt((points_array - origin) ** 2)
near_points = points_array[distance_array < max_distance]
print('near points', near_points)

I receive 

(6, 2)
near points [2 2 3 4 2]

which looks a bit weird because I didn't expect odd number of elements even though if I get [(2, 2), (3, 4), 2=?]

Upvotes: 1

Views: 1642

Answers (3)

oekopez
oekopez

Reputation: 1090

Fast calculations for lots of data are best handled with numpy. Instead of creating shapely objects from the coordinates and use the built in distance function, it is much easier (for points!, not for polygons) to calculate the distance using a numpy array. If you want to do the same for lines and polygons and use the nearest distances, shapely is your friend.

import numpy as np
points_array = np.array(points)
print(points_array.shape)  # prints [n,2] where n is the number of points
max_distance = 100
origin = np.array([3, 3])

Now calculate the euclidean distances as the square sum of the differences over axis 1 (the coordinates)

distance_array = np.sqrt(np.sum((points_array - origin) ** 2, 1))

And retrieve the points where the distance is smaller than max_distance

near_points = points_array[distance_array < max_distance]

To compare the numpy solutions to the other answers in terms of speed I timed the answers for the same set of 1e6 random points:

  • The code above takes 49 ms
  • the optimized solution by Peter Collingridge: 44ms
  • List solution by vurmax (using list comprehension, see below): 2.88s (60x slower)
  • The list solution with Peter Collingridge's optimization: 2.48s
  • Toy shapely solution by Christian Sloper: 15.2s (300x slower)
points_vurmax = [(x,y) for x,y in points_array 
                 if math.sqrt((x - origin[0])**2 + (y - origin[1])**2) < max_dist]
points_vurmax2 = [(x,y) for x,y in points_array 
                  if (x - origin[0])**2 + (y - origin[1])**2 < max_dist ** 2]

Upvotes: 5

Christian Sloper
Christian Sloper

Reputation: 7510

Not going to be as fast as the numpy solutions, but here is one in shapely.

from shapely.geometry import MultiPoint

m = MultiPoint(points)
circle = Point(3,3).buffer(N)
[p for p in m if circle.covers(p) ]

returns the points as shapely points.

Upvotes: 2

vurmux
vurmux

Reputation: 10020

You don't need to use shapely. The distance function is just:

d = math.sqrt((x2 - x1)**2 + (y2 - y1)**2)

You don't even need numpy or any another library except std math

from math import sqrt

points = [(1, 1), (2, 2), (-1, 3)]
target = (3, 4)
distance = 3
for x, y in points:
    if math.sqrt((x - target[0])**2 + (y - target[1])**2) < distance:
        print(x, y)

(Of course you can use numpy, if you have the huge amount of points and need to process them really fast)

Upvotes: 0

Related Questions