8

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?

Fields and Roads

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?

2 Answers 2

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

1
  • 4
    Just remember that you're finding the bounding boxes nearest each other. If you need an exact result you'll want to follow up the rtree query with a shapely object.distance(other) routine
    – mikewatt
    Nov 19, 2019 at 19:11
1

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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.