RTree: Count points in the neighbourhoods within each point of another set of points

Question

Why is this not returning a count of number of points in each neighbourhoods (bounding box)?

import geopandas as gpd

def radius(points_neighbour, points_center, new_field_name, r):
    """
    :param points_neighbour:
    :param points_center:
    :param new_field_name: new field_name attached to points_center
    :param r: radius around points_center
    :return:
    """
    sindex = points_neighbour.sindex
    pts_in_neighbour = []
    for i, pt_center in points_center.iterrows():
        nearest_index = list(sindex.intersection((pt_center.LATITUDE-r, pt_center.LONGITUDE-r, pt_center.LATITUDE+r, pt_center.LONGITUDE+r)))
        pts_in_this_neighbour = points_neighbour[nearest_index]
        pts_in_neighbour.append(len(pts_in_this_neighbour))
    points_center[new_field_name] = gpd.GeoSeries(pts_in_neighbour)

Every loop gives the same result.

Second question, how can I find k-th nearest neighbour?

More information about the problem itself:

• We are doing it at a very small scale e.g. Washington State, US or British Columbia, Canada

• We hope to utilize geopandas as much as possible since it's similar to pandas and supports spatial indexing: RTree

• For example, sindex here has method nearest, intersection, etc.

Please comment if you need more information. This is the code in class GeoPandasBase

@property
def sindex(self):
    if not self._sindex_generated:
        self._generate_sindex()
    return self._sindex

I tried Richard's example but it didn't work

def radius(points_neighbour, points_center, new_field_name, r):
    """
    :param points_neighbour:
    :param points_center:
    :param new_field_name: new field_name attached to points_center
    :param r: radius around points_center
    :return:
    """
    sindex = points_neighbour.sindex
    pts_in_neighbour = []
    for i, pt_center in points_center.iterrows():
        pts_in_this_neighbour = 0
        for n in sindex.intersection(((pt_center.LATITUDE-r, pt_center.LONGITUDE-r, pt_center.LATITUDE+r, pt_center.LONGITUDE+r))):
            dist = pt_center.distance(points_neighbour['geometry'][n])
            if dist < radius:
                pts_in_this_neighbour = pts_in_this_neighbour + 1
        pts_in_neighbour.append(pts_in_this_neighbour)
    points_center[new_field_name] = gpd.GeoSeries(pts_in_neighbour)

To download the shape file, goto https://catalogue.data.gov.bc.ca/dataset/hellobc-activities-and-attractions-listing and choose ArcView to download

Answer 1

I've attached code which should, with some minor modifications, do what you want.

I think your problem arose for one of two reasons:

1. You were not correctly constructing the spatial index. Your responses to my comments suggested that you weren't wholly aware of how the spatial index was getting made.

2. The bounding box for your spatial query was not built correctly.

I'll discuss both possibilities below.

Constructing the spatial index

As it turns out, the spatial index is constructed simply by typing:

sindex = gpd_df.sindex

Magic.

But from whence does gpd_df.sindex get its data? It assumes that the data is stored in a column called geometry in a shapely format. If you have not added data to such a column, it will raise a warning.

A correct initialization of the data frame would look like so:

#Generate random points throughout Oregon
x = np.random.uniform(low=oregon_xmin, high=oregon_xmax, size=10000)
y = np.random.uniform(low=oregon_ymin, high=oregon_ymax, size=10000)

#Turn the lat-long points into a geodataframe
gpd_df = gpd.GeoDataFrame(data={'x':x, 'y':y})
#Set up point geometries so that we can index the data frame
#Note that I am using x-y points!
gpd_df['geometry'] = gpd_df.apply(lambda row: shapely.geometry.Point((row['x'], row['y'])), axis=1)

#Automagically constructs a spatial index from the `geometry` column
gpd_df.sindex 

Seeing the foregoing sort of example code in your question would have been helpful in diagnosing your problem and getting going on solving it.

Since you did not get the extremely obvious warning geopandas raises when a geometry column is missing:

AttributeError: No geometry data set yet (expected in column 'geometry'.

I think you've probably done this part right.

Constructing the bounding box

In your question, you form a bounding box like so:

nearest_index = list(sindex.intersection((pt_center.LATITUDE-r, pt_center.LONGITUDE-r, pt_center.LATITUDE+r, pt_center.LONGITUDE+r)))

As it turns out, bounding boxes have the form:

(West, South, East, North)

At least, they do for X-Y styled-points, e.g. shapely.geometry.Point(Lon,Lat)

In my code, I use the following:

bbox = (cpt.x-radius, cpt.y-radius, cpt.x+radius, cpt.y+radius)

Working example

Putting the above together leads me to this working example. Note that I also demonstrate how to sort points by distance, answering your second question.

#!/usr/bin/env python3

import numpy as np
import numpy.random
import geopandas as gpd
import shapely.geometry
import operator

oregon_xmin = -124.5664
oregon_xmax = -116.4633
oregon_ymin = 41.9920
oregon_ymax = 46.2938

def radius(gpd_df, cpt, radius):
  """
  :param gpd_df: Geopandas dataframe in which to search for points
  :param cpt:    Point about which to search for neighbouring points
  :param radius: Radius about which to search for neighbours
  :return:       List of point indices around the central point, sorted by
                 distance in ascending order
  """
  #Spatial index
  sindex = gpd_df.sindex
  #Bounding box of rtree search (West, South, East, North)
  bbox = (cpt.x-radius, cpt.y-radius, cpt.x+radius, cpt.y+radius)
  #Potential neighbours
  good = []
  for n in sindex.intersection(bbox):
    dist = cpt.distance(gpd_df['geometry'][n])
    if dist<radius:
      good.append((dist,n))
  #Sort list in ascending order by `dist`, then `n`
  good.sort() 
  #Return only the neighbour indices, sorted by distance in ascending order
  return [x[1] for x in good]

#Generate random points throughout Oregon
x = np.random.uniform(low=oregon_xmin, high=oregon_xmax, size=10000)
y = np.random.uniform(low=oregon_ymin, high=oregon_ymax, size=10000)

#Turn the lat-long points into a geodataframe
gpd_df = gpd.GeoDataFrame(data={'x':x, 'y':y})
#Set up point geometries so that we can index the data frame
gpd_df['geometry'] = gpd_df.apply(lambda row: shapely.geometry.Point((row['x'], row['y'])), axis=1)

#The 'x' and 'y' columns are now stored as part of the geometry, so we remove
#their columns in order to save space
del gpd_df['x']
del gpd_df['y']

for i, row in gpd_df.iterrows():
  neighbours = radius(gpd_df,row['geometry'],0.5)
  print(neighbours)
  #Use len(neighbours) here to construct a new row for the data frame

(What I had been requesting in the comments is code that looks like the foregoing, but which exemplifies your problem. Note the use of random to succinctly generate a dataset for experimentation.)