Creating multi-ring buffers around points in a single layer and for each point selecting neighbors within a specified radius in R

Question

I have a CSV file of data on a set of 3000 points (locations) with respective coordinates and a number of attributes. I want to perform a regression analysis that takes cognizance of any spatial relationships among the points(locations). As such, I want to arrange the data in a dataframe in a way such that for each well, I will have the following:

1. A list of points (location IDs) within the specified radius of the buffers around each point(location), where the specified radius are for example: 0.5, 1, 2, 3, 4, 5 miles.

2. The sum and/or average of the values of given attributes(eg. population_size) of the point (location) that fall within each of these specified radii of each point (location).

I want to write the resulting table to a CSV file.

Below is what I tried based on some online resources.

library(raster)  # for handling geographic raster data
library(dplyr)
library(ggplot2)
library(stringr) # for working with strings (pattern matching)
library(rgdal) #For checking available CRSs interactively. 
library(lwgeom) #Needed for distances
library(geosphere)
library(data.table)
#library(maps)
library(reshape) 


#Load Data:
data<- read.csv("Loc.csv")

#Set as sf object
datatosf=st_as_sf(data,coords = c("lon", "lat"))

#Then set appropriate geometric CRS as
datatosf_geo = st_set_crs(datatosf, 4326)

#Then project into UTM
data_projected = st_transform(datatosf_geo, 26913) #26913 is for NAD UTM Zone 13

# create all possible pairs of origin-destination in a long format

newdata <- expand.grid.df(data_projected ,data_projected)  #Uisng the sf object data or projected data here gives me erro below:
Error in lapply(x[i], as.numeric) : 
 'list' object cannot be coerced to type 'double'

#But when I use the original csv data, data, I do not get this error

newdata <- expand.grid.df(data ,data)
names(newdata )[28:29] <- c("lat_dest","lon_dest")

# calculate distances in miles:
setDT(dtt)[ , dist_km := distGeo(matrix(c(lon, lat), ncol = 2), 
                                 matrix(c(lon_dest, lat_dest), ncol = 2))/1.609]

#Write results to a csv file
write.csv(newdata,'newdata.csv')

Now, clearly I did not create buffers. Instead what I have is a table with each point(location) and the distance from it to all other point s(locations) (including itself) -- that is a 9,000,000 rows of data. And my CSV file could return or load only about a million rows.

Additionally, since I couldn't use the projected data, I am not even sure if the resulting distances are correct. Apart from zero which makes sense for distance between a point (location) and itself, the least distance is around 277. I

I am new to using R, especially to do spatial analysis. I understand I need to create a distance matrix with the buffers around the points(locations). How can I achieve my aim using R?

Answer 1

Not entirely sure if this solves your question. See code below and maybe adjust to your needs (e.g. projection, miles as distance measurement).

library(sf)
library(sp)

data    <- data.frame(id = paste0("ID_", 1:10), lon = runif(10, -10, 10), lat = runif(10, -10, 10), pop_size = runif(10, 0, 2000))
data_sf <- st_as_sf(data,coords = c("lon", "lat"), crs = 4326) %>% st_transform(CRS("+proj=laea"))

# plot(data_sf$geometry)

bufferR <- seq(150000, 450000, length.out = 5) # sequence of radii

out <- do.call("rbind", lapply(1:length(bufferR), function(y) {
    bfr <- data_sf %>% st_buffer(bufferR[y]) ## create Buffer
    ## minus the next smaller buffer
    if(y>1) {
      inters <- suppressWarnings(st_difference(bfr, data_sf %>% st_buffer(bufferR[y-1])))    
      bfr <- inters[which(inters$id == inters$id.1),]
    }

    # get ids that intersect with buffer
    inters <- bfr %>% st_intersects(data_sf) 


    do.call("rbind", lapply(which(sapply(inters, length)>0), 
         function(z) data.frame(orig = data_sf[z,]$id, radius = bufferR[y],
             incl = data_sf[unlist(inters[z]),]$id, 
                pop_size = data_sf[unlist(inters[z]),]$pop_size)))
}))



subset(out, orig == unique(out$orig)[1])

#    orig radius  incl  pop_size
# 1  ID_1 150000  ID_1 1704.4946
# 21 ID_1 225000  ID_4  980.9457
# 45 ID_1 450000 ID_12  530.4216

So in this case, the smallest buffer only includes the point around the buffer (ID_1). With a radius of 225000m, ID_4 is included and ID_12 joins with a radius of 450000 m.