Pharmacies
Ishi Agrawal
2023-03-29
Packages needed
#read in an excel file of data
library(readxl)
#To make shiny apps
library(shiny)
#manipulate and clean our datasets
library(tidyverse)## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
## ✔ ggplot2 3.4.0 ✔ purrr 0.3.4
## ✔ tibble 3.1.8 ✔ dplyr 1.0.10
## ✔ tidyr 1.2.1 ✔ stringr 1.4.1
## ✔ readr 2.1.2 ✔ forcats 0.5.2## Warning: package 'ggplot2' was built under R version 4.2.2## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()#for leaflet maps
library(leaflet)## Warning: package 'leaflet' was built under R version 4.2.2#GIS specific packages
library(sf)## Warning: package 'sf' was built under R version 4.2.2## Linking to GEOS 3.9.3, GDAL 3.5.2, PROJ 8.2.1; sf_use_s2() is TRUElibrary(tmap)## Warning: package 'tmap' was built under R version 4.2.2library(tmaptools)## Warning: package 'tmaptools' was built under R version 4.2.2library(terra)## Warning: package 'terra' was built under R version 4.2.2## terra 1.7.3
##
## Attaching package: 'terra'
##
## The following object is masked from 'package:tidyr':
##
## extractlibrary(rgdal)## Loading required package: sp
## Please note that rgdal will be retired by the end of 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
##
## rgdal: version: 1.5-32, (SVN revision 1176)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.4.3, released 2022/04/22
## Path to GDAL shared files: C:/Users/Ishi Agrawal/AppData/Local/R/win-library/4.2/rgdal/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ runtime: Rel. 7.2.1, January 1st, 2021, [PJ_VERSION: 721]
## Path to PROJ shared files: C:/Users/Ishi Agrawal/AppData/Local/R/win-library/4.2/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.5-0
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.
##
## Attaching package: 'rgdal'
##
## The following object is masked from 'package:terra':
##
## projectlibrary(raster)## Warning: package 'raster' was built under R version 4.2.2##
## Attaching package: 'raster'
##
## The following object is masked from 'package:dplyr':
##
## selectReading in and cleaning Washington County data
waCounties <- st_read("county10.shp")## Reading layer `county10' from data source
## `C:\Users\Ishi Agrawal\OneDrive\QGIS\pharmaciesGIS_Project\county10.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 39 features and 39 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 579552.7 ymin: 81835.08 xmax: 2551107 ymax: 1355593
## Projected CRS: NAD83(HARN) / Washington South (ftUS)waCounties <- waCounties %>% dplyr:: rename(county = NAME10)Reading in and cleaning Pharmacy data
I retrieved the pharmacy data from Department of Health
pharmNew = st_read("Pharmacies.gdb")## Reading layer `Pharmacies' from data source
## `C:\Users\Ishi Agrawal\OneDrive\QGIS\pharmaciesGIS_Project\Pharmacies.gdb'
## using driver `OpenFileGDB'
## Simple feature collection with 1876 features and 47 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -13847450 ymin: 5712994 xmax: -13023350 ymax: 6273697
## Projected CRS: WGS 84 / Pseudo-MercatorI cleaned up the dataset by selecting the relavent columns (those pertaining to the location of the pharmacy), and then filtered to make sure I did not have rows where the pharmacy is NOT in washington
I’m going to make sure the CRS of pharmacies matches that of wa counties
pharmacies <- st_transform(pharmNew, crs = st_crs(waCounties))
pharmaciesClean <- pharmacies %>% dplyr::select(inFacility, inAddress1, inCity, inState, inCounty, inZip, acAddress1, Av_score, Geolevel, X_coords, Y_coords, Shape) %>% dplyr::filter(inState == "WA") %>% dplyr::mutate(inState = NULL) %>% dplyr::rename(lat = X_coords, long = Y_coords) %>% stats::na.omit()
head(pharmaciesClean)## Simple feature collection with 6 features and 10 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 805383.9 ymin: 617920 xmax: 2401208 ymax: 924872.5
## Projected CRS: NAD83(HARN) / Washington South (ftUS)
## inFacility inAddress1
## 1 Accredo Health Group Inc 22623 68th Ave S
## 2 Credena Health Pharmacy Sacred Heart 101 W 8th Ave Ste 207653
## 3 Fred Meyer Pharmacy 457 21045 Bothell Everett Hwy
## 4 Grays Harbor Community Hospital Pharmacy West 915 Anderson Dr
## 5 Harrison Medical Center Pharmacy Silverdale Campus 1780 NW Myhre Rd Ste 2160
## 6 Invara 507 Stevens Ave Unit AAA
## inCity inCounty inZip acAddress1 Av_score
## 1 KENT King 98032 22623 68TH AVE S 100
## 2 SPOKANE Spokane 99204 101 W 8TH AVE STE 207653 100
## 3 BOTHELL Snohomish 98021 21045 BOTHELL EVERETT HWY 100
## 4 ABERDEEN Grays Harbor 98520 915 ANDERSON DR 100
## 5 SILVERDALE Kitsap 98383 1780 NW MYHRE RD STE 2160 100
## 6 SULTAN Snohomish 98294 507 W STEVENS AVE UNIT AAA 91
## Geolevel lat long Shape
## 1 CENSUS TRACT -122.24945 47.39916 POINT (1207176 758175.6)
## 2 CENSUS TRACT -117.41303 47.64833 POINT (2401208 859153.4)
## 3 CENSUS TRACT -122.20751 47.80723 POINT (1220783 906776.6)
## 4 CENSUS TRACT -123.84639 46.97922 POINT (805383.9 617920)
## 5 CENSUS TRACT -122.67686 47.65521 POINT (1103931 854177.6)
## 6 CENSUS TRACT -121.81523 47.86190 POINT (1317512 924872.5)Reading in and cleaning median household income data
Since I want to compare pharmacy locations to the income of the surrounding areas, I will read in median household income data.
Let’s remove NA values.
Since the county boundaries are for 2010, I will use relatively recent data, so 2010 to 2019 data.
To clean the data, I only select the years I need, and rename them to include “income” in their column names for easier future reference
median_household_income_estimates <- read_excel("median_household_income_estimates.xlsx",
skip = 3)## New names:
## • `` -> `...1`incomes <- median_household_income_estimates %>%
stats:: na.omit() %>%
dplyr::select(...1, "2010", "2011", "2012", "2013", "2014", "2015", "2016", "2017", "2018", "2019") %>%
rename(
county = ...1,
income2010 = "2010",
income2011 = "2011",
income2012 = "2012",
income2013 = "2013",
income2014 = "2014",
income2015 = "2015",
income2016 = "2016",
income2017 = "2017",
income2018 = "2018",
income2019 = "2019",
) %>%
dplyr::filter(county != "Washington")
head(incomes)## # A tibble: 6 × 11
## county incom…¹ incom…² incom…³ incom…⁴ incom…⁵ incom…⁶ incom…⁷ incom…⁸ incom…⁹
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Adams 40656. 41068. 42354. 43541. 45712. 47646. 49501. 48849. 52870.
## 2 Asotin 39820. 40171. 41703. 42174. 43368. 46107. 47020. 51767. 50746.
## 3 Benton 60070. 60608. 62739. 63710. 63157. 62071. 62282. 63502. 67912.
## 4 Chelan 45478. 46275. 47265. 51713. 50825. 53068. 55109. 60791. 60747.
## 5 Clall… 38397. 38886. 41887. 44824. 45454. 46241. 48187. 47767. 55664.
## 6 Clark 54581. 54951. 56054. 57852. 61711. 63639. 66782. 71922. 71659.
## # … with 1 more variable: income2019 <dbl>, and abbreviated variable names
## # ¹income2010, ²income2011, ³income2012, ⁴income2013, ⁵income2014,
## # ⁶income2015, ⁷income2016, ⁸income2017, ⁹income2018Number of pharmacies in WA
Now let’s see the exact number of pharmacies in WA state by county
We can join the pharmacy and county data now since they have the same CRS now. Now the data has a “point” geometry which represents each pharmacy location, which can be overlayed onto a map of the county boundaries
pharmAndWaCounties <- st_join(pharmaciesClean, waCounties)
head(pharmAndWaCounties)## Simple feature collection with 6 features and 49 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 805383.9 ymin: 617920 xmax: 2401208 ymax: 924872.5
## Projected CRS: NAD83(HARN) / Washington South (ftUS)
## inFacility inAddress1
## 1 Accredo Health Group Inc 22623 68th Ave S
## 2 Credena Health Pharmacy Sacred Heart 101 W 8th Ave Ste 207653
## 3 Fred Meyer Pharmacy 457 21045 Bothell Everett Hwy
## 4 Grays Harbor Community Hospital Pharmacy West 915 Anderson Dr
## 5 Harrison Medical Center Pharmacy Silverdale Campus 1780 NW Myhre Rd Ste 2160
## 6 Invara 507 Stevens Ave Unit AAA
## inCity inCounty inZip acAddress1 Av_score
## 1 KENT King 98032 22623 68TH AVE S 100
## 2 SPOKANE Spokane 99204 101 W 8TH AVE STE 207653 100
## 3 BOTHELL Snohomish 98021 21045 BOTHELL EVERETT HWY 100
## 4 ABERDEEN Grays Harbor 98520 915 ANDERSON DR 100
## 5 SILVERDALE Kitsap 98383 1780 NW MYHRE RD STE 2160 100
## 6 SULTAN Snohomish 98294 507 W STEVENS AVE UNIT AAA 91
## Geolevel lat long STATEFP10 COUNTYFP10 COUNTYNS10 GEOID10
## 1 CENSUS TRACT -122.24945 47.39916 53 033 01531933 53033
## 2 CENSUS TRACT -117.41303 47.64833 53 063 01529225 53063
## 3 CENSUS TRACT -122.20751 47.80723 53 061 01529222 53061
## 4 CENSUS TRACT -123.84639 46.97922 53 027 01531823 53027
## 5 CENSUS TRACT -122.67686 47.65521 53 035 01529223 53035
## 6 CENSUS TRACT -121.81523 47.86190 53 061 01529222 53061
## county NAMELSAD10 LSAD10 CLASSFP10 MTFCC10 CSAFP10 CBSAFP10
## 1 King King County 06 H1 G4020 500 42660
## 2 Spokane Spokane County 06 H1 G4020 <NA> 44060
## 3 Snohomish Snohomish County 06 H1 G4020 500 42660
## 4 Grays Harbor Grays Harbor County 06 H1 G4020 <NA> 10140
## 5 Kitsap Kitsap County 06 H1 G4020 500 14740
## 6 Snohomish Snohomish County 06 H1 G4020 500 42660
## METDIVFP10 FUNCSTAT10 INTPTLON10 INTPTLAT10 ALANDM AWATERM ALANDMI
## 1 42644 A -121.8324 47.49355 5479291809 495471264 2115.566
## 2 <NA> A -117.4044 47.62038 4568198649 43802705 1763.791
## 3 42644 A -121.7664 48.05491 5406013264 282371494 2087.273
## 4 <NA> A -123.8270 47.14279 4926225833 834957671 1902.027
## 5 <NA> A -122.6496 47.63969 1022893152 442829766 394.941
## 6 42644 A -121.7664 48.05491 5406013264 282371494 2087.273
## AWATERMI POP10 HHP10 GQ10 HU10 OHU10 POPWHITE POPBLACK POPAIAN
## 1 191.303 1931249 1894118 37131 851261 789232 1325845 119801 16147
## 2 16.912 471221 456529 14692 201434 187167 420275 8056 7295
## 3 109.024 713335 702938 10397 286659 268325 559011 18168 9793
## 4 322.379 72797 70066 2731 35166 28579 61825 803 3325
## 5 170.978 251133 242411 8722 107367 97220 207369 6650 4040
## 6 109.024 713335 702938 10397 286659 268325 559011 18168 9793
## POPASIAN POPNHOPI POPOTH POPTWO POPHISP POPWHITE2 POPBLACK2 POPAIAN2
## 1 282075 14486 76096 96799 172378 1251300 116326 12931
## 2 9957 1902 5880 17856 21260 408629 7714 6478
## 3 63385 3135 27121 32722 64249 529814 17299 8451
## 4 1032 182 2817 2813 6272 59282 762 3005
## 5 12396 2310 3919 14449 15686 198745 6329 3524
## 6 63385 3135 27121 32722 64249 529814 17299 8451
## POPASIAN2 POPNHOPI2 POPOTH2 POPTWO2 Shape
## 1 280029 14068 4688 79529 POINT (1207176 758175.6)
## 2 9799 1817 610 14914 POINT (2401208 859153.4)
## 3 62760 3002 1220 26540 POINT (1220783 906776.6)
## 4 995 177 72 2232 POINT (805383.9 617920)
## 5 12082 2177 423 12167 POINT (1103931 854177.6)
## 6 62760 3002 1220 26540 POINT (1317512 924872.5)Then we use geom_sf on both data sets and theme_void to remove the lattitute and longitude lines
Running this we can see clusters of pharmacies in some areas, probably those of higher populations, whereas other counties have more dispersed pharmacy locations.
ggplot() +
geom_sf(data = waCounties) +
geom_sf(data = pharmAndWaCounties) +
theme_void()
## Looking at the pharmacy locations in wa counties with respect to
income
How does this distribution compare against income levels?
First, let’s quantify income levels into groups for a specific year We will use “cut” to break up each level of income data to represent “income brackets”
incomeLevel <- quantile(incomes$income2010, probs = seq(0, 1, 1/3))
incomeLevel## 0% 33.33333% 66.66667% 100%
## 31061.78 40845.31 47946.63 65383.49incomesWithLevels <- incomes %>% dplyr::mutate(
incomeBracket10 = cut(incomes$income2010, breaks = quantile(incomes$income2010, breaks = quantile(incomes$income2010, probs = seq(0, 1, 1/3), include.lowest = T)), dig.lab = 7),
incomeBracket11 = cut(incomes$income2011, breaks = quantile(incomes$income2011, breaks = quantile(incomes$income2011, probs = seq(0, 1, 1/3), include.lowest = T)), dig.lab = 7),
incomeBracket12 = cut(incomes$income2012, breaks = quantile(incomes$income2012, breaks = quantile(incomes$income2012, probs = seq(0, 1, 1/3), include.lowest = T)), dig.lab = 7),
incomeBracket13 = cut(incomes$income2013, breaks = quantile(incomes$income2010, breaks = quantile(incomes$income2013, probs = seq(0, 1, 1/3), include.lowest = T)), dig.lab = 7),
incomeBracket14 = cut(incomes$income2014, breaks = quantile(incomes$income2010, breaks = quantile(incomes$income2014, probs = seq(0, 1, 1/3), include.lowest = T)), dig.lab = 7),
incomeBracket15 = cut(incomes$income2015, breaks = quantile(incomes$income2010, breaks = quantile(incomes$income2015, probs = seq(0, 1, 1/3), include.lowest = T)), dig.lab = 7),
incomeBracket16 = cut(incomes$income2016, breaks = quantile(incomes$income2010, breaks = quantile(incomes$income2016, probs = seq(0, 1, 1/3), include.lowest = T)), dig.lab = 7),
incomeBracket17 = cut(incomes$income2017, breaks = quantile(incomes$income2010, breaks = quantile(incomes$income2017, probs = seq(0, 1, 1/3), include.lowest = T)), dig.lab = 7),
incomeBracket18 = cut(incomes$income2018, breaks = quantile(incomes$income2010, breaks = quantile(incomes$income2018, probs = seq(0, 1, 1/3), include.lowest = T)), dig.lab = 7),
incomeBracket19 = cut(incomes$income2019, breaks = quantile(incomes$income2010, breaks = quantile(incomes$income2019, probs = seq(0, 1, 1/3), include.lowest = T)), dig.lab = 7)
)
head(incomesWithLevels)## # A tibble: 6 × 21
## county incom…¹ incom…² incom…³ incom…⁴ incom…⁵ incom…⁶ incom…⁷ incom…⁸ incom…⁹
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Adams 40656. 41068. 42354. 43541. 45712. 47646. 49501. 48849. 52870.
## 2 Asotin 39820. 40171. 41703. 42174. 43368. 46107. 47020. 51767. 50746.
## 3 Benton 60070. 60608. 62739. 63710. 63157. 62071. 62282. 63502. 67912.
## 4 Chelan 45478. 46275. 47265. 51713. 50825. 53068. 55109. 60791. 60747.
## 5 Clall… 38397. 38886. 41887. 44824. 45454. 46241. 48187. 47767. 55664.
## 6 Clark 54581. 54951. 56054. 57852. 61711. 63639. 66782. 71922. 71659.
## # … with 11 more variables: income2019 <dbl>, incomeBracket10 <fct>,
## # incomeBracket11 <fct>, incomeBracket12 <fct>, incomeBracket13 <fct>,
## # incomeBracket14 <fct>, incomeBracket15 <fct>, incomeBracket16 <fct>,
## # incomeBracket17 <fct>, incomeBracket18 <fct>, incomeBracket19 <fct>, and
## # abbreviated variable names ¹income2010, ²income2011, ³income2012,
## # ⁴income2013, ⁵income2014, ⁶income2015, ⁷income2016, ⁸income2017,
## # ⁹income2018We can also extract specifically 2010 information
incomesWithLevels10 <-incomesWithLevels %>% dplyr::select(county, income2010, incomeBracket10)
incomesWithLevels10## # A tibble: 39 × 3
## county income2010 incomeBracket10
## <chr> <dbl> <fct>
## 1 Adams 40656. (39913.64,43915.28]
## 2 Asotin 39820. (31061.78,39913.64]
## 3 Benton 60070. (53198.07,65383.49]
## 4 Chelan 45478. (43915.28,53198.07]
## 5 Clallam 38397. (31061.78,39913.64]
## 6 Clark 54581. (53198.07,65383.49]
## 7 Columbia 38474. (31061.78,39913.64]
## 8 Cowlitz 40867. (39913.64,43915.28]
## 9 Douglas 46159. (43915.28,53198.07]
## 10 Ferry 36712. (31061.78,39913.64]
## # … with 29 more rowsLet’s start with a general map of pharmacies in Washington state
First we join wa county data with income data so that each row of our data has county geometry AND the income brackets associated with each county for EVERY year
incomeAndWaCounties <- waCounties %>%
dplyr::left_join(incomesWithLevels, by= "county")
head(incomeAndWaCounties)## Simple feature collection with 6 features and 59 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 579552.7 ymin: 81835.08 xmax: 2551107 ymax: 1182898
## Projected CRS: NAD83(HARN) / Washington South (ftUS)
## STATEFP10 COUNTYFP10 COUNTYNS10 GEOID10 county NAMELSAD10 LSAD10
## 1 53 001 01531601 53001 Adams Adams County 06
## 2 53 003 01533502 53003 Asotin Asotin County 06
## 3 53 005 01513302 53005 Benton Benton County 06
## 4 53 007 01531932 53007 Chelan Chelan County 06
## 5 53 009 01531341 53009 Clallam Clallam County 06
## 6 53 011 01531820 53011 Clark Clark County 06
## CLASSFP10 MTFCC10 CSAFP10 CBSAFP10 METDIVFP10 FUNCSTAT10 INTPTLON10
## 1 H1 G4020 <NA> <NA> <NA> A -118.5333
## 2 H1 G4020 <NA> 30300 <NA> A -117.2278
## 3 H1 G4020 <NA> 28420 <NA> A -119.5169
## 4 H1 G4020 <NA> 48300 <NA> A -120.6185
## 5 H1 G4020 <NA> 38820 <NA> A -123.9306
## 6 H1 G4020 <NA> 38900 <NA> A -122.4859
## INTPTLAT10 ALANDM AWATERM ALANDMI AWATERMI POP10 HHP10 GQ10 HU10
## 1 47.00484 4985668796 12623299 1924.978 4.874 18728 18566 162 6242
## 2 46.18186 1647784012 11475136 636.213 4.431 21623 21449 174 9872
## 3 46.22807 4403961281 154412831 1700.379 59.619 175177 173751 1426 68618
## 4 47.85989 7564124657 189992686 2920.525 73.357 72453 71525 928 35465
## 5 48.11301 4502251295 2414668004 1738.329 932.309 71404 69505 1899 35582
## 6 45.77167 1629112686 70523141 629.004 27.229 425363 422153 3210 167413
## OHU10 POPWHITE POPBLACK POPAIAN POPASIAN POPNHOPI POPOTH POPTWO POPHISP
## 1 5720 11703 109 356 125 4 5899 532 11099
## 2 9236 20392 92 302 117 36 163 521 643
## 3 65304 144418 2221 1574 4691 253 15798 6222 32696
## 4 27827 57484 236 700 588 100 11355 1990 18713
## 5 31329 62092 596 3630 1007 94 1269 2716 3627
## 6 158099 363397 8426 3624 17504 2708 12485 17219 32166
## POPWHITE2 POPBLACK2 POPAIAN2 POPASIAN2 POPNHOPI2 POPOTH2 POPTWO2 income2010
## 1 7262 47 76 99 3 18 124 40655.90
## 2 20026 89 272 114 36 12 431 39819.69
## 3 130437 2031 1280 4621 221 260 3631 60070.36
## 4 51202 177 514 570 89 76 1112 45477.70
## 5 60400 558 3326 993 86 114 2300 38396.84
## 6 347793 8025 2970 17276 2589 639 13905 54580.91
## income2011 income2012 income2013 income2014 income2015 income2016 income2017
## 1 41068.44 42353.54 43541.08 45711.82 47645.73 49501.39 48848.79
## 2 40171.03 41703.24 42173.78 43367.86 46107.27 47020.11 51766.88
## 3 60608.28 62739.22 63709.96 63157.40 62071.02 62282.05 63501.71
## 4 46274.71 47265.40 51712.87 50825.03 53067.96 55109.14 60790.91
## 5 38885.97 41886.68 44823.81 45453.82 46240.92 48187.11 47766.59
## 6 54950.94 56054.15 57852.10 61710.84 63639.18 66782.17 71922.31
## income2018 income2019 incomeBracket10 incomeBracket11
## 1 52870.06 53535 (39913.64,43915.28] (40226.43,44606.01]
## 2 50746.37 54776 (31061.78,39913.64] (31396.19,40226.43]
## 3 67912.30 72847 (53198.07,65383.49] (53779.94,66293.71]
## 4 60747.40 59838 (43915.28,53198.07] (44606.01,53779.94]
## 5 55664.09 57571 (31061.78,39913.64] (31396.19,40226.43]
## 6 71658.81 80407 (53198.07,65383.49] (53779.94,66293.71]
## incomeBracket12 incomeBracket13 incomeBracket14
## 1 (41794.96,45689.63] (39913.64,43915.28] (43915.28,53198.07]
## 2 (32569.72,41794.96] (39913.64,43915.28] (39913.64,43915.28]
## 3 (55057.98,68312.79] (53198.07,65383.49] (53198.07,65383.49]
## 4 (45689.63,55057.98] (43915.28,53198.07] (43915.28,53198.07]
## 5 (41794.96,45689.63] (43915.28,53198.07] (43915.28,53198.07]
## 6 (55057.98,68312.79] (53198.07,65383.49] (53198.07,65383.49]
## incomeBracket15 incomeBracket16 incomeBracket17
## 1 (43915.28,53198.07] (43915.28,53198.07] (43915.28,53198.07]
## 2 (43915.28,53198.07] (43915.28,53198.07] (43915.28,53198.07]
## 3 (53198.07,65383.49] (53198.07,65383.49] (53198.07,65383.49]
## 4 (43915.28,53198.07] (53198.07,65383.49] (53198.07,65383.49]
## 5 (43915.28,53198.07] (43915.28,53198.07] (43915.28,53198.07]
## 6 (53198.07,65383.49] <NA> <NA>
## incomeBracket18 incomeBracket19 geometry
## 1 (43915.28,53198.07] (53198.07,65383.49] POLYGON ((2019638 579100.4,...
## 2 (43915.28,53198.07] (53198.07,65383.49] POLYGON ((2419090 359524.8,...
## 3 <NA> <NA> POLYGON ((1801305 290122.6,...
## 4 (53198.07,65383.49] (53198.07,65383.49] POLYGON ((1530000 780445.8,...
## 5 (53198.07,65383.49] (53198.07,65383.49] POLYGON ((904078.9 1070275,...
## 6 <NA> <NA> POLYGON ((1061154 171892.5,...Now we can use this new data set to plot for 2010 The first geom_sf creates a chlorpleth for WA counties based on 2010 income brackets using the “POLYGON” geometry representing a whole county When we overlay the second geom_sf, then it plots the “POINT” geometry of each individual pharmacy location
PharmaciesByIncomeCountyMap <- ggplot() +
geom_sf(data = incomeAndWaCounties, aes(fill = incomeBracket10)) +
scale_fill_brewer(name = expression("Income levels"),
palette = "Greens") +
theme_void() +
geom_sf(data = pharmAndWaCounties)
PharmaciesByIncomeCountyMap
What if we just want to look at a map of wa counties and see the number of pharmacies, not the individual pharmacies for the county?
We can try using leaflet
To set up the leaflet data, we need to group the pharmacy and county data by county and sum up the number of pharmacies in each county. Then we drop the geometry since we do not need the specific point geometry of the pharmacy, but rather just the polygon geometry of the county.
This dataset is then joined with waCounties so that we can add a column with the number of pharmacies for each county while having the county polygon data so that boundaries can be drawn later.
pharmaciesGroupedByCounty <- pharmAndWaCounties %>% dplyr::select(inCounty) %>% dplyr::group_by(inCounty) %>% dplyr::summarise(numOfPharm = n()) %>% st_drop_geometry()
head(pharmaciesGroupedByCounty)## # A tibble: 6 × 2
## inCounty numOfPharm
## <chr> <int>
## 1 " " 2
## 2 "Adams" 7
## 3 "Asotin" 19
## 4 "Benton" 50
## 5 "Chelan" 35
## 6 "Clallam" 31head(waCounties)## Simple feature collection with 6 features and 39 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 579552.7 ymin: 81835.08 xmax: 2551107 ymax: 1182898
## Projected CRS: NAD83(HARN) / Washington South (ftUS)
## STATEFP10 COUNTYFP10 COUNTYNS10 GEOID10 county NAMELSAD10 LSAD10
## 1 53 001 01531601 53001 Adams Adams County 06
## 2 53 003 01533502 53003 Asotin Asotin County 06
## 3 53 005 01513302 53005 Benton Benton County 06
## 4 53 007 01531932 53007 Chelan Chelan County 06
## 5 53 009 01531341 53009 Clallam Clallam County 06
## 6 53 011 01531820 53011 Clark Clark County 06
## CLASSFP10 MTFCC10 CSAFP10 CBSAFP10 METDIVFP10 FUNCSTAT10 INTPTLON10
## 1 H1 G4020 <NA> <NA> <NA> A -118.5333
## 2 H1 G4020 <NA> 30300 <NA> A -117.2278
## 3 H1 G4020 <NA> 28420 <NA> A -119.5169
## 4 H1 G4020 <NA> 48300 <NA> A -120.6185
## 5 H1 G4020 <NA> 38820 <NA> A -123.9306
## 6 H1 G4020 <NA> 38900 <NA> A -122.4859
## INTPTLAT10 ALANDM AWATERM ALANDMI AWATERMI POP10 HHP10 GQ10 HU10
## 1 47.00484 4985668796 12623299 1924.978 4.874 18728 18566 162 6242
## 2 46.18186 1647784012 11475136 636.213 4.431 21623 21449 174 9872
## 3 46.22807 4403961281 154412831 1700.379 59.619 175177 173751 1426 68618
## 4 47.85989 7564124657 189992686 2920.525 73.357 72453 71525 928 35465
## 5 48.11301 4502251295 2414668004 1738.329 932.309 71404 69505 1899 35582
## 6 45.77167 1629112686 70523141 629.004 27.229 425363 422153 3210 167413
## OHU10 POPWHITE POPBLACK POPAIAN POPASIAN POPNHOPI POPOTH POPTWO POPHISP
## 1 5720 11703 109 356 125 4 5899 532 11099
## 2 9236 20392 92 302 117 36 163 521 643
## 3 65304 144418 2221 1574 4691 253 15798 6222 32696
## 4 27827 57484 236 700 588 100 11355 1990 18713
## 5 31329 62092 596 3630 1007 94 1269 2716 3627
## 6 158099 363397 8426 3624 17504 2708 12485 17219 32166
## POPWHITE2 POPBLACK2 POPAIAN2 POPASIAN2 POPNHOPI2 POPOTH2 POPTWO2
## 1 7262 47 76 99 3 18 124
## 2 20026 89 272 114 36 12 431
## 3 130437 2031 1280 4621 221 260 3631
## 4 51202 177 514 570 89 76 1112
## 5 60400 558 3326 993 86 114 2300
## 6 347793 8025 2970 17276 2589 639 13905
## geometry
## 1 POLYGON ((2019638 579100.4,...
## 2 POLYGON ((2419090 359524.8,...
## 3 POLYGON ((1801305 290122.6,...
## 4 POLYGON ((1530000 780445.8,...
## 5 POLYGON ((904078.9 1070275,...
## 6 POLYGON ((1061154 171892.5,...pharmaciesGroupedByCountyWithGeo <- waCounties %>%
dplyr::left_join(pharmaciesGroupedByCounty, by= c("county" = "inCounty"))
head(pharmaciesGroupedByCountyWithGeo)## Simple feature collection with 6 features and 40 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 579552.7 ymin: 81835.08 xmax: 2551107 ymax: 1182898
## Projected CRS: NAD83(HARN) / Washington South (ftUS)
## STATEFP10 COUNTYFP10 COUNTYNS10 GEOID10 county NAMELSAD10 LSAD10
## 1 53 001 01531601 53001 Adams Adams County 06
## 2 53 003 01533502 53003 Asotin Asotin County 06
## 3 53 005 01513302 53005 Benton Benton County 06
## 4 53 007 01531932 53007 Chelan Chelan County 06
## 5 53 009 01531341 53009 Clallam Clallam County 06
## 6 53 011 01531820 53011 Clark Clark County 06
## CLASSFP10 MTFCC10 CSAFP10 CBSAFP10 METDIVFP10 FUNCSTAT10 INTPTLON10
## 1 H1 G4020 <NA> <NA> <NA> A -118.5333
## 2 H1 G4020 <NA> 30300 <NA> A -117.2278
## 3 H1 G4020 <NA> 28420 <NA> A -119.5169
## 4 H1 G4020 <NA> 48300 <NA> A -120.6185
## 5 H1 G4020 <NA> 38820 <NA> A -123.9306
## 6 H1 G4020 <NA> 38900 <NA> A -122.4859
## INTPTLAT10 ALANDM AWATERM ALANDMI AWATERMI POP10 HHP10 GQ10 HU10
## 1 47.00484 4985668796 12623299 1924.978 4.874 18728 18566 162 6242
## 2 46.18186 1647784012 11475136 636.213 4.431 21623 21449 174 9872
## 3 46.22807 4403961281 154412831 1700.379 59.619 175177 173751 1426 68618
## 4 47.85989 7564124657 189992686 2920.525 73.357 72453 71525 928 35465
## 5 48.11301 4502251295 2414668004 1738.329 932.309 71404 69505 1899 35582
## 6 45.77167 1629112686 70523141 629.004 27.229 425363 422153 3210 167413
## OHU10 POPWHITE POPBLACK POPAIAN POPASIAN POPNHOPI POPOTH POPTWO POPHISP
## 1 5720 11703 109 356 125 4 5899 532 11099
## 2 9236 20392 92 302 117 36 163 521 643
## 3 65304 144418 2221 1574 4691 253 15798 6222 32696
## 4 27827 57484 236 700 588 100 11355 1990 18713
## 5 31329 62092 596 3630 1007 94 1269 2716 3627
## 6 158099 363397 8426 3624 17504 2708 12485 17219 32166
## POPWHITE2 POPBLACK2 POPAIAN2 POPASIAN2 POPNHOPI2 POPOTH2 POPTWO2 numOfPharm
## 1 7262 47 76 99 3 18 124 7
## 2 20026 89 272 114 36 12 431 19
## 3 130437 2031 1280 4621 221 260 3631 50
## 4 51202 177 514 570 89 76 1112 35
## 5 60400 558 3326 993 86 114 2300 31
## 6 347793 8025 2970 17276 2589 639 13905 82
## geometry
## 1 POLYGON ((2019638 579100.4,...
## 2 POLYGON ((2419090 359524.8,...
## 3 POLYGON ((1801305 290122.6,...
## 4 POLYGON ((1530000 780445.8,...
## 5 POLYGON ((904078.9 1070275,...
## 6 POLYGON ((1061154 171892.5,...PharmaciesPerCountyLeafletMap <- leaflet() %>% addTiles() %>%
addCircleMarkers(
data = pharmaciesGroupedByCountyWithGeo,
lng = ~INTPTLON10,
lat = ~INTPTLAT10,
radius = pharmaciesGroupedByCountyWithGeo$numOfPharm / 10,
label = paste(pharmaciesGroupedByCountyWithGeo$county, "county:", pharmaciesGroupedByCountyWithGeo$numOfPharm)
)
PharmaciesPerCountyLeafletMapNow I want to see this data along with the income bracketthat the county falls under
To do this, we need the dataset with the number of pharmacies per county to be joined with the dataset with the income brackets each county falls under.
pharmaciesGroupedByCountyIncomesWithGeo <- pharmaciesGroupedByCountyWithGeo %>%
dplyr::left_join(incomesWithLevels10, by= c("county" = "county"))
levels(pharmaciesGroupedByCountyIncomesWithGeo$incomeBracket10) <- c('very low', 'medium', 'high', 'highest')incomeColors <- colorFactor("inferno", pharmaciesGroupedByCountyIncomesWithGeo$incomeBracket10)
PharmaciesPerCountyColoredIncomeLeafletMap <- leaflet() %>% addTiles() %>%
addCircleMarkers(
data = pharmaciesGroupedByCountyIncomesWithGeo,
lng = ~INTPTLON10,
lat = ~INTPTLAT10,
radius = pharmaciesGroupedByCountyIncomesWithGeo$numOfPharm/4,
label = paste(pharmaciesGroupedByCountyIncomesWithGeo$county, "county:", pharmaciesGroupedByCountyIncomesWithGeo$numOfPharm, pharmaciesGroupedByCountyIncomesWithGeo$incomeBracket10),
color = ~incomeColors(pharmaciesGroupedByCountyIncomesWithGeo$incomeBracket10),
weight = 3
) %>%
addLegend("bottomright", pal = incomeColors, values = pharmaciesGroupedByCountyIncomesWithGeo$incomeBracket10,
title = "Median Income (2010)",
opacity = 1
)
PharmaciesPerCountyColoredIncomeLeafletMapNEW TOPIC: Adding interactivity with Shiny
Now what if I want to look at the pharmacies by different income levels where I select a specific income level
One option is to manually make a new separate map for a specific income, like “very low”
A better option is to make it interactive where we dynamically choose the income level desired
This can be done with SHINY
# Define UI ----
ui <- fluidPage(
titlePanel("Pharmacies By Income Level"),
selectInput("IncomeLevelSelect", label = h3("View Number of Pharmacies by Income Levels"),
choices = levels(pharmaciesGroupedByCountyIncomesWithGeo$incomeBracket10)),
leafletOutput("PharmaciesPerCountyColoredIncomeLeafletMap")
)
# Define server logic ----
server <- function(input, output) {
filtered <- reactive({
pharmaciesGroupedByCountyIncomesWithGeo %>% dplyr::filter(pharmaciesGroupedByCountyIncomesWithGeo$incomeBracket10 == input$IncomeLevelSelect) %>%
st_drop_geometry() %>%
st_as_sf(coords = c('INTPTLON10', 'INTPTLAT10'))
})
output$PharmaciesPerCountyColoredIncomeLeafletMap <- renderLeaflet({
leaflet() %>% addTiles() %>%
addCircleMarkers(
data = filtered(),
radius = filtered()$numOfPharm,
label = paste(filtered()$NAMELSAD10, "Pharmacies: ", filtered()$numOfPharm, "Pharmacies, Income: ", round(filtered()$income2010, digits = 0)),
weight = 3,
color = ~incomeColors(filtered()$incomeBracket10),
) %>%
addLegend("bottomright", pal = incomeColors, values = filtered()$incomeBracket10,
title = "Median Income (2010)",
opacity = 1
)
})
}
# Run the app ----
shinyApp(ui = ui, server = server)Shiny applications not supported in static R Markdown documents
What if we want to just see which counties have some minimum number of pharmacies
# Define UI ----
ui <- fluidPage(
titlePanel("Find Counties with some number of pharmacies"),
sliderInput("slider", label = h3("Counties with a minimum number of pharmacies"), min = min(pharmaciesGroupedByCountyIncomesWithGeo$numOfPharm),
max = max(pharmaciesGroupedByCountyIncomesWithGeo$numOfPharm), value = 50),
leafletOutput("PharmaciesLeafletMap")
)
# Define server logic ----
server <- function(input, output) {
filtered2 <- reactive({
pharmaciesGroupedByCountyIncomesWithGeo %>% dplyr::filter(pharmaciesGroupedByCountyIncomesWithGeo$numOfPharm >= input$slider) %>%
st_drop_geometry() %>%
st_as_sf(coords = c('INTPTLON10', 'INTPTLAT10'))
})
output$PharmaciesLeafletMap <- renderLeaflet({
leaflet() %>% addTiles() %>%
addCircleMarkers(
data = filtered2(),
radius = filtered2()$numOfPharm/5,
label = paste(filtered2()$NAMELSAD10, "Pharmacies: ", filtered2()$numOfPharm, "Pharmacies, Income: ", round(filtered2()$income2010, digits = 0)),
weight = 3,
color = ~incomeColors(filtered2()$incomeBracket10),
) %>%
addLegend("bottomright", pal = incomeColors, values = filtered2()$incomeBracket10,
title = "Median Income (2010)",
opacity = 1
)
})
}
# Run the app ----
shinyApp(ui = ui, server = server)Shiny applications not supported in static R Markdown documents