Create polygon from list of points only if they are nearby

Question

I have a list of points (longitude and latitude), as well as their associated point geometries in a geodataframe. All of the points should be able to be subdivided into individual polygons, as the points are generally clustered in several areas. What I would like to do is have some sort of algorithm that loops over the points and checks the the distance between the previous and current point. If the distance is sufficiently small, it would group those points together. This process would continue until the current point is too far away. It would make a polygon out of those close points, and then continue the process with the next group of points.

gdf
longitude   latitude    geometry
0   -76.575249  21.157229   POINT (-76.57525 21.15723)
1   -76.575035  21.157453   POINT (-76.57503 21.15745)
2   -76.575255  21.157678   POINT (-76.57526 21.15768)
3   -76.575470  21.157454   POINT (-76.57547 21.15745)
5   -112.973177 31.317333   POINT (-112.97318 31.31733)
... ... ... ...
2222    -113.492501 47.645914   POINT (-113.49250 47.64591)
2223    -113.492996 47.643609   POINT (-113.49300 47.64361)
2225    -113.492379 47.643557   POINT (-113.49238 47.64356)
2227    -113.487443 47.643142   POINT (-113.48744 47.64314)
2230    -105.022627 48.585669   POINT (-105.02263 48.58567)

So in the data above, the first 4 points would be grouped together and turned into a polygon. Then, it would move onto the next group, and so forth. Each group of points is not evenly spaced, i.e., the next group might be 7 pairs of points, and the following could be 3. Ideally, the final output would be another geodataframe that is just a bunch of polygons.

please provide more data - the approach depends critically on how distinct the clusters are — anon01, Mar 06 '21 at 19:16
How do you want the polygons created, by convex hull or along original sequential points? — swatchai, Mar 07 '21 at 16:03

forgetso · Accepted Answer · 2021-03-02T14:04:38.417

You can try DBSCAN clustering as it will automatically find the best number of clusters and you can specify a maximum distance between points ( ε ).

Using your example, the algorithm identifies two clusters.

import pandas as pd
from sklearn.cluster import DBSCAN

df = pd.DataFrame(
    [
        [-76.575249, 21.157229, (-76., 21.15723)],
        [-76.575035, 21.157453, (-76.57503, 21.15745)],
        [-76.575255, 21.157678, (-76.57526, 21.15768)],
        [-76.575470, 21.157454, (-76.57547, 21.15745)],
        [-112.973177, 31.317333, (-112.97318, 31.31733)],
        [-113.492501, 47.645914, (-113.49250, 47.64591)],
        [-113.492996, 47.643609, (-113.49300, 47.64361)],
        [-113.492379, 47.643557, (-113.49238, 47.64356)],
        [-113.487443, 47.643142, (-113.48744, 47.64314)],
        [-105.022627, 48.585669, (-105.02263, 48.58567)]
    ], columns=["longitude", "latitude", "geometry"])

clustering = DBSCAN(eps=0.3, min_samples=4).fit(df[['longitude','latitude']].values)
gdf = pd.concat([df, pd.Series(clustering.labels_, name='label')], axis=1)
print(gdf)
gdf.plot.scatter(x='longitude', y='latitude', c='label')

    longitude   latitude                geometry  label
0  -76.575249  21.157229       (-76.0, 21.15723)      0
1  -76.575035  21.157453   (-76.57503, 21.15745)      0
2  -76.575255  21.157678   (-76.57526, 21.15768)      0
3  -76.575470  21.157454   (-76.57547, 21.15745)      0
4 -112.973177  31.317333  (-112.97318, 31.31733)     -1 # not in cluster
5 -113.492501  47.645914   (-113.4925, 47.64591)      1
6 -113.492996  47.643609    (-113.493, 47.64361)      1
7 -113.492379  47.643557  (-113.49238, 47.64356)      1
8 -113.487443  47.643142  (-113.48744, 47.64314)      1
9 -105.022627  48.585669  (-105.02263, 48.58567)     -1 # not in cluster

If we add random data to your data set, run the clustering algorithm, and filter out those data points not in clusters, you get a clearer idea of how it's working.

import numpy as np
rng = np.random.default_rng(seed=42)
arr2 = pd.DataFrame(rng.random((3000, 2)) * 100, columns=['latitude', 'longitude'])
randdf = pd.concat([df[['latitude', 'longitude']], arr2]).reset_index()
clustering = DBSCAN(eps=1, min_samples=4).fit(randdf[['longitude','latitude']].values)
labels =  pd.Series(clustering.labels_, name='label')
gdf = pd.concat([randdf[['latitude', 'longitude']], labels], axis=1)

subgdf = gdf[gdf['label']> -1]
subgdf.plot.scatter(x='longitude', y='latitude', c='label', colormap='viridis', figsize=(20,10))

print(gdf['label'].value_counts())
-1     2527
 16      10
 3        8
 10       8
 50       8
       ... 
 57       4
 64       4
 61       4
 17       4
 0        4
Name: label, Length: 99, dtype: int64

Getting the clustered points from this dataframe would be relatively simple. Something like this:

subgdf['point'] = subgdf.apply(lambda x: (x['latitude'], x['longitude']), axis=1)
subgdf.groupby(['label'])['point'].apply(list)

label
0     [(21.157229, -76.575249), (21.157453, -76.5750...
1     [(47.645914, -113.492501), (47.643609, -113.49...
2     [(46.67210037270342, 4.380376578722878), (46.5...
3     [(85.34030732681661, 23.393948586534073), (86....
4     [(81.40203846660347, 16.697291990770392), (82....
                            ...                        
93    [(61.419880354359925, 23.25522624430636), (61....
94    [(50.893415175135424, 90.70863269095085), (52....
95    [(88.80586950148697, 81.17523712192651), (88.6...
96    [(34.23624333000541, 40.8156668231013), (35.86...
97    [(16.10456828199399, 67.41443008931344), (15.9...
Name: point, Length: 98, dtype: object

Although you'd probably need to do some kind of sorting to make sure you were connecting the closest points when drawing the polygons.

Similar SO question

DBSCAN from sklearn

CarmeloAlvarez · Answer 2 · 2021-03-02T17:14:49.877

Haversine Formula in Python (Bearing and Distance between two GPS points)

https://gis.stackexchange.com/questions/121256/creating-a-circle-with-radius-in-metres

You may be able to use the haversine formula to group points within a distance. Create polygons for each point (function below) with the formula then filter points inside from the master list of points and repeat until there are no more points.

#import modules
import numpy as np
import pandas as pd
import geopandas as gpd
from geopandas import GeoDataFrame, GeoSeries
from shapely import geometry
from shapely.geometry import Polygon, Point
from functools import partial
import pyproj
from shapely.ops import transform

#function to create polygons on radius
def polycir(lat, lon, radius):
    
    local_azimuthal_projection = """+proj=aeqd +R=6371000 +units=m +lat_0={} +lon_0= 
                                  {}""".format(lat, lon)
    
    wgs84_to_aeqd = partial(
        pyproj.transform,
        pyproj.Proj("+proj=longlat +datum=WGS84 +no_defs"),
        pyproj.Proj(local_azimuthal_projection),
    )
    
    aeqd_to_wgs84 = partial(
        pyproj.transform,
        pyproj.Proj(local_azimuthal_projection),
        pyproj.Proj("+proj=longlat +datum=WGS84 +no_defs"),
    )

    center = Point(float(lon), float(lat))
    point_transformed = transform(wgs84_to_aeqd, center)
    buffer = point_transformed.buffer(radius)
    # Get the polygon with lat lon coordinates
    circle_poly = transform(aeqd_to_wgs84, buffer)
    
    return circle_poly

#Convert df to gdf
gdf = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df.longitude, df.latitude))

#Create circle polygons col
gdf['polycir'] = [polycir(x, y, <'Radius in Meters'>) for x, y in zip(gdf.latitude, 
                  gdf.longitude)]

gdf.set_geometry('polycir', inplace=True)


#You should be able to loop through the polygons and find the geometries that overlap with 
# gdf_filtered = gdf[gdf.polycir.within(gdf.iloc[0,4])]

score -1 · Answer 3 · answered Mar 02 '21 at 13:38

-1

Looks like a job for k-means clustering.

You may need to be careful to how you define your distance (actual disctance "through" earth, or shortest path around?)

Turning each cluster into a polygon depends on what you want to do... just chain them or look for their convex enveloppe...

answered Mar 02 '21 at 13:38

shermelin

1
1

2

KMeans requires you to know how many clusters you want to detect. – forgetso Mar 02 '21 at 14:00
2

It is usual to run k-means for several values of k and use some criterion to select some optimal k value. Of course, you need to define that criterion for "optimal". The OP does not clarify when clusters should be split, which size... which leaves open what is this measure. – shermelin Mar 04 '21 at 13:48

Create polygon from list of points only if they are nearby

3 Answers3