St_Hexagongrid Geom Vector to Find All Points

Caculate point 50 miles away (North, 45% NE, 45% SW)

Try combining ST_Project with a CTE - adjust the values of radians to the azimuth you need.

WITH j AS (
  SELECT poi::geography AS poi FROM t
)
SELECT 
  ST_AsText(ST_Project(j.poi, 80467.2, radians(90.0)),2),
  ST_AsText(ST_Project(j.poi, 80467.2, radians(45.0)),2),
  ST_AsText(ST_Project(j.poi, 80467.2, radians(180.0)),2),
  ST_AsText(ST_Project(j.poi, 80467.2, radians(135.0)),2),
  ST_AsText(ST_Project(j.poi, 80467.2, radians(270.0)),2),
  ST_AsText(ST_Project(j.poi, 80467.2, radians(225.0)),2),
  ST_AsText(ST_Project(j.poi, 80467.2, radians(360.0)),2),
  ST_AsText(ST_Project(j.poi, 80467.2, radians(315.0)),2)
FROM j;

      st_astext      |      st_astext      |    st_astext     |     st_astext      |      st_astext      |     st_astext      |    st_astext     |      st_astext      
---------------------+---------------------+------------------+--------------------+---------------------+--------------------+------------------+---------------------
 POINT(-73.05 40.71) | POINT(-73.32 41.22) | POINT(-74 39.99) | POINT(-73.33 40.2) | POINT(-74.95 40.71) | POINT(-74.67 40.2) | POINT(-74 41.43) | POINT(-74.68 41.22)
(1 Zeile)

Note: The buffer (circle) in the image is just for illustration.

Find neighboring polygons with maximum of 3 other polygons

Your hint with ST_Touches is correct, however to get the amount of neighbor cells from one column related to other records in the same table you either need to run a subquery or call the table twice in the FROM clause.

Given the following grid on a table called tb ..

.. you can filter the cells with three neighbor cells or less like this:

SELECT * FROM tb q1
WHERE (
  SELECT count(*) FROM tb q2
  WHERE ST_Touches(q2.geom,q1.geom)) <=3;

If you want to also list which are the neighbor cells you might wanna first join the cells that touch in the WHERE clause and in a subquery or CTE count the results:

WITH j AS (
  SELECT 
    q1.poly_name AS p1,q2.poly_name p2, 
    COUNT(*) OVER (PARTITION BY q1.poly_name) AS qt
  FROM tb q1, tb q2
  WHERE ST_Touches(q2.geom,q1.geom))
SELECT * FROM j
WHERE qt <= 3;

Demo: db<>fiddle

Further reading:

• Create Hexagons (maybe relevant for your project)
• Window Functions

How to write id of nearby areas in PostgreSql?

I am not quite sure about your tablenames and about geometry. But if you replace the tablename and <distance_function(c.geom, c2.geom)> to something giving you the distance between the cities, this might help:

UPDATE cities c
SET d1_nearest = (
    SELECT c2.id
    FROM cities c2
    WHERE c2.d1 = 1
    ORDER BY <distance_function(c.geom, c2.geom)>
    LIMIT 1
) WHERE d1 = 0;

For the other two columns it is analogous. Of course you need to create the new columns in advance.

Find Minimum number of Circles with 50-Mile Radius that Covers All Points

Updated answer based on comments:

You have many options.

Here are 3 differents ways of doing that:

1. With scipy.CKDTree:

Pros :

• This will be fast

Cons :

• less accurate because the computed distance is euclidean
• and the radius will be the same as your inputs, so here in degrees

I would go with a cKDTree and a radius query to find all points in radius, remove theses points from list, and continue with remaining points. This is not optimal but can be a good basis.

from scipy.spatial import cKDTree

points = [(41.81014,-72.550028), (41.995833,-72.581525), (41.377211,-72.150307), (41.710626,-72.763862), (41.55254,-72.815454), (41.415022,-73.401914), (41.0554,-73.54142), (41.660572,-72.725673), (41.350949,-72.871673), (41.280278,-72.987515), (41.23354,-73.151677), (41.235174,-73.038092), (41.58254,-73.034321), (41.89121,-72.6521), (41.340446,-73.078943), (41.81886,-73.0755), (41.228735,-73.225326), (41.839019,-71.883778), (41.585192,-71.99693), (41.611472,-72.901357), (41.783976,-72.748229), (43.634242,-70.347774), (44.842191,-68.74156), (43.934038,-69.985271), (43.474,-70.5141), (44.312403,-69.804993), (42.552616,-70.937616), (41.877743,-71.068577), (41.940344,-71.351931), (42.399035,-71.071855), (42.168221,-72.642232), (42.518609,-71.135461), (42.160827,-71.498868), (42.481583,-71.024154), (42.305328,-71.398387), (42.29247,-71.7751), (41.796058,-71.321145), (42.376695,-71.090028), (42.364178,-71.156462), (41.971125,-70.716858), (42.280435,-71.655929), (42.359487,-71.607159), (42.503468,-70.919421), (42.194395,-71.774687), (42.357311,-72.547241), (42.328872,-71.062845), (42.033714,-71.310581), (42.39976,-71.000326), (42.527193,-71.71374), (42.495264,-73.206116), (41.63729,-71.003268), (42.110519,-70.927683), (42.152383,-71.073541), (42.02714,-71.1438), (42.740784,-71.161323), (41.773672,-70.745562), (42.788072,-71.115959), (42.623622,-71.318304), (42.137401,-70.83883), (42.348748,-71.504967), (41.749066,-71.207427), (42.2045,-71.1553), (42.22142,-71.021844), (42.589718,-71.159895), (42.344172,-71.099961), (42.364561,-71.102575), (42.2882972,-71.1267483), (42.350679,-71.114022), (42.494932,-71.103401), (42.42072,-71.09902), (42.388648,-71.118659), (42.484104,-71.186185), (41.666927,-70.294616), (42.275401,-71.029299), (42.299241,-71.062748), (42.361045,-71.0626), (42.764475,-71.215039), (43.2189,-71.485199), (42.702771,-71.437791), (43.045615,-71.461202), (42.79899,-71.53679), (42.941002,-71.473513), (42.928188,-72.301906), (43.235048,-70.884519), (43.048951,-70.818587), (43.633682,-72.322002), (44.466154,-73.18226)]

# Radius of circle. Note that the unit is the same as in your list (here, degrees.)
radius = 1

num_circles = 0

list_is_no_empty = True

while(list_is_no_empty):

    # Take the first point in order to find all points within distance radius
    start_point = points[0]

    # Create a KDTree
    tree = cKDTree(points)

    # Find indexes of all points in radius
    indexes_of_points_in_radius = tree.query_ball_point(start_point, radius)

    # Create the list of points to remove (points that were found within distance radius)
    points_to_remove = [points[i] for i in indexes_of_points_in_radius]

    # Remove these points
    points = list(set(points) - set(points_to_remove))

    # Increment the number of circles
    num_circles += 1

    # If no points remain, exit loop
    if points == []:
        list_is_no_empty = False

print("Number of circles:", num_circles)

2. With sklearn.neighbors.BallTree:

Pros:

• This will be more accurate because the computed distance here is the great-circle distance (Haversine formula).

Cons:

• Like the cKDTree, the radius will be the same as your inputs, so here in degrees.
• A little slower than scipy.cKDTree (2 times slower when I tested)

Note too that I found that some recommend to convert your inputs in radians because that is required for the Haversine formula (https://scikit-learn.org/stable/modules/generated/sklearn.neighbors.DistanceMetric.html). But from my testings, with scikit-learn 1.0.1, that wasn't needed. But just in case, you would do:

from math import radians
points = [tuple(map(radians, point)) for point in points]
start_point = (radians(start_point[0]), radians(start_point[1]))
radius = radians(radius)

Code with BallTree:

from sklearn.neighbors import BallTree
import numpy as np

num_circles = 0

list_is_no_empty = True

while(list_is_no_empty):

    # Take the first point in order to find all points within distance radius
    start_point = np.array([points[0]])

    # Create a BallTree, and chose the Haversine formula (great circle distance)
    tree = BallTree(points, metric="haversine")

    # Find indexes of all points in radius
    indexes_of_points_in_radius = tree.query_radius(start_point, r=radius)[0]

    # Create the list of points to remove (points that were found within distance radius)
    points_to_remove = [points[i] for i in indexes_of_points_in_radius]

    # Remove these points
    points = list(set(points) - set(points_to_remove))

    # Increment the number of circles
    num_circles += 1

    # If no points remain, exit loop
    if points == []:
        list_is_no_empty = False

print("Number of circles:", num_circles)

3. With sklearn.neighbors.BallTree, using a user-defined distance function:

Pros :

• We will be able here to use a very accurate distance
• We will be able to specify this distance in miles (or meters)

Cons:

• Way slower than cKDTree (10 times when I tested)

from pyproj import Geod
from sklearn.neighbors import BallTree
import numpy as np

# Create a WGS84 ellipsoid
geod = Geod(ellps='WGS84')

# User-defined function for BallTree
# We use the "inv" method of pyproj in order to get the distance in meters between 2 points
# It computes the geodesic distance using the wonderful Karney's algorithm
def geodedsic_distance(point_01, point_02):
    lat1,lon1 = point_01
    lat2,lon2 = point_02
    az12,az21,distance_in_meters = geod.inv(lon1,lat1,lon2,lat2)
    return distance_in_meters

def miles_to_meters(miles):
    return miles * 1609.344

# Radius in miles
radius_in_miles = 50

radius_in_meters = miles_to_meters(50)

num_circles = 0

list_is_no_empty = True

while(list_is_no_empty):

    # Take the first point in order to find all points within distance radius
    start_point = np.array([points[0]])

    # Create a BallTree, and chose our custom function
    tree = BallTree(points, metric=geodedsic_distance)

    # Find indexes of all points in radius, specified in meters
    indexes_of_points_in_radius = tree.query_radius(start_point, r=radius_in_meters)[0]

    # Create the list of points to remove (points that were found within distance radius)
    points_to_remove = [points[i] for i in indexes_of_points_in_radius]

    # Remove these points
    points = list(set(points) - set(points_to_remove))

    # Increment the number of circles
    num_circles += 1

    # If no points remain, exit loop
    if points == []:
        list_is_no_empty = False

print("Number of circles:", num_circles)

If you want to learn more about miles to degrees conversion (and why, in fact, we can't) and computing distances on earth:

Is the Haversine Formula or the Vincenty's Formula better for calculating distance?

https://gis.stackexchange.com/questions/84885/difference-between-vincenty-and-great-circle-distance-calculations

---

Related Topics