2

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.

3

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

enter image description here

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.