This is a follow up to my earlier question (Spatial join in R - Adding points to polygons with multiple corresponding points).

I have successfully joined a spatial points file to a polygon file in R using the st_join function within the sf package with more than one point being assigned to a polygon if necessary, duplicating rows but keeping all points which fall within a polygon.

st_join(polygons, points)

However I also need to join points which fall outside of the polygons but within 500m of a polygon to their nearest polygon. Points which are >500m away from a polygon can be discarded.

I thought that combining the above with st_nn from the nngeo package should work using the following:

st_join(polygons, points, join = st_nn, maxdist = 500)

However in this case only 1 point is assigned to a polygon, even if more than one point falls within a polygon or within 500m of a polygon. i.e. the rows are not duplicated.

Here is a screenshot of a sample of points and polygons:

enter image description here

And here is table showing how the points should be assigned to the polygons and how they have been assigned in the respective methods:

enter image description here

I find it a little strange that the second method does not keep the duplicates, even though it is based on the same function. Can anyone tell me what I'm doing wrong here?

Edit: I tried adjusting the k parameter but this simply joins the first points within the given distance up to the max number given and therefore can assign 1 point to 2 polygons. e.g.

st_join(polygons, points, join = st_nn, k = 10, maxdist = 500)

returns 5 points for polygon 89028 as there are 5 points within 500m, when in fact only 1 point should be returned (011-05-0529) as the other 4 points are already assigned to the polygon in which they fall. A point should only be assigned to one polygon.

  • 1
    Can you use st_distance to create a distance matrix - returns 0 if point is in polygon otherwise returns distance to polygon...
    – Spacedman
    Apr 25, 2019 at 14:49
  • @Spacedman, I'm not sure this would make any difference. 'join = st_nn, maxdist = 500' already finds all the points within 500m of the polygons. The problem is with how points within and within 500m of polygons are assigned when more than one point is within or within 500m of a polygon. If you disagree, could you please elaborate on your answer?
    – Adam G
    May 2, 2019 at 12:02

2 Answers 2


If I understood correctly, you find the containing polygon of each point, or else the nearest polygon (up to 500m) if the point is not contained inside any polygon.

If so, the following expression, where the order of x and y is reversed, should work -

st_join(points, polygons, join = st_nn, k = 1, maxdist = 500)

The function will look for the nearest polygon from each point. The containing polygon, if any, is always considered to be nearest since its distance from the point is zero. If no containing polygon is found, the function will look for the nearest polygon, up to a maximal distance of 500m.

  • You understand correctly but reversing the order of x and y results in a point file. The resulting file should be a polygon file. This will have duplicate rows as one polygon can be associated with many points but each point can only be associated with one polygon. Using st_join(polygons, points) provides the correct output for those points found within the polygons. Adding st_join(polygons, points, join = st_nn, maxdist = 500) changes the output so that only one point is assigned to each polygon. Increasing k results in points being added to more than one polygon.
    – Adam G
    May 3, 2019 at 9:34
  • You are right that the result will be a point layer. But the important thing is that the layer will have the point IDs and polygon IDs that you want to match. Afterwards, the point geometry can be deleted (which leaves an ordinary data.frame) and the polygon geometry can be joined back using the polygon ID, to get the final sf object. May 3, 2019 at 14:57
  • If you can post sample data and the expected result I will prepare an example. May 3, 2019 at 14:58
  • I realised that this would be a possible solution but I was still hoping that there would be an easier way. I'll work through the process tomorrow and repost to comfirm if it worked. If not, I'll post some sample data for you to take a look at. Thanks.
    – Adam G
    May 6, 2019 at 16:43
  • Is there a quicker process than this? I know we're doing Geo-processing here, but often the purpose of doing this in R is to save time AND effort. Thanks
    – obrunt
    May 24 at 12:55

I used a combination of @Michael's answer plus some additional manipulation to get the correct format. The resulting file is a polygon file with no duplicate polygons. If a polygon has >1 associated point, then the point columns from the join are repeated until every associated point is included.


#Load files
Poly <-st_read("Path/Poly.shp")
Pts <- st_read("Path/Pts.shp")

names(Pts) #Get list of names for selecting required columns
col_interest <- c("col1", "col2") #add column names here

Join Pts to Poly resulting in a pts file with the ID of the nearest polygon within 500m attached in the polygon attributes
Poly_Pts_pts <- st_join(Pts, Poly, join = st_nn, maxdist = 500)

#convert to data.table
Poly_Pts_pts_DT <- as.data.table(Poly_Pts_pts)

#add a new column with running number for each individual Point within each polygon ID
Poly_Pts_pts_DT <-  Poly_Pts_pts_DT[, New_ID := seq_len(.N), by = ID]

#Cast into wide format
Poly_Pts_pts_wide <- dcast.data.table(Poly_Pts_pts_DT, ID ~ New_ID, value.var = col_interest)#output is data.table

#Join Pts wide format to original polygons on ID column
Poly_Pts <- merge(Poly, Poly_Pts_pts_wide, by = "ID", all.x = TRUE)

#Write to disk
st_write(Poly_Pts, "Path/Poly_Pts.shp")```

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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