Find nearest polygon (from GeoSeries) to point (from GeoSeries)

Question

I have two GeoSeries, consisting of points and polygons. I want to find the polygon in dataframe B that is closest to each point in dataframe A. The polygons are rooftops from https://github.com/Microsoft/USBuildingFootprints, which I already geocoded using https://github.com/Bonsanto/polygon-geohasher.

I'm currently computing the 7 digit geohash of each point, and merging on buildings in neighboring 7 digit geohashes using geotools.expand. This is better than doing a full outer merge, but relies on explode. My general approach was to minimize calls of distance, since computing the distance from a point to a polygon is expensive.

The code is a bit slow (~20 minutes to match 100k rows), and so I'm trying to make it faster. My searching points to r-trees, but the sklearn implementation seems to be geared towards identifying the nearest point, rather than the nearest polygon. I'm interested in the left join rather than the right join.

Code below:

import pandas as pd
import numpy as np
import geopandas
import geohash
from shapely.geometry import Point

def match_func(df):
    point = Point(df.iloc[0,:][['lat', 'long']])
    df.loc[:, 'dist'] = geopandas.GeoSeries(df.geometry).distance(point)
    df = df.sort_values('dist')
    return(df.head(1))

def main(file):

    x           = import_points()
    rooftop_df  = import_rooftops()

    x['id'] = range(1, len(x) + 1)

    def neighbor_fun(lat,long):
         return(geohash.encode(lat,long,precision=7))

    func1 = np.vectorize(neighbor_fun)

    x['g7_neighbor'] = func1(x['lat'], x['long'])
    x = x.explode('g7_neighbor')
    x = x.merge(rooftop_df, left_on='g7_neighbor', right_on='geo7')

    xg = x.groupby('id')
    xout = pd.concat([match_fun2(df_group) for group_name, df_group in xg])
    return(xout)

Answer 1

Following Stefan's suggestion (Find nearest polygon (from GeoSeries) to point (from GeoSeries)), I wrote up the following using STRtree from Shapely.

Profiling indicates that, even ignoring the overhead associated with creating the tree, each nearest call is multiple seconds, so that it's slower than my old method. I suspect this is because the tree contains every rooftop, but have to think more about how to subset the points and rooftops data intelligently to get smaller trees.

import pandas as pd
import numpy as np
import geopandas
import geohash
from shapely.geometry import Point
from shapely.strtree import STRtree

def main(file):

    x           = import_points()
    rooftop_df  = import_rooftops()

    x['row'] = range(1, len(x) + 1)

    tree      = STRtree(rooftop_df.geometry)
    row_ dict = pd.Series(rooftop_df.row.values,index=rooftop_df.geometry.apply(id)).to_dict()

    def func1(x, y):
        return(accpt_dict[id(tree.nearest(Point(x,y)))])

    func1_v = np.vectorize(func1)

    x['accpt_id'] = func1_v(x.lat, x.long)
    return(x)

Answer 2

This will only work if you use a projected CRS since it's computing the euclidean distance but it's worth giving a shot if you can reproject.

import numpy as np
from shapely.geometry import Polygon, Point, MultiPolygon

def prep_polygons_asarr(gs):
    def get_pts(poly):
        if isinstance(poly, Polygon):
            coords = np.array(poly.exterior.coords)
        elif isinstance(poly, MultiPolygon):
            coords = np.concatenate([get_pts(sp) for sp in poly.geoms])
        return coords

    return [get_pts(poly) for poly in gs]


def get_nearest_poly(pt, polys):
    pt = np.array(pt.coords)
    dists = np.array([np.abs(np.linalg.norm(poly - pt, axis=1)).min() for poly in polys])
    return dists.argmin()


pt = Point(0, 1)

p1 = Polygon([( 7, 10), ( 7,  3), (10,  4), ( 7, 10)])
p2 = Polygon([(10, 10), (11, 11), (13, 10), (11, 15), (10, 10)])
p3 = Polygon([(13, 10), (14, 10), (15, 11), (12, 22), (13, 10)])

polys = prep_polygons_asarr([p1, p2, p3])

get_nearest_poly(pt, polys)

I tested on a dataset with 44k polygons and each call of get_nearest_poly() took ~.9-1.5s. Prepping the polygons took ~5-6s.