Data Preparation

Data source is a list of all Hospitals that have been registered with Medicare. See website: https://data.medicare.gov/Hospital-Compare/Hospital-General-Information/xubh-q36u

Get the latitude and longitude for the ZIP codes.

So let’s start with getting a count of hospital by ZIP, then merging our ZIP codes data with those counts. Make sure to have the packages at the top installed.

library(zipcode)
library(ggplot2)
library(maps)
library(viridis)
## Loading required package: viridisLite
library(viridisLite)
#data
Hos<-read.csv("C:/Work/R/Google density map/Hospital in US/Hospital_General_Information.csv",stringsAsFactors = F)
data(zipcode)
Hos$ZIPCode<- clean.zipcodes(Hos$ZIPCode)
#size by zip
Hos.ZIPCode<-aggregate(data.frame(count=Hos$ProviderID),list(zip=Hos$ZIPCode,county=Hos$CountyName),length)
Hos<- merge(Hos.ZIPCode, zipcode, by='zip')

Now we have a dataframe that contains the ZIP code, a count of the hospitals, the city, state latitude and longitude.So we want to map this, to show all of our hospitals by ZIP, and we will use the viridis package for the color scale.

us<-map_data('state')

ggplot(Hos,aes(longitude,latitude)) +
  geom_polygon(data=us,aes(x=long,y=lat,group=group),color='gray',fill=NA,alpha=.35)+
  geom_point(aes(color = count),size=1,alpha=0.25) + scale_colour_gradient(high="dark blue", low="light blue") +
  theme_classic()+
  xlim(-125,-65)+ylim(20,50)
## Warning: Removed 86 rows containing missing values (geom_point).

#counties
county_df <- map_data("county")
names(county_df) <- c("long", "lat", "group", "order", "state_name", "county")
county_df$state <- state.abb[match(county_df$state_name, tolower(state.name))]
county_df$state_name <- NULL
Hos$county<-tolower(Hos$county)
state_df <- map_data("state")

choropleth <- merge(county_df, Hos, by = c("state", "county"))
choropleth <- choropleth[order(choropleth$order), ]

Upgrade to Heatmap

#
ggplot(choropleth, aes(long, lat, group = group)) +
  geom_polygon(data=county_df,aes(long,lat,group=group),fill='#ffffb3',colour = alpha("gray", 1 / 2))+
  geom_polygon(aes(fill = count), colour = alpha("gray", 1 / 2), size = 0.2) +
  geom_polygon(data = state_df, colour = "white", fill = NA) +
  theme(axis.line = element_blank(), axis.text = element_blank(),
        axis.ticks = element_blank(), axis.title = element_blank()) +
  scale_fill_gradientn(colours = rev(rainbow(7)),
                         breaks = c(0:6),
                         trans = "log10") +theme_minimal()

 # scale_fill_viridis(option="magma",direction=-1)+theme_minimal()