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(stringr) # for working with strings (pattern matching)
library(rgdal) #For checking available CRSs interactively. 
library(lwgeom) #Needed for distances

#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

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?

  • If you use the 64bit (x64) version of R you can get around the row limitation, I have a feeling that you are using the 32bit version. So, do you actually need the buffers are are you under the assumption that this is how you would get nearest neighbors? You could use the knn function in spatialEco to return the kNN for each point. It will take a minimum/maximum distance argument. You could also follow this approach: gis.stackexchange.com/questions/163287/… – Jeffrey Evans Sep 26 at 17:23

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


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.

| improve this answer | |
  • Thank you @SimeonL for your answer. I think what you suggested is a single buffer around each point. I need to create multi-ring buffers around each point. How do I achieve that and then be able to get the index of points within each multi-ring buffer so that if I have 5 different buffers (eg. 1, 2, 3, 4, 5) I get 5 columns like incl1, incl2, incl3, incl4, and incl5 where points within buffer 2 does not include points within buffer 1, etc? – Lomnewton Apr 14 at 2:04
  • Not sure how you want to store the indices of the points within each buffer (minus the points of the next smaller buffer) in one column per point per buffer. This could be a string of variable lenght? I think it depends on what you want to do with the indices. Maybe a few more words about the application would help. – SimeonL Apr 14 at 6:09
  • I intend to do a number of things. For example I will be estimating a regression model of the form Y_it = beta_1*sum of pop_size of point i's neighbors within a specified radius. So my thinking is that, if I have separate columns for the definition point i's neighbors based on the respective buffers, I could run the model for each of them separate and be able to compare the estimated value of beta_1. I may be wrong on the arrangement of the data for such an analysis. But I hope you get an idea of what I am trying to do. – Lomnewton Apr 14 at 19:27
  • I think it is best to remain a long form with all information about radius/buffer, the ID of interest (orig) and the including id's per buffer (incl). Then you can filter/subset according to your question. I edited the answer. – SimeonL Apr 15 at 8:34
  • Can @SimeonL or anyone else help resolve this error I am getting : Error in CPL_geos_op2(op, st_geometry(x), st_geometry(y)) : Evaluation error: std::bad_alloc. The code runs well when I implement it on round a 1000 location points. But beyond a 1000 points, R throws out the above error. – Lomnewton May 30 at 6:52

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.