Objectives

Using the spatial visualization techniques, explore this data set on Pennsylvania hospitals (http://www.arcgis.com/home/item.html?id=eccee5dfe01e4c4283c9be0cfc596882). Create a series of 5 maps that highlight spatial differences in hospital service coverage for the state of PA.

To help you in getting the data imported into R, I have included the code below:

To import the data I use the foreign package, if you do not have it than be sure to install it prior to testing the code.

library(foreign)
library(leaflet)
## Warning: package 'leaflet' was built under R version 3.5.3
library(data.table)
library(dplyr)
## Warning: package 'dplyr' was built under R version 3.5.2
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:data.table':
## 
##     between, first, last
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(ggmap)
## Warning: package 'ggmap' was built under R version 3.5.3
## Loading required package: ggplot2
## Warning: package 'ggplot2' was built under R version 3.5.3
## Google's Terms of Service: https://cloud.google.com/maps-platform/terms/.
## Please cite ggmap if you use it! See citation("ggmap") for details.
library(plyr)
## Warning: package 'plyr' was built under R version 3.5.3
## -------------------------------------------------------------------------
## You have loaded plyr after dplyr - this is likely to cause problems.
## If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
## library(plyr); library(dplyr)
## -------------------------------------------------------------------------
## 
## Attaching package: 'plyr'
## The following objects are masked from 'package:dplyr':
## 
##     arrange, count, desc, failwith, id, mutate, rename, summarise,
##     summarize
dat <- read.dbf("C:/Users/Priya/Downloads/pennsylv/pennsylv.dbf")

The dataset contains a number of variables about each hospital, many of them are clear and straight forward.

names(dat)
##   [1] "acc_trauma" "air_amb"    "als"        "arc_street" "arc_zone"  
##   [6] "bas_ls"     "bassinets"  "bb_id"      "bc_beds"    "bc_sus_bed"
##  [11] "beds_sus"   "birthing_r" "bone_marro" "burn_car"   "burn_care" 
##  [16] "card_beds"  "card_surge" "card_sus_b" "cardiac"    "cardiac_ca"
##  [21] "cardio_reh" "chemo"      "city"       "clin_lab"   "clin_psych"
##  [26] "county"     "countyname" "ct_scan"    "cty_key"    "cystoscopi"
##  [31] "deliv_rms"  "dental"     "detox_alc_" "diag_radio" "diag_xray" 
##  [36] "doh_hosp"   "doh_phone"  "emer_dept"  "endoscopie" "fac_id"    
##  [41] "facility"   "flu_old"    "fred_con_1" "fred_conta" "fred_email"
##  [46] "fred_fax"   "fred_hosp"  "fred_pager" "fred_phone" "gamma_knif"
##  [51] "gen_outpat" "gene_counc" "heart_tran" "helipad"    "hemodial_c"
##  [56] "hemodial_m" "hosp_id"    "hospice"    "hyper_cham" "icu"       
##  [61] "icu_beds"   "icu_sus_be" "inpat_flu_" "inpat_pneu" "kidney_tra"
##  [66] "labor_rms"  "lic_beds"   "lic_dent"   "lic_dos"    "lic_mds"   
##  [71] "lic_pod"    "linear_acc" "lithotrips" "liver_tran" "loc_method"
##  [76] "ltc"        "mcd"        "mcd_key"    "mcd_name"   "medical"   
##  [81] "mob_ccu"    "mob_icu"    "mri"        "ms1"        "neo2_beds" 
##  [86] "neo2_sus_b" "neo3_beds"  "neo3_sus_b" "neuro_surg" "neurology" 
##  [91] "obs_gyn"    "occ_ther"   "optometry"  "organ_bank" "ped_trauma"
##  [96] "pediatric"  "pet"        "pharmacy"   "phys_med"   "phys_ther" 
## [101] "podiatry"   "providerid" "psych"      "psych_inpa" "reg_trauma"
## [106] "resp_ther"  "so_flu_65u" "social_wor" "speech_pat" "street"    
## [111] "surgical"   "surgical_s" "thera_radi" "typ_org"    "typ_serv"  
## [116] "ultrasound" "x"          "y"          "zip"

Now create 5 maps, including descriptions, that highlight the spatial distribution of hosptical services in the state of PA. Upload these maps as a document to rpubs.com and submit that link the Moodle assignment.

qmplot(x,y,data=dat,colour=I('blue'),size=I(3),darken=.3, main = "All Hospitals in Pennsylvania", xlab = "Longitude", ylab = "Latitude")
## Using zoom = 8...
## Source : http://tile.stamen.com/terrain/8/70/94.png
## Source : http://tile.stamen.com/terrain/8/71/94.png
## Source : http://tile.stamen.com/terrain/8/72/94.png
## Source : http://tile.stamen.com/terrain/8/73/94.png
## Source : http://tile.stamen.com/terrain/8/74/94.png
## Source : http://tile.stamen.com/terrain/8/70/95.png
## Source : http://tile.stamen.com/terrain/8/71/95.png
## Source : http://tile.stamen.com/terrain/8/72/95.png
## Source : http://tile.stamen.com/terrain/8/73/95.png
## Source : http://tile.stamen.com/terrain/8/74/95.png
## Source : http://tile.stamen.com/terrain/8/70/96.png
## Source : http://tile.stamen.com/terrain/8/71/96.png
## Source : http://tile.stamen.com/terrain/8/72/96.png
## Source : http://tile.stamen.com/terrain/8/73/96.png
## Source : http://tile.stamen.com/terrain/8/74/96.png
## Source : http://tile.stamen.com/terrain/8/70/97.png
## Source : http://tile.stamen.com/terrain/8/71/97.png
## Source : http://tile.stamen.com/terrain/8/72/97.png
## Source : http://tile.stamen.com/terrain/8/73/97.png
## Source : http://tile.stamen.com/terrain/8/74/97.png

# From Map 1, we can see that the major cities of Pittsburgh & Philadelphia has a higher concentration of hospitals in Pennsylvania.
icu = subset(dat, icu == "Y")
icuPHLMO = subset(icu, county == "Philadelphia" | county =="Montgomery")
qmplot(x,y, data = icuPHLMO, colour = I('red'),size = I(3) , main = "Comparison Between Philadelphia & Montgomery Hospitals with ICU", xlab = "Longitude", ylab = "Latitude")
## Using zoom = 12...
## Source : http://tile.stamen.com/terrain/12/1190/1548.png
## Source : http://tile.stamen.com/terrain/12/1191/1548.png
## Source : http://tile.stamen.com/terrain/12/1192/1548.png
## Source : http://tile.stamen.com/terrain/12/1193/1548.png
## Source : http://tile.stamen.com/terrain/12/1194/1548.png
## Source : http://tile.stamen.com/terrain/12/1190/1549.png
## Source : http://tile.stamen.com/terrain/12/1191/1549.png
## Source : http://tile.stamen.com/terrain/12/1192/1549.png
## Source : http://tile.stamen.com/terrain/12/1193/1549.png
## Source : http://tile.stamen.com/terrain/12/1194/1549.png
## Source : http://tile.stamen.com/terrain/12/1190/1550.png
## Source : http://tile.stamen.com/terrain/12/1191/1550.png
## Source : http://tile.stamen.com/terrain/12/1192/1550.png
## Source : http://tile.stamen.com/terrain/12/1193/1550.png
## Source : http://tile.stamen.com/terrain/12/1194/1550.png
## Source : http://tile.stamen.com/terrain/12/1190/1551.png
## Source : http://tile.stamen.com/terrain/12/1191/1551.png
## Source : http://tile.stamen.com/terrain/12/1192/1551.png
## Source : http://tile.stamen.com/terrain/12/1193/1551.png
## Source : http://tile.stamen.com/terrain/12/1194/1551.png

# Philadelphia county has higher number of ICU hospitals than Montgomery county.
organbank = subset(dat, organ_bank=="Y")
organbankPI = subset(organbank, city == "Pittsburgh")
qmplot(x,y, data = organbankPI, colour = I('blue'), size = I(3), main = "Philadelphia Hospitals with Organ Bank", xlab = "Longitude", ylab = "Latitude")
## Using zoom = 13...
## Source : http://tile.stamen.com/terrain/13/2274/3083.png
## Source : http://tile.stamen.com/terrain/13/2275/3083.png
## Source : http://tile.stamen.com/terrain/13/2276/3083.png
## Source : http://tile.stamen.com/terrain/13/2274/3084.png
## Source : http://tile.stamen.com/terrain/13/2275/3084.png
## Source : http://tile.stamen.com/terrain/13/2276/3084.png
## Source : http://tile.stamen.com/terrain/13/2274/3085.png
## Source : http://tile.stamen.com/terrain/13/2275/3085.png
## Source : http://tile.stamen.com/terrain/13/2276/3085.png
## Source : http://tile.stamen.com/terrain/13/2274/3086.png
## Source : http://tile.stamen.com/terrain/13/2275/3086.png
## Source : http://tile.stamen.com/terrain/13/2276/3086.png
## Source : http://tile.stamen.com/terrain/13/2274/3087.png
## Source : http://tile.stamen.com/terrain/13/2275/3087.png
## Source : http://tile.stamen.com/terrain/13/2276/3087.png
## Source : http://tile.stamen.com/terrain/13/2274/3088.png
## Source : http://tile.stamen.com/terrain/13/2275/3088.png
## Source : http://tile.stamen.com/terrain/13/2276/3088.png

# There are 4 hospitals with an organ bank facility in Pittsburgh. 
organbank = subset(dat, organ_bank=="Y")
organbankPI = subset(organbank, city == "Philadelphia")
qmplot(x,y, data = organbankPI, colour = I('blue'), size = I(3), main = "Philadelphia Hospitals with Organ Bank", xlab = "Longitude", ylab = "Latitude")
## Using zoom = 15...
## Source : http://tile.stamen.com/terrain/15/9537/12408.png
## Source : http://tile.stamen.com/terrain/15/9538/12408.png
## Source : http://tile.stamen.com/terrain/15/9539/12408.png
## Source : http://tile.stamen.com/terrain/15/9540/12408.png
## Source : http://tile.stamen.com/terrain/15/9541/12408.png
## Source : http://tile.stamen.com/terrain/15/9542/12408.png
## Source : http://tile.stamen.com/terrain/15/9543/12408.png
## Source : http://tile.stamen.com/terrain/15/9537/12409.png
## Source : http://tile.stamen.com/terrain/15/9538/12409.png
## Source : http://tile.stamen.com/terrain/15/9539/12409.png
## Source : http://tile.stamen.com/terrain/15/9540/12409.png
## Source : http://tile.stamen.com/terrain/15/9541/12409.png
## Source : http://tile.stamen.com/terrain/15/9542/12409.png
## Source : http://tile.stamen.com/terrain/15/9543/12409.png
## Source : http://tile.stamen.com/terrain/15/9537/12410.png
## Source : http://tile.stamen.com/terrain/15/9538/12410.png
## Source : http://tile.stamen.com/terrain/15/9539/12410.png
## Source : http://tile.stamen.com/terrain/15/9540/12410.png
## Source : http://tile.stamen.com/terrain/15/9541/12410.png
## Source : http://tile.stamen.com/terrain/15/9542/12410.png
## Source : http://tile.stamen.com/terrain/15/9543/12410.png

# Currently there are 4 hospitals in Philadelphia with an organ bank facility.
helipad = subset(dat, helipad == "Y")
qmplot(x,y,data=helipad,colour=I('white'),size=I(1), main = "Helipad Density In PA", xlab = "Longitude", ylab = "Latitude") +
stat_density2d(aes(x = x, y = y, fill=..level..), data=helipad,geom="polygon", alpha=0.2) +
scale_fill_gradient(low = "green", high = "red")
## Using zoom = 8...

# The map shows that the helipads are mostly concentrated around two regions in PA which happen to be also the biggest cities in the state(Pittsburgh and Philadelphia).