ANLY 512 - Problem Set 5

Spatial Visualization

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)
library(ggmap)

## Loading required package: ggplot2

## Warning: package 'ggplot2' was built under R version 3.6.2

## 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(sf)

## Warning: package 'sf' was built under R version 3.6.2

## Linking to GEOS 3.7.2, GDAL 2.4.2, PROJ 5.2.0

library(rgdal)

## Warning: package 'rgdal' was built under R version 3.6.2

## Loading required package: sp

## Warning: package 'sp' was built under R version 3.6.2

## rgdal: version: 1.5-16, (SVN revision 1050)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 2.4.2, released 2019/06/28
## Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.6/Resources/library/rgdal/gdal
## GDAL binary built with GEOS: FALSE 
## Loaded PROJ runtime: Rel. 5.2.0, September 15th, 2018, [PJ_VERSION: 520]
## Path to PROJ shared files: /Library/Frameworks/R.framework/Versions/3.6/Resources/library/rgdal/proj
## Linking to sp version:1.4-2
## Overwritten PROJ_LIB was /Library/Frameworks/R.framework/Versions/3.6/Resources/library/rgdal/proj

library(dplyr)

## Warning: package 'dplyr' was built under R version 3.6.2

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(ggrepel)
library(reshape2)

## Warning: package 'reshape2' was built under R version 3.6.2

library(plotly)

## Warning: package 'plotly' was built under R version 3.6.2

## 
## Attaching package: 'plotly'

## The following object is masked from 'package:ggmap':
## 
##     wind

## The following object is masked from 'package:ggplot2':
## 
##     last_plot

## The following object is masked from 'package:stats':
## 
##     filter

## The following object is masked from 'package:graphics':
## 
##     layout

dat <- read.dbf("pennsylv.dbf")  # change this to connect to your version
latlong <- dat[,c(117, 118)]

# county <- read.dbf("tl_2017_us_county.dbf")
# county$INTPTLAT <- as.numeric(as.character(county$INTPTLAT))
# county$INTPTLON <- as.numeric(as.character(county$INTPTLON))
# countylatlon <- county[,c(16,17)]
# 
# penncounties <- subset(county, county$STATEFP == 42)
# 
# penncountylatlon <-penncounties[,c(16,17)]

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

# register_google(key = "XX-XX", write = TRUE)
penn <- geocode("pennsylvania", source = "google")

## Source : https://maps.googleapis.com/maps/api/geocode/json?address=pennsylvania&key=xxx-lTHO5gsOu6t9X6bOFjyoIl8S6I

centerpenn <- as.numeric(penn)

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

Map 1 - Hospitals in Pennsylvania

map1 <- ggmap(get_googlemap(center = centerpenn, zoom = 7, maptype = "terrain")) #markers = df, 

## Source : https://maps.googleapis.com/maps/api/staticmap?center=41.203322,-77.194525&zoom=7&size=640x640&scale=2&maptype=terrain&key=xxx-lTHO5gsOu6t9X6bOFjyoIl8S6I

map1 +
  geom_point(aes(x=x,y=y), data=latlong, alpha=0.5) +
  labs(x="Longitude",y="Latitude",title="Hospitals in Pennsylvania")

Map 2 - Pennsylvania & Hospitals in Pennsylvania

counties <- readOGR(dsn="tl_2017_us_county", layer="tl_2017_us_county")

## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/anjanasasi/Documents/HarrisburgU/OneDrive - Harrisburg University/0ANLY512 Data Visualization/Assisnments-ProblemSets/ProblemSet5/tl_2017_us_county", layer: "tl_2017_us_county"
## with 3233 features
## It has 17 fields
## Integer64 fields read as strings:  ALAND AWATER

counties$GEOID <- as.numeric(as.character(counties$GEOID))
counties$INTPTLAT <- as.numeric(as.character(counties$INTPTLAT))
counties$INTPTLON <- as.numeric(as.character(counties$INTPTLON))
PAcounties <- counties[which(counties$STATEFP == "42"),]
pac <- fortify(PAcounties)

## Regions defined for each Polygons

ggplot() + geom_polygon(data=pac, aes(x=long,y=lat, group=group)) +
  geom_point(aes(x=x,y=y), data=latlong, alpha=0.7, size=3, color="orange", position = "jitter") +
  labs(x="Longitude",y="Latitude",title="Hospitals in Pennsylvania")

Map 3 - Allegheny County and Hospitals

agh_sf <- subset(PAcounties, NAME=="Allegheny")
agh <- fortify(agh_sf)

## Regions defined for each Polygons

aghHosp <- subset(dat, dat$county == "Allegheny")

map2 <- ggmap(get_googlemap("Allegheny, Pennsylvania", zoom = 10, maptype = "roadmap")) 

## Source : https://maps.googleapis.com/maps/api/staticmap?center=Allegheny,%20Pennsylvania&zoom=10&size=640x640&scale=2&maptype=roadmap&key=xxx-lTHO5gsOu6t9X6bOFjyoIl8S6I

## Source : https://maps.googleapis.com/maps/api/geocode/json?address=Allegheny,+Pennsylvania&key=xxx-lTHO5gsOu6t9X6bOFjyoIl8S6I

map2 + 
  geom_polygon(data=agh, aes(x=long,y=lat, group=group), alpha=0.3) +
  geom_point(aes(x=x,y=y), data=aghHosp, alpha=0.5, size=3, color="#E21313") +
  labs(x="Longitude",y="Latitude",title = "Hospitals in Allegheny County")

Map 4 - Density of Hospitals in Pennsylvania

dataa <- dat %>%
  group_by(county) %>%
  summarise(zip = n(), x=mean(x), y=mean(y))

## `summarise()` ungrouping output (override with `.groups` argument)

ggplot() + 
  geom_polygon(data=pac, aes(x=long,y=lat, group=group), fill="white", color="grey") +
  geom_point(aes(x=x,y=y, size=zip), data=dataa, alpha=0.5, color="darkgreen") +
  scale_size(range = c(1,15)) + 
  labs(x="Longitude",y="Latitude",size="# of hospitals",title="Density of Hospitals in Pennsylvania") +
  geom_text_repel(data = dataa, aes(x=x,y=y,label = county), size=2.5) +
  theme(panel.background = element_blank())

Map 5 - Dental and Heart Transplant departments and Organ Banks

depts <- dat[,c(26,32,53,94,117,118)]
colnames(depts) <- c("County","Dental","Heart Transplant","Organ Bank","x","y")

departments <- melt(depts, id.vars = c("County","x","y"), measure.vars=c("Dental","Heart Transplant", "Organ Bank"))
departments <- subset(departments, value == "Y")

map5 <- ggmap(get_googlemap(center = centerpenn, zoom = 7, maptype = "terrain")) #markers = df, 

## Source : https://maps.googleapis.com/maps/api/staticmap?center=41.203322,-77.194525&zoom=7&size=640x640&scale=2&maptype=terrain&key=xxx-lTHO5gsOu6t9X6bOFjyoIl8S6I

map5 +
  geom_point(aes(x=x,y=y), data=departments, alpha=0.5, size=3, color="orange", position = "jitter") +
  labs(x="Longitude",y="Latitude") + facet_wrap(~ variable) +
  geom_polygon(data=pac, aes(x=long,y=lat, group=group), alpha=0.2)