Libraries

# install.packages("devtools")
#devtools::install_github("UrbanInstitute/urbnmapr")
#install.packages("remotes")
#remotes::install_github("UrbanInstitute/urbnthemes", build_vignettes = TRUE)
library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──
## ✓ ggplot2 3.3.5     ✓ purrr   0.3.4
## ✓ tibble  3.1.6     ✓ dplyr   1.0.7
## ✓ tidyr   1.1.4     ✓ stringr 1.4.0
## ✓ readr   2.1.1     ✓ forcats 0.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(dplyr)

Read CSV

bm202004 <- read.table(file = "bm2004.csv",
                       sep = ",",
                       header = TRUE,
                       colClasses = "character")

Filter PA VSOs (Subsection 19)

SUBSECCD==“19” is Post or Organization of Past or Present Members of the Armed Forces (Form 990 or 990-EZ or Form 1024) SUBSECCD==“23” is Veterans Organizations Created Before 1880 OUTNCCS = “IN” is NCCS “Out of Scope Flag”, for US nonprofit sector analysis, often foreign-based entity filing with IRS or operating in US territories or overseas NTEECC==“W30” is Military & Veteran Organizations

filter(SUBSECCD==“23”, OUTNCCS==“OUT”) returns 35 rows, largely Puerto Rico, U.S. Virgin Islands, Guam, Germany, Japan but also Springfield, VA, Duncansville, PA, Greybull, WY and APO, AE (American Post Office in Europe, Middle East, Africa or Canada)

filter(SUBSECCD==“23”, OUTNCCS==“OUT”) returns 0 rows

filter(STATE==“PA”, SUBSECCD==“19”, OUTNCCS==“OUT”) returns 1 row, Duncansville, PA VFW ???

filter(STATE==“PA”, NTEECC==“W30”) returns 2156 rows filter(STATE==“PA”, SUBSECCD==“19”) returns 1300 rows filter(STATE==“PA”, SUBSECCD==“23”) returns 0 rows filter(SUBSECCD==“23”) returns 2 rows (Mutual Aid Associations) in Arlington and Reston VA filter(STATE==“PA”, NTEECC==“W30”, SUBSECCD==“19”) returns 1187 rows

pa_vso_bm202004 <-
  bm202004 %>%
    filter(STATE=="PA", NTEECC=="W30")
pa_vso_bm202004

Examine Centre County (FIPS 42027)

centre_pa_vso_bm202004 <-
  pa_vso_bm202004 %>%
    filter(FIPS=="42027.0")
centre_pa_vso_bm202004

Summarize

fips_vso_sum <- 
  bm202004 %>%
    #filter(STATE=="PA", SUBSECCD=="19") %>%
    #filter(SUBSECCD=="19") %>%
    filter(NTEECC=="W30") %>%
    group_by(FIPS) %>%
    tally(sort = TRUE, name = "vsos") %>%
    rename(county_fips = FIPS)

fips_vso_sum <- drop_na(fips_vso_sum)
fips_vso_sum$county_fips <- gsub("\\..*","",fips_vso_sum$county_fips)
fips_vso_sum$county_fips <- sprintf("%05s", fips_vso_sum$county_fips)

Map

library(sf)
## Linking to GEOS 3.10.2, GDAL 3.4.1, PROJ 8.2.1; sf_use_s2() is TRUE
library(urbnmapr)
library(urbnthemes)
## Setting Mac/Linux options...
## 
## Attaching package: 'urbnthemes'
## The following objects are masked from 'package:ggplot2':
## 
##     geom_bar, geom_col, scale_color_discrete, scale_color_gradientn,
##     scale_color_ordinal, scale_colour_discrete, scale_colour_gradientn,
##     scale_colour_ordinal, scale_fill_discrete, scale_fill_gradientn,
##     scale_fill_ordinal
#set_urbn_defaults(style = "map")
counties_sf <- get_urbn_map("counties", sf = TRUE) %>%
#filter(state_name %in% c("Pennsylvania"))  
filter(state_name %in% c("Pennsylvania", "New York", "New Jersey"))
## old-style crs object detected; please recreate object with a recent sf::st_crs()
bm_map_data <- st_as_sf(left_join(fips_vso_sum, counties_sf, by = "county_fips"))

bm_map_data %>%
  ggplot() +
  geom_sf(aes(
    # Color in states by the chip_pct variable
    fill = vsos
  )) + 
  # scale_fill_gradientn adds colors with more interpolation and reverses color scale
  scale_fill_gradientn(
    # Legend format
    labels = scales::comma_format(),
    # Make legend title more readable
    name = "VSO Counts",
    # Manually add 0 to lower limit to include it in legend. NA=use maximum value in data
    limits = c(0, NA),
    # Set number of breaks on legend = 3
    n.breaks = 7
    )

library(mapview)
mapview(bm_map_data, zcol = "vsos")