knitr::opts_chunk$set(echo = TRUE)
library(maps)
library(leaflet)
library(foreign)
library(ggthemes)
library(ggplot2)
library(raster)
## Loading required package: sp
library(shapefiles)
## 
## Attaching package: 'shapefiles'
## The following objects are masked from 'package:foreign':
## 
##     read.dbf, write.dbf
library(rgdal)
## rgdal: version: 1.2-16, (SVN revision 701)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.2.0, released 2017/04/28
##  Path to GDAL shared files: C:/Users/prerna/Documents/R/win-library/3.3/rgdal/gdal
##  GDAL binary built with GEOS: TRUE 
##  Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
##  Path to PROJ.4 shared files: C:/Users/prerna/Documents/R/win-library/3.3/rgdal/proj
##  Linking to sp version: 1.2-6
library(ggmap)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:raster':
## 
##     intersect, select, union
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
getwd()
## [1] "D:/MyRWork/Data512/Hospital"
setwd("D:/MyRWork/Data512/Hospital")
HP <- shapefile("DOH_Hospitals201709.shp")
HP1 <- data.frame(HP)
HPmap <- ggmap(get_map("Pennsylvania",maptype = "roadmap",source = "google", zoom = 7), extent = "device", legend = "top")
## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=Pennsylvania&zoom=7&size=640x640&scale=2&maptype=roadmap&language=en-EN&sensor=false
## Information from URL : http://maps.googleapis.com/maps/api/geocode/json?address=Pennsylvania&sensor=false
## Warning: `panel.margin` is deprecated. Please use `panel.spacing` property
## instead

MAP 1: Plot Hospitals in Pennsylvania

In this map we have retrieved all the hospitals in Pennsylvania. It shows the total spread of hospitals in the state.

HPmap +
geom_point(data = HP1 , aes(x = LONGITUDE, y = LATITUDE), col = HP1$AREA_CODE, size = 1 ) + labs(title = "Hospitals in  Pennsylvania") + 
  theme_map()

MAP 2 Map tells us the high concentration of hospitals in particular areas which are Philadelphia, Pittsburg, Harrisburg and Allentown

ggmap(get_map("Pennsylvania", maptype = "roadmap", zoom = 7, source = "google"), legend = "topright") + stat_bin2d( aes(x = LONGITUDE, y = LATITUDE, fill = cut(..count.., c(0, 3, 6, 9, 12, Inf))), size = 1,  alpha = 0.8, 
data = HP1) + labs(fill = "No. of Hospitals in Area", title = "Hospitals in Pennsylvania Area" ) + theme()
## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=Pennsylvania&zoom=7&size=640x640&scale=2&maptype=roadmap&language=en-EN&sensor=false
## Information from URL : http://maps.googleapis.com/maps/api/geocode/json?address=Pennsylvania&sensor=false

MAP 3

Density of hospitals is high in Philadelphia area. Buck county has a bigger cluster of hospitals.

HPmap +
  geom_density2d( aes(x = LONGITUDE, y = LATITUDE, alpha = ..level..) ,size =2, bins = 5, data = HP1, color = "maroon") +labs(title = "Density map for Hospitals in State of Pennsylvania") + theme_map()

MAP 4

Below are the top 6 cities with high density of hospitals

penn_city <- head(HP1 %>% group_by(CITY_TOWN) %>% count() %>% arrange(desc(n)), 6)
penncitydf = subset(HP1, CITY_TOWN %in% penn_city$CITY_TOWN)

HPmap +
geom_point(data = penncitydf, aes(x = LONGITUDE, y = LATITUDE), colour = "maroon") + 
 facet_wrap(~ CITY_TOWN) +
 labs(title = "Top 6 cities in State of Pennsylvania by number of hospitals") +
  theme_map()

MAP 5

Below maps tells us 7 zipcodes with higher number of hospitals in Pennsylvania. We can see the namnes of the hospitals via popup’s.

Pennhigherdensityzip <- head(HP1 %>% group_by(ZIP_CODE) %>% count() %>% arrange(desc(n)), 7)

Pennhigherdensityzip_df = subset(HP1, ZIP_CODE %in% Pennhigherdensityzip$ZIP_CODE)

leaflet() %>%
addTiles() %>%
addMarkers(data = Pennhigherdensityzip_df, ~LONGITUDE, ~LATITUDE) %>%
addPopups(Pennhigherdensityzip_df$LONGITUDE, Pennhigherdensityzip_df$LATITUDE,  Pennhigherdensityzip_df$FACILITY_N)