In kmcd39/divM: What the Package Does (One Line, Title Case)

knitr::opts_chunk$set(echo = T
                      ,message = F
                      ,warning = F
                      ,dpi=300
                      ,fig.width=6.5
                      ,fig.height=6)
library(tidyverse)
library(sf)

options(tigris_use_cache = TRUE)

rm(list = ls())

library(mapview)

Setup

#' notes:
#'
#' Also remember source for additional metadata:
#' https://www.protectedlands.net/help/
#'
#' The "Technical How-Tos" in that site has "Detailed Descriptions of Layer
#' Files" which is helpful, as well as "Tips for GIS Users"; in particular ---
#'

# CT analysis -------------------------------------------------------------

## load generated measures -------------------------------------------------
# FIPS for selected states
#selfps <- c('09', '25', '44')
selfps <- '09'
read.dir <- '~/R/local-data/usgs/generated-PAD-measurues/'
readpth <- paste0(read.dir, 'statefp-', selfps, '.csv')
dir.exists(read.dir)
padm <- vroom::vroom(readpth)


## also load PAD data ------------------------------------------------------

ddir <- '~/R/local-data/usgs/PADUS3_0Geodatabase/'
gdb <- ddir %>% list.files(pattern = 'gdb$', full.names = T)

lyrs <- gdb %>% st_layers()
lyrs


fcts <- tigris::tracts( year = 2021
                       ,state = '09'
                       ) %>%
  rename_with(tolower) %>%
  st_transform(5070)

# see https://github.com/pauldzy/USGS_Albers_Equal_Area_Projections for EPSG.

wktf <- fcts %>% st_bbox() %>% st_as_sfc() %>% st_as_text()

# read
pad <- st_read(
  gdb
  , layer = lyrs$name[9]
  , wkt_filter = wktf # over sample area
)
# i get a warning/GDAL error reading..

# cast to polygons (get rid of multisurface), then turn to tibble for
# non-spatial processing
pad <- pad %>%
  rename(geometry = SHAPE) %>%
  rename_with(tolower) %>%
  st_cast("MULTIPOLYGON") %>%
  tibble()

## bring in metadata -------------------------------------------------------

# map column name in the data to metadata name; can skip 8, which is just state
# names. Note Agency name and type match either Owner or Manager in the spatial
# data.
metadata2colm <- tibble(
  metadata.nm = lyrs$name[1:7]
  ,colm.nm = Hmisc::Cs(pub_access, featclass, des_tp, gap_sts, iucn_cat, agency_name, agency_type)
)

# create an organized list of dataframes for every metadata layer
metadata <- lyrs$name[1:7] %>%
  map2( metadata2colm$colm.nm
        , ~{st_read(gdb, layer = .x) %>%
            select(!!.y := 1, !!.x := 2)
        }
  ) %>%
  set_names(lyrs$name[1:7])

metadata$Agency_Name
metadata$Agency_Type

# setup for some visuals --------------------------------------------------

# epsg for good CT projection
ctcrs <- 2234
ctcrs %>% st_crs()

# get geos and tract attrs for CT --------------------------------------------------

# get state
ctstate <- tigris::states(year = 2021) %>%
  rename_with(  tolower  ) %>%
  filter(stusps == 'CT') %>%
  st_transform(ctcrs)

# ct counties
ctcos <- tigris::counties(state = 'CT'
                          ,year = 2021
                          ) %>%
  rename_with(  tolower  ) %>%
  st_transform(ctcrs)

ctcos

cattrs <- censusrx::get.tract.attrs(
                      state = '09'
                     ,cofps = ctcos$countyfp
                     ,year = 2021
                     ,geo = 'tract'
                   )

cattrs

ctsf <- tigris::tracts(state = 'CT'
                       ,year = 2021) %>%
  rename_with(tolower) %>%
  st_transform(ctcrs) %>%
  filter(aland > 0)

## also create a subset of raw PAD data for a smaller area ----------------------

# ctcos %>% mapview(zcol = 'geoid')

# Hartford box
selcos <- ctcos %>%
  filter(geoid == '09003') %>%
  st_transform(5070)

# cropped pad
cpad <- pad %>%
  st_sf() %>%
  st_crop(selcos) %>%
  tibble()

Initial Demographic Maps

Population density and median household income, for some preliminary context:

# med hh inc and pop dens
cattrs %>%
  left_join(ctsf)%>%
  st_sf() %>%
  ggplot() +
  geom_sf(aes(fill = med.hhinc)
          ,color = NA) +
  scale_fill_viridis_b(
    name = 'median household income'
    ,labels = ~paste0('$',scales::comma(.x) )
  ) +
  theme_void()


cattrs %>%
  left_join(ctsf)%>%
  st_sf() %>%
  mutate(pop_dens.capped =
           visaux::cap.at.quantile(pop_dens ,.9)
         ) %>%
  ggplot() +
  geom_sf(aes(fill = as.numeric(pop_dens.capped))
          ,color = NA) +
  scale_fill_viridis_b(
     name = 'Pop density\n(persons/acre)'
    ,labels = visaux::ggcapped.labels
  ) +
   theme_void() +
  labs(caption = 'scale capped at 90th percential of census tracts in connecticut (CTs in CT...)')

Subseting protected areas data

Features of the Protected Areas Data (PAD) will overlap and intersect. In general, this will happen by "feature class," which describes how the protected areas are defined. The documentation helps clarify what these are.

Based on looking at the data, it seems a set of filters based on these feature classes are appropriate for trimming the data before any analysis. In particular, we should focus on Fees and Easements, and remove Designations and Proclamations. The latter two classes will frequently overlap with the former two. Fees and easements seem to correspond very well with legal or ownership boundaries, while Designations and Proclamations are more "contextual" or illustrative. Marine features are mostly water, but also comprise protected wetlands and similar, so I think we should include those.

# exploratory visuals of PAD measures & data -----------------------------------------

#' join to aland to calculate n protected acres. remember units of ALAND are sq
#' meters
#'
#' acres column will be protected acreage by type
padm <- padm %>%
  left_join(
    tibble(ctsf)[,c('geoid', 'aland')]
    ) %>%
  mutate(acres =
           tract.perc.area *
           units::set_units(
             units::set_units(aland, 'meters^2')
             ,'acres'           )
  )

padm$acres <- as.numeric(padm$acres)

## analysis by land managers -----------------------------------------------

mgmt <- padm %>%
  filter(grepl('^mang',
               protected.area.descriptor)) %>%
  group_by(protected.area.descriptor,
           protected.area.value) %>%
  summarise(count = n()
            ,acres = sum(acres)) %>%
  arrange(protected.area.descriptor, desc(acres)) %>%
  ungroup()

# add long name
metamanagers <-
  metadata[grepl('^Agency', names(metadata))] %>%
  imap( ~rename(.x, protected.area.value = 1,
                label = 2 )) %>%
  imap_dfr( ~mutate(.x,
                    protected.area.descriptor =
                      paste0('mang_',
                             str_extract( tolower(.y)
                                         ,'name$|type$'
                                         ))
                    ,.before = everything()
                    )
            )
# mgmt %>%
#   left_join(metamanagers) %>%
#   # filter( protected.area.descriptor ==
#   #           'mang_type') %>%
#   mutate( label =
#             if_else( row_number() >= 7
#                      ,'Unknown or other'
#                      ,label
#             ))  %>% View()


mgmt <- mgmt %>%
  left_join(metamanagers) %>%
  filter( protected.area.descriptor ==
            'mang_type') %>%
  group_by(protected.area.descriptor) %>%
  mutate( label =
            if_else( row_number() >= 7
                              ,'Unknown or other (recode)'
                              ,label
                     )) %>%
  group_by(label) %>%
  summarise(across(c(count, acres)
                   ,~sum(.x, na.rm = T) # NAs are 0s there
                   )) %>%
  arrange(desc(acres)) %>%
  mutate(label = factor(label
                        ,levels = .$label))


### map by those categories -------------------------------------------------


# tract-level by mgmt type
ct.mgmt <- padm %>%
  filter( protected.area.descriptor ==
            'mang_type') %>%
  left_join(metamanagers)
# %>%
#   mutate(label =
#            if_else(label %in% mgmt$label
#                    , label
#                    , 'Unknown or other (recode)'))

# (i haven't re-aggregated to recode yet.. because there are overlaps.)


# facet plot by mgmr type
ct.mgmt %>%
  full_join(ctsf[c('geoid', 'geometry')]) %>%
  st_sf() %>%
  ggplot() +
  geom_sf(data = ctstate
          ,fill = 'grey70') +
  geom_sf( aes(fill = tract.perc.area)
           ,color = NA
           ,linewidth = 0
           ) +
    scale_fill_viridis_c() +
  facet_wrap(vars(label)) +
  theme_void() +
  theme(strip.background =
          element_rect(fill = visaux::jewel.pal()[4] )
        ,strip.text = element_text(
          color = 'white', face = 'bold'
        )) +
  labs(
    subtitle = 
      str_wrap('We can see how the Federal protected areas covers so much. This is due to a Proclamation, one of the feature classes by which we should filter.'
  )) 

# pad %>%
#   filter(mang_type == 'FED') %>%
#   st_sf() %>%
#   mapview::mapview()


#' what comprise the designations by mgmt type? -- for the non-Federal common
#' mgmt types.
#'
#' Federal managed areas dwarf those w all other managers, and are entirely
#' "outer continenal shelf areas" (water) or "Approved or Proclamation Boundary"
#'
#' Note metadata on those proclamation boundaries (from
#' https://www.protectedlands.net/pad-us-technical-how-tos/#feature-classes-in-pad-us-2)
#' --
#'
#' Proclamations – Boundaries of the administrative area of National Forests,
#' Parks, Wildlife Refuges and other lands. These are not ownership lines but
#' are used for agency administration purposes (note that some commercial
#' mapping providers and others incorrectly use these boundaries to show
#' protected areas and in doing so often show large areas of private lands as
#' part of public lands). Some proclamation boundaries can cover extremely large
#' areas of private land, while others are close to ownership boundaries. In any
#' case, it’s important not to consider proclamation boundaries as reflecting
#' actual land ownership – they’re just for agency administration purposes


tmp <-
  pad %>%
  filter(mang_type %in%
           filter(metamanagers #
                  ,protected.area.descriptor == 'mang_type' &
                    label %in% mgmt$label[2:6])$protected.area.value
                    #label %in% mgmt$label[1:6])$protected.area.value

           ) %>%
  group_by(mang_type,
           des_tp) %>%
  summarise( count = n()
            ,acres = sum(gis_acres)) %>%
  arrange(mang_type, desc(acres)) %>%
  left_join(
    rename(metadata$Designation_Type
           ,label = 2)
            ) %>%
  ungroup()

ordr <- tmp %>%
  group_by(des_tp, label) %>%
  summarise(acres = sum(acres)) %>%
  arrange(desc(acres))

tmp$label <- factor(tmp$label
                    ,levels = ordr$label)

tmp %>%
  left_join(metadata$Agency_Type
            ,by = c('mang_type' = 'agency_type')
            ) %>%
  ggplot(
    aes(y = label
        ,x = acres
        ,fill = Agency_Type)
  ) +
  geom_col() +
  scale_y_discrete(
    name = NULL
    ,limits = rev
  ) +
  scale_fill_manual(
    name = 'Managing entity'
    ,values = visaux::jewel.pal()
  ) +
  theme_minimal() +
  visaux::upper.legend.box() +
  labs(
    title = 'Acreage of Protected Lands by Purpose and Managing Entity in CT'
    ,subtitle = str_wrap('With Federal mgmt removed')
  )

#' "Proclamations" with federal managers dwarf all other types of protected
#' lands, but are omitted from the plot -- because proclamations have very
#' different interpretaions.


# analysis by designation ------------------------------------------------

#' recode.designations
#'
#' Function to recode labels for des_tp to some combined categories. Assumes a
#' label column from metadata.
#'
#'
recode.designations <- function(x) {

  x %>%
    mutate(recode =
             case_when(
               grepl('Park|Recreation', label) ~ 'Park or recreation'
               ,grepl(
                 'Watershed Protection|Resource Management|Forest Stewardship|Wildlife Refuge|National Forest'
                 , label) ~
                 str_wrap(width = 45,'Watershed protection, resource management, or forest stewardship')
               ,grepl('Conservation', label) ~ 'Conservation'
               ,grepl('Military|Research', label) ~ 'Military or research area'
               ,grepl('Native American', label) ~ 'Native American land area'
               ,grepl('Recreation|Cultural|Historic', label) ~ 'Other recreational/historic/cultural'
               ,grepl('Agricultural|Ranch', label) ~ 'Agricultural or ranch'
               ,grepl('Marine|Continental Shelf', label) ~ 'Marine or continental shelf area'
               ,TRUE ~ 'Unkown or Other'
             ) )

}

padm %>%
  filter(protected.area.descriptor ==
           'des_tp') %>%
  left_join(
    rename(metadata$Designation_Type,
           protected.area.value = 1, label = 2)
  ) %>%
  recode.designations()

pad$gis_acres %>% quantile(seq(0, 1, .1))


tmp <-
  pad %>%
    filter(! featclass %in% c( 'Proclamation'
                              ,'Marine')) %>%
  left_join(
    rename(metadata$Designation_Type,
           des_tp = 1, label = 2)
  ) %>%
  recode.designations()

pad %>%
  group_by(featclass) %>%
  summarise(count = n(),
            acres = sum(gis_acres)) %>%
  arrange(desc(acres))

pad %>%
  group_by(featclass) %>%
  taux::quantiles.across.groups('gis_acres')

pad$gis_acres %>% quantile(seq(0,1,.1))



# exploratory maps of subarea ---------------------------------------------


#' it makes sense to really just retain Easements and Fees.
#'
#' Designations tend to overlap with the other cat
tmp <- cpad %>%
  filter(! featclass %in% c( 'Proclamation'
                             ,'Designation')) %>%
  filter( des_tp != 'OCS' # (outer continental shelf)
          ) %>%
  left_join(
    rename(metadata$Designation_Type,
           des_tp = 1, label = 2)
  ) %>%
  recode.designations()

tmp$gis_acres %>% quantile(seq(0,1,.1))

#tmp %>% st_sf() %>% mapview(zcol = 'featclass')


# making sure i understood Marine correctly (it also contains some swamps etc.)
pad %>%
  filter( featclass %in% c( 'Marine'
                            )) %>%
  left_join(
    rename(metadata$Designation_Type,
           des_tp = 1, label = 2)
  ) %>%
  recode.designations() %>%
  st_sf() %>%
  mapview(zcol = 'mang_name')



# ct analysis with revised types ------------------------------------------

trimmed.pad <- pad %>%
  filter(! featclass %in% c( 'Proclamation'
                             ,'Designation') &
           des_tp != 'OCS' # (outer continental shelf)
         )

# check overlaps from within this subset...
#overlaps <- st_intersection(trimmed.pad$geometry)

devtools::load_all()

padm <- trimmed.pad %>%
  protected.area.by.type.by.tract(
        cts = ctsf
        ,category.colm = 'des_tp'
      )


# with labels and recodes
padmr <- padm %>%
  left_join(
    rename(metadata$Designation_Type,
           des_tp = 1, label = 2)
  ) %>%
  recode.designations()

padmr %>% select(des_tp, label, recode) %>% distinct() %>% arrange(recode)



# add acres
padmra <- padmr %>%
  left_join(
    tibble(ctsf)[,c('geoid', 'aland')]
  ) %>%
  mutate(acres =
           perc.area *
           units::set_units(
             units::set_units(aland, 'meters^2')
             ,'acres'           )
  ) %>%
  group_by(geoid, recode) %>%
  summarise(across(c(perc.area, acres)
                   ,sum
                   ))

padmra$perc.area %>% quantile(seq(0,1,.05))

# plurality map of designation.
padmra %>%
  full_join(ctsf[c('geoid', 'geometry')]) %>%
  st_sf() %>%
  mutate(capped_perc.area =
           visaux::cap.at.quantile(perc.area
                                   ,.9)
         ) %>%
  ggplot() +
  geom_sf(data = ctstate
          ,fill = 'grey90') +
  geom_sf( aes(fill = capped_perc.area)
           ,color = NA
           ,linewidth = 0
  ) +
  scale_fill_viridis_c(
     '% of tract'
    ,labels = visaux::ggcapped.labels
  ) +
  facet_wrap(vars(recode)) +
  theme_void() +
  theme(strip.background =
          element_rect(fill = visaux::jewel.pal()[4] )
        ,strip.text = element_text(
          color = 'white', face = 'bold'
        )) +
  visaux::upper.legend.box()

# selected designations within that filters.. just drop agriculatural/ranch,
# native american lands, and military/research lands
padmr %>% select(des_tp, label, recode) %>% distinct() %>% arrange(recode) %>% View()

padsel <- padmr %>%
  filter( !grepl('Agricultural|Military|Native American Land Area',
                 label) )

padsel <- padsel %>%
  left_join(
    tibble(ctsf)[,c('geoid', 'aland')]
  ) %>%
  mutate(acres =
           perc.area *
           units::set_units(
             units::set_units(aland, 'meters^2')
             ,'acres'           )
  ) %>%
  group_by(geoid) %>%
  summarise(across(c(perc.area, acres)
                   ,sum
  ))


padsel %>% select(perc.area, acres) %>% map( ~quantile(.x, seq(0,1,.1)))



padsel %>%
  full_join(ctsf[c('geoid', 'geometry')]) %>%
  st_sf() %>%
  mutate(capped_perc.area =
           visaux::cap.at.quantile(perc.area
                                   ,.9)
  ) %>%
  ggplot() +
  geom_sf(data = ctstate
          ,fill = 'grey90') +
  geom_sf( aes(fill = capped_perc.area)
           ,color = NA
           ,linewidth = 0
  ) +
  scale_fill_viridis_c(
    '% of tract'
    ,labels = visaux::ggcapped.labels
  ) +
  theme_void() +
  theme(strip.background =
          element_rect(fill = visaux::jewel.pal()[4] )
        ,strip.text = element_text(
          color = 'white', face = 'bold'
        )) +
  visaux::upper.legend.box()+
  labs(
    title = 'Protected area coverage for selected types'
    ,subtitle =
      str_wrap(width = 80,
'Filters are: Not Proclamation or Designation feature class (only Fees, Easements, and Marine); not continental shelf; not Indian land; not agricultural or military.')
  )


# bivariate plots ------------------------------------------------------

#install.packages("biscale", dependencies = TRUE)
library(biscale)
library(cowplot)

##  w pop dens and protected land -----------------------------------------

tmp <- padsel %>%
  full_join(cattrs)

tmp %>%
  select(geoid, pop, pop_dens,
         aland.acre, perc.area, med.hhinc) %>%
  filter(is.na(perc.area) |
           is.na(med.hhinc))

tmp$pop_dens %>% as.numeric()

btmp <-
  tmp %>%
  mutate(pop_dens = as.numeric(pop_dens)
         ) %>%
  filter(!is.nan(pop_dens)
         ) %>%
  bi_class( x =
              pop_dens
           , y = perc.area
           , style = "quantile",
           dim = 3)


# bi_pal('DkCyan')
bi.map <-
  btmp %>%
  full_join(ctsf[c('geoid', 'geometry')]) %>%
  st_sf() %>%
  ggplot() +
  geom_sf(aes(fill = bi_class),
          color = NA,
          size = 0.0,
          show.legend = FALSE) +
  bi_scale_fill(pal = "DkCyan",
                dim = 3) +
  bi_theme()

#?bi_pal()
bi.lgnd <-
  bi_legend(pal = "DkCyan",
            dim = 3,
            xlab = "More population ",
            ylab = "More protected area",
            size = 8)

cowplot::ggdraw() +
  draw_plot(bi.map, 0, 0, 1, 1) +
  draw_plot(bi.lgnd, .7, .1, 0.25, 0.25)

##  w hh income and protected land -----------------------------------------




tmp <- padsel %>%
  full_join(cattrs)

tmp %>%
  select(geoid, pop, aland.acre, perc.area, med.hhinc) %>%
  filter(is.na(perc.area) |
           is.na(med.hhinc))

cattrs %>%
  select(geoid, pop, n.hh, aland.acre, med.hhinc) %>%
  filter(is.na(med.hhinc))

cattrs$n.hh
btmp <-
  bi_class(tmp
           , x =
             med.hhinc
           , y = perc.area
           , style = "quantile", dim = 3)


bi_pal('PurpleGrn')
bi.map <-
  btmp %>%
  full_join(ctsf[c('geoid', 'geometry')]) %>%
  st_sf() %>%
  ggplot() +
  geom_sf(aes(fill = bi_class),
          color = NA,
          size = 0.0,
          show.legend = FALSE) +
  bi_scale_fill(pal = "PurpleGrn",
                dim = 3) +
  bi_theme()

bi.lgnd <-
  bi_legend(pal = "PurpleGrn",
            dim = 3,
            xlab = "Higher income ",
            ylab = "More protected area",
            size = 8)

cowplot::ggdraw() +
  draw_plot(bi.map, 0, 0, 1, 1) +
  draw_plot(bi.lgnd, .7, .1, 0.25, 0.25)
#