# Create grid points and compute distance within regions

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`.

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"
type = "populated_places",
returnclass = "sf") %>%

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