I. Introduction

[add background and motivations]

..

read in insts

(insts <- read_csv('input_data/temple x - Sheet1.csv', col_select = -lat_long) |> 
   drop_na() |> # can't have missing coords
   rowid_to_column(var = 'id_tx') |> 
   st_as_sf(coords = c('long', 'lat'), crs = st_crs(4326)) |> 
   st_transform(crs = st_crs(patches)) |> 
   select(starts_with('id'), name, type, address) # cosmetic re-order
 )
## Rows: 361 Columns: 6
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (3): type, name, address
## dbl (3): id_orig, lat, long
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
## drop_na: removed 2 rows (1%), 359 rows remaining
## 
## select: columns reordered (id_tx, id_orig, name, type, address, …)
## Simple feature collection with 359 features and 5 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 1369280 ymin: 567927.2 xmax: 1448926 ymax: 625099.2
## Projected CRS: NAD83 / Maryland (ftUS)
## # A tibble: 359 × 6
##    id_tx id_orig name                    type  address                 geometry
##    <int>   <dbl> <chr>                   <chr> <chr>   <POINT [US_survey_foot]>
##  1     1       1 New Covenant Church     fait… 110 Re…       (1393251 618347.5)
##  2     2       2 Love Alive Ministries   fait… 5105 G…       (1394979 591793.4)
##  3     3       3 The Broken Walls Proje… fait… 181 N …       (1395922 588917.1)
##  4     4       4 Concord Baptist         fait… 5204 L…       (1396611 607310.3)
##  5     5       5 Miracle Seventh Day Ad… fait… 100 S …       (1397077 588739.1)
##  6     6       6 StillMeadow Community … fait… 5110 F…       (1397357 588260.4)
##  7     7       7 New Miracle CCC-48      fait… 4802 L…       (1397947 606686.7)
##  8     8       8 Bible Way Christian Fe… fait… 4412 M…       (1399121 605782.8)
##  9     9       9 Park Community United … fait… 3604 M…       (1399168 606443.5)
## 10    10      10 Parklane Baptist Church fait… 3606 M…       (1399171 606591.6)
## # … with 349 more rows
# 2x checks
# FIXME - what to do about spatial duplicates? Throw an iterative distinct on there?
(insts_dups_ids <- 
  insts |> 
  st_drop_geometry() |> 
  group_by(id_orig) |> 
  count() |> 
  filter(n > 1) |> 
  arrange(desc(n)))
## group_by: one grouping variable (id_orig)
## count: now 210 rows and 2 columns, one group variable remaining (id_orig)
## filter (grouped): removed 148 rows (70%), 62 rows remaining
## # A tibble: 62 × 2
## # Groups:   id_orig [62]
##    id_orig     n
##      <dbl> <int>
##  1      51    11
##  2     406     9
##  3      21     8
##  4     142     8
##  5     164     8
##  6      67     7
##  7     217     7
##  8     239     7
##  9     348     7
## 10      81     6
## # … with 52 more rows
(insts_dups_adds <- 
  insts |> 
  st_drop_geometry() |> 
  group_by(address) |> 
  count() |> 
  filter(n > 1) |> 
  arrange(desc(n)))
## group_by: one grouping variable (address)
## count: now 273 rows and 2 columns, one group variable remaining (address)
## filter (grouped): removed 237 rows (87%), 36 rows remaining
## # A tibble: 36 × 2
## # Groups:   address [36]
##    address                                     n
##    <chr>                                   <int>
##  1 3400 Ellerslie Ave, Baltimore, MD 21218     9
##  2 3701 Eldorado Ave, Baltimore, MD 21207      8
##  3 2400 Round Rd, Baltimore, MD 21225          7
##  4 2400 Windsor Ave, Baltimore, MD 21216       7
##  5 1201 Cambria St, Baltimore, MD 21225        6
##  6 128 W Franklin St, Baltimore, MD 21201      6
##  7 1406 N Ellamont St, Baltimore, MD 21216     6
##  8 3301 Carlisle Ave, Baltimore, MD 21216      6
##  9 1900 Edgewood St, Baltimore, MD 21216       5
## 10 100 Kane St, Baltimore MD 21224             4
## # … with 26 more rows
# take a look, on the map
# insts |> mapview(zcol = 'type') # TODO exclude those outside of the city?

map of patches and institutions

# combo map
mapview(patches, col.regions = 'darkgreen') + 
  mapview(insts, zcol = 'type')

how many faith based insts overlap with pathces?

insts |> st_intersection(patches) # NONE!
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
## Simple feature collection with 0 features and 8 fields
## Bounding box:  xmin: NA ymin: NA xmax: NA ymax: NA
## Projected CRS: NAD83 / Maryland (ftUS)
## # A tibble: 0 × 9
## # … with 9 variables: id_tx <int>, id_orig <dbl>, name <chr>, type <chr>,
## #   address <chr>, Id <dbl>, gridcode <dbl>, type.1 <chr>,
## #   geometry <GEOMETRY [US_survey_foot]>
# within 1000 feet?
insts_1000_feet <- 
  insts |> 
  st_buffer(1000) %>%
  mutate(area = st_area(.)) |> 
  st_intersection(patches) |> 
  group_by(id_tx) |> count() 
## mutate: new variable 'area' (double) with 172 unique values and 0% NA
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
## group_by: one grouping variable (id_tx)
## count: now 195 rows and 3 columns, ungrouped
((insts |> nrow()) - (insts_1000_feet |> nrow())) # 164!
## [1] 164
((insts_1000_feet |> nrow()) / (insts |> nrow())) * 100 # 54%
## [1] 54.31755
# combo map with buffers - cartoonish, ew
mapview(patches, col.regions = 'darkgreen') + 
  mapview(insts |> st_buffer(1000), zcol = 'type')

make isochrones

walk_times <- seq(5, 15, 5)

# # create walksheds
# tic(); insts_iso <-
#   insts |>
#   mb_isochrone(
#       profile = 'walking'
#     , time = walk_times # 5, for testing 
#     , id_column = 'id_tx'
#   ) |>
#   st_transform(st_crs(patches)) |>
#   rename(id_tx = id); toc(); beepr::beep() # ~90 seconds with 5 min walk
# 
# # write out
# insts_iso |>
#   st_write(paste0('input_data/walksheds_5_10_15/walksheds_5_10_15_', Sys.Date(), '.gpkg'))

# read in
insts_iso <- 
  st_read('input_data/walksheds_5_10_15/walksheds_5_10_15_2023-01-17.gpkg')
## Reading layer `walksheds_5_10_15_2023-01-17' from data source 
##   `/Users/dlocke/patches/BGS/temple_x/input_data/walksheds_5_10_15/walksheds_5_10_15_2023-01-17.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 1077 features and 2 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 1365852 ymin: 564809.3 xmax: 1452914 ymax: 627304.9
## Projected CRS: NAD83 / Maryland (ftUS)
mapview(patches, col.regions = 'darkgreen') + 
  mapview(insts_iso, alpha.regions = .1, zcol = 'time') + 
  mapview(insts, cex = .75)
insts_iso
## Simple feature collection with 1077 features and 2 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 1365852 ymin: 564809.3 xmax: 1452914 ymax: 627304.9
## Projected CRS: NAD83 / Maryland (ftUS)
## First 10 features:
##    time id_tx                           geom
## 1    15     1 POLYGON ((1391262 621711.7,...
## 2    10     1 POLYGON ((1392679 620569.1,...
## 3     5     1 POLYGON ((1393248 619459.6,...
## 4    15     2 POLYGON ((1397514 595452, 1...
## 5    10     2 POLYGON ((1396669 594248, 1...
## 6     5     2 POLYGON ((1395541 593083.9,...
## 7    15     3 POLYGON ((1397046 591280.4,...
## 8    10     3 POLYGON ((1394785 590414.6,...
## 9     5     3 POLYGON ((1396769 589370.2,...
## 10   15     4 POLYGON ((1398580 610723.7,...

create walkability surface - argh

# per five minute walk, which ones overlap with patches?
(insts_iso_5 <- 
  insts_iso |> 
  filter(time == 5) |> 
  st_intersection(patches) |> 
  st_drop_geometry() |> 
  group_by(id_tx) |> 
  count())
## filter: removed 718 rows (67%), 359 rows remaining
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
## group_by: one grouping variable (id_tx)
## count: now 205 rows and 2 columns, one group variable remaining (id_tx)
## # A tibble: 205 × 2
## # Groups:   id_tx [205]
##    id_tx     n
##    <int> <int>
##  1     3     8
##  2     4     3
##  3     5     5
##  4     6    10
##  5     7     2
##  6    12     1
##  7    14     2
##  8    17     4
##  9    18     3
## 10    19     3
## # … with 195 more rows
insts_iso_5
## # A tibble: 205 × 2
## # Groups:   id_tx [205]
##    id_tx     n
##    <int> <int>
##  1     3     8
##  2     4     3
##  3     5     5
##  4     6    10
##  5     7     2
##  6    12     1
##  7    14     2
##  8    17     4
##  9    18     3
## 10    19     3
## # … with 195 more rows
((insts |> nrow()) - (insts_iso_5 |> nrow())) # 154!
## [1] 154
((insts_iso_5 |> nrow()) / (insts |> nrow())) * 100 # 57%
## [1] 57.10306
ten_cols <- grDevices::hcl.colors(10, palette = "viridis")

# join the counts back to polygons (either locations or isos) and map with patches
mapview(patches, col.regions = 'darkgreen') +
# insts |> 
insts_iso |> 
  filter(time == 5) |> 
  left_join(insts_iso_5, by = 'id_tx') |> 
  # mutate(n = ifelse(is.na(n), 0, n)) |> 
  # filter(n > 0) |> 
  mapview(zcol = 'n', alpha.regions = .2
          , color.regions = ten_cols
          )
## filter: removed 718 rows (67%), 359 rows remaining
## left_join: added one column (n)
##            > rows only in x   154
##            > rows only in y  (  0)
##            > matched rows     205
##            >                 =====
##            > rows total       359
(insts_iso_time <- 
  insts_iso |> 
  # filter(time == 5) |> 
  st_intersection(patches) |> 
  st_drop_geometry() |> 
  group_by(id_tx, time) |> 
  count())
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
## group_by: 2 grouping variables (id_tx, time)
## count: now 863 rows and 3 columns, 2 group variables remaining (id_tx, time)
## # A tibble: 863 × 3
## # Groups:   id_tx, time [863]
##    id_tx  time     n
##    <int> <int> <int>
##  1     1    10     1
##  2     1    15     5
##  3     2    10     6
##  4     2    15    22
##  5     3     5     8
##  6     3    10    18
##  7     3    15    35
##  8     4     5     3
##  9     4    10    16
## 10     4    15    27
## # … with 853 more rows
((insts |> nrow()) - (insts_iso_time |> filter(time == 10) |> nrow())) # 47!
## filter (grouped): removed 551 rows (64%), 312 rows remaining
## [1] 47
((insts_iso_time |> filter(time == 10) |> nrow()) / (insts |> nrow())) * 100 # 86%
## filter (grouped): removed 551 rows (64%), 312 rows remaining
## [1] 86.90808
((insts |> nrow()) - (insts_iso_time |> filter(time == 15) |> nrow())) # 13!
## filter (grouped): removed 517 rows (60%), 346 rows remaining
## [1] 13
((insts_iso_time |> filter(time == 15) |> nrow()) / (insts |> nrow())) * 100 # 96%
## filter (grouped): removed 517 rows (60%), 346 rows remaining
## [1] 96.37883

Analyses prepared by Dexter H. Locke, PhD USDA Forest Service

Questions?

compiled on 2023-01-18 07:44:44