Combining st_join with st_nn to include all points within, as well as within a given distance of a polygon

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.

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.

library(sf)
library(data.table)
library(nngeo)

#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")```