DV Lab 7.1

Arizona Hospitals

Author

Aaron Devine

library(ggplot2)
library(ggthemes)
library(socviz)
library(maps)
library(mapproj)
library(viridis)

GET ARIZONA POPULATION DATA BY COUNTY DIRECTLY FROM THE US CENSUS VIA API

head(as_tibble(az_pop))
# A tibble: 6 × 6
  GEOID NAME                   variable estimate   moe                  geometry
  <chr> <chr>                  <chr>       <dbl> <dbl>        <MULTIPOLYGON [°]>
1 04012 La Paz County, Arizona B01003_…    16681    NA (((-114.7312 33.30404, -…
2 04001 Apache County, Arizona B01003_…    66054    NA (((-110.0007 36.99798, -…
3 04013 Maricopa County, Ariz… B01003_…  4430871    NA (((-113.3351 33.55195, -…
4 04015 Mohave County, Arizona B01003_…   214229    NA (((-114.7532 36.08951, -…
5 04005 Coconino County, Ariz… B01003_…   144705    NA (((-113.3541 36.04209, -…
6 04007 Gila County, Arizona   B01003_…    53419    NA (((-111.7207 34.1627, -1…

Here’s a look at our Arizona County data.

CREATE A CHOLOPLETH USING LEAFLET FOR ARIZONA COUNTIES

MapPalette <- colorQuantile(palette = "inferno", domain = az_pop$estimate, n = 20)

library(sf)
az_pop %>% 
  st_transform(crs = "+proj=longlat +datum=WGS84") %>% 
  leaflet(width = "100%", height = 500) %>% 
  addProviderTiles(provider = "Esri.WorldPhysical") %>% 
  addPolygons(popup = ~NAME,
              stroke = FALSE,
              smoothFactor = 0,
              fillOpacity = 0.7,
              color = ~ MapPalette(estimate)) %>% 
  addLegend("bottomright", 
            pal = MapPalette,
            values = ~ estimate,
            title = "Population Percentiles",
            opacity = 1) 

Here are the population percentiles plotted by county using the ‘inferno’ palette from Viridis. Felt like this palette matched the Arizona desert.

library(readxl)
Hospitals <- read_excel("C:/Users/legol/OneDrive - John Brown University/Data Viz/Hospitals.xlsx")

HospitalsAZ <- Hospitals %>% filter(STATE=="AZ")

HospitalsAZ %>% leaflet(width = "100%") %>% 
             addTiles() %>% 
             setView(-111.0937, 34.0489, zoom = 6) %>% 
             addMarkers(lat = ~Y, 
                                 lng = ~X, 
                                 popup = HospitalsAZ$NAME)

Here are Arizona hospitals with popups displaying the hospital name.

MapPalette <- colorQuantile(palette = "inferno", domain = az_pop$estimate, n = 20)


az_pop %>% 
  st_transform(crs = "+proj=longlat +datum=WGS84") %>% 
  leaflet(width = "100%", height = 500) %>% 
  addProviderTiles(provider = "Esri.WorldPhysical") %>% 
  addPolygons(popup = ~NAME,
              stroke = FALSE,
              smoothFactor = 0,
              fillOpacity = 0.7,
              color = ~ MapPalette(estimate)) %>% 
  addLegend("bottomright", 
            pal = MapPalette,
            values = ~ estimate,
            title = "Population Percentiles",
            opacity = 1) %>% 
addCircleMarkers(data = HospitalsAZ, 
                   lat = HospitalsAZ$Y,
                   lng = HospitalsAZ$X,
                   popup = HospitalsAZ$NAME,
                   weight = 1,
                   radius=4,
                   color = "blue", 
                   opacity = 1)

AZ hospitals overlaying the Choropleth. We see a huge concentration in the Phoneiz area which is to be expected.

az_pop %>% 
  st_transform(crs = "+proj=longlat +datum=WGS84") %>% 
  leaflet(width = "100%", height = 500) %>% 
  addProviderTiles(provider = "OpenStreetMap") %>% 
  addPolygons(popup = ~NAME,
              stroke = FALSE,
              smoothFactor = 0,
              fillOpacity = 0.7,
              color = ~MapPalette(estimate)) %>% 
  addLegend("bottomright", 
            pal = MapPalette,
            values = ~estimate,
            title = "Population Percentiles",
            opacity = 1) %>% 
  addCircleMarkers(data = HospitalsAZ, 
                   lat = ~Y,
                   lng = ~X,
                   popup = ~paste("Hospital:", NAME, "<br>",
                                  "Beds:", BEDS, "<br>",
                                  "Owner:", OWNER, "<br>",
                                  "City:", CITY, "<br>",
                                  "County:", COUNTY, "<br>",
                                  "Trauma Level:", TRAUMA, "<br>",
                                  "Helipad:", HELIPAD),
                   weight = 1,
                   radius = ~sqrt(BEDS) * 0.22,  
                   color = ~OWNER,  # Assign colors based on the OWNER variable
                   fillColor = ~colorFactor(palette = c("green", "blue","red", "cyan"), 
                                             domain = c("Government - District/Authority", 
                                                        "Non-profit", 
                                                        "Proprietary","Not Available"))(as.character(OWNER)),  # Convert OWNER to character
                   fillOpacity = 0.8,  # Increase fill opacity
                   opacity = 1,  # Set opacity to 1
                   group = "hospital_markers") %>%  # Add a group name for the main markers
  addLegend("bottomleft",
            pal = colorFactor(palette = c("green", "blue","red", "cyan"), 
                              domain = c("Government - District/Authority", 
                                         "Non-profit", 
                                         "Proprietary",
                                         "Not Available")),
            values = HospitalsAZ$OWNER,
            title = "Hospital Owner",
            opacity = 1)
#Esri.WorldStreetMap") %>% 
#addProviderTiles(provider = "OpenStreetMap") %>% 
#addProviderTiles(provider = "Esri.WorldPhysical") %>% 
#addProviderTiles(provider = "Esri.WorldImagery") %>% 
#addProviderTiles(provider = "Esri.WorldTopoMap") %>% 

Added additional dimensions with the bubble size indicating the hospital size by number of beds and the color indicating the type of owner that the hospital has. We see many more Non-profit and Propriety hospitals than government hospitals. We also see more and larger hospitals in high population areas.The popup also tells the user if the hospital has a helipad and what level of trauma it can treat. Chose legend colors that would be very easy to differentiate from the inferno palette. Chose the openstreetmap provider because it had city names and roads.