Create grid points and compute distance within regions

Question

The spatial analysis problem is as follows. Suppose you want to compute to average distance a consumer within each postal code has to travel to get the nearest bank. I have coordinates of bank locations (using the sf package, the geometry is point). And I have a shape file with postal codes (this is given as a multipolygon geometry type).

At this stage I want to create a grid with very fine points and overlay it on top of the postal code shape file. Using these points, I need to find the mean distance each point as to travel to get to the nearest bank location. The idea then is the to compute the average by zip code.

I have it conceptually worked out, but I'm stuck on the part where I have to make a fine grid. I'm not even sure how to approach this given I have two difference data sets with different geometries. I'm using the sf package and the tidyverse.

Answer 1

See here a reprex, the key functions are st_make_grid, st_nearest_feature and st_distance. One thing, to make all this work, I need to project the shapes. I used for this example EPSG:3857 for no particular reason, but you have to pay attention to the units on st_crs(yourshape), as the grid would use the same unit.

For the reprex I used countries as zipcodes and cities as bank branches.

library(tidyverse)
library(sf)

#Dummy data
library(rnaturalearth)
zipcode = ne_countries(
  country = c("Germany", "Denmark", "Poland", "Czech Republic", "Austria"),
  returnclass = "sf"
) %>% select(adm0_a3)
banks = ne_download(category = "cultural",
                    type = "populated_places",
                    returnclass = "sf") %>%
  filter(ADM0_A3 %in% zipcode$adm0_a3) %>% select(NAME)


#Need to project shapes, for this example 3857  (m)
zipcode = st_transform(zipcode, 3857)
banks = st_transform(banks, 3857)

st_crs(3857)

# Coordinate Reference System:
#  EPSG: 3857 
#  proj4string: "+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 # # +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"


#Fine grid - adjust as desired, in this case 50,000 m (50kms),
# see st_crs +units
densgrid = 50000

grid = st_make_grid(zipcode, cellsize = densgrid, what = "centers") %>%
  st_intersection(zipcode)

#Get nearest
grid.df = st_sf(indexnearest = st_nearest_feature(grid, banks),
                geom = grid)
#Get distances
disttopoint = st_distance(grid.df, banks, sparse = FALSE)

mindist = lapply(1:nrow(distances), function(x)
  distances[x, grid.df[x,]$indexnearest]) %>% unlist()


grid.df$disttopoint = disttopoint
banks$indexnearest = 1:nrow(banks)

#Bring zipcode, in example is NAME
grid.df = grid.df %>% left_join(st_drop_geometry(banks) %>%
                                  select(indexnearest, NAME))

grid.df %>% st_drop_geometry() %>%
  group_by(NAME) %>%
  summarise(meandist = mean(disttopoint))

# A tibble: 5 x 2
#  NAME       meandist
#  <chr>           [m]
# 1 Berlin     756496.0
# 2 København 1015690.8
# 3 Prague     778736.7
# 4 Vienna     777320.7
# 5 Warsaw     789605.3

#For coloring
nearestbank = factor(grid.df$indexnearest)

ggplot() +
  geom_sf(data = zipcode) +
  geom_point(aes(colour = nearestbank),
             x = st_coordinates(grid.df)[, 1]
             ,
             y = st_coordinates(grid.df)[, 2]) +
  geom_sf(data = banks, col = "red")