Finding the closest line to polygon using Python Rtree

Question

I have a series of fields and roads. I am trying to find the nearest road to each field. I want to add field to each of my polygons with the index/name etc. of the closest road.

Would this be the correct method to find the nearest road for each field?

Below is the code I am trying

So I have loaded my shapes and created an empty index.

fieldlyr = fiona.open(fields, 'r')
roadlyr = fiona.open(roads, 'r')
# search radius, i've changed sizes to make it large enough
areaft = 2500
index = rtree.index.Index()

I add my fields to the index.

for fid, feat in fieldlyr.items():
  geometry = shape(feat['geometry'])
  index.insert(fid, geometry.bounds)

I then buffer my roads

for feat in roadlyr:
  geometry = shape(feat['geometry'])
  geometry_buffered = geometry.buffer(areaft)

Then I believe this will give me the lists of intersections.

fids = [int(i) for i in index.intersection(geometry_buffered.bounds)]

How would I get road "id" A is nearest to polygon "id" 2 and so on.

Is this even the correct method?

Answer 1

So it seems i was over complicating things.

this appears to work how i need it to. Just thought i would share in case anyone has a similar need or perhaps a better method

fieldlyr = fiona.open(fields, 'r')
roadlyr = fiona.open(roads, 'r')

# dont need the search radius/ buffer for near
index = rtree.index.Index()

Correct on adding to index, this should be the road layer since i want to know which road is closest to an area

for fid, feat in roadlyr.items():
  geom = shape(feat['geometry'])
  index.insert(fid, geom.bounds)

So i had the roads and area reveresed in my original question

for afeat in fieldlyr:
  geometry = shape(afeat['geometry'])

  fids = [int(i) for i in idx.nearest(geometry.bounds)]

And to view the matches

print (str(afeat['id']) + " area is nearest to line "+ str(fids))

As if now, this seems to work

Answer 2

If you want to use the rtree to speed up the spatial join, i would use the GeoPandas built-in function ".sindex".

import geopandas as gpd

# set crs
gdf_lines = gdf_lines.set_crs(YOURCRS, allow_override=True)
gdf_polygons = gdf_polygons.set_crs(YOURCRS, allow_override=True)

# index the data
gdf_lines.sindex
gdf_polygons.sindex

# spatial join
gdf_sjoin_nearest = gpd.sjoin_nearest(gdf_polygons, 
                                        gdf_lines, 
                                        how="left", 
                                        lsuffix='_left', 
                                        rsuffix='_right')

Note that according to the Geopandas Documentation GeoPandas offers built-in support for spatial indexing using an R-Tree algorithm. Depending on the ability to import pygeos, GeoPandas will either use pygeos.STRtree or rtree.index.Index.