Adding density to a spatial map in R

Question

Objective: Adding density to a map

Tested Packages: ggplot, ggmap, sp, rgdal

Based on the rdgal package, I have a map of a specific state.

library(rgdal)
us<-getData('GADM', country='USA', level=2)  
WV <-subset(us, NAME_1=="West Virginia")
plot(WV)

I want to plot some prescription data on the map. Based on the geocode() method, specified in the following post: R convert zipcode or lat/long to county, I converted the county name to a set of latitude and longitude coordinates.

library(ggmap)
WV.opioid <- opioid.dat[opioid.dat$State.Name == "West Virginia",]
WV.coordinates <- geocode(as.character(WV.opioid$County.Name))
WV.opioid <- cbind(WV.opioid, WV.coordinates)
WV.opioid <- WV.opioid[which(!is.na(WV.opioid$lat)),]
rownames(WV.opioid) <- NULL

I have tried to plot one column of the dataset on the map:

 library(ggplot2)
 #opioid data from following link: https://www.datafiles.samhsa.gov/study-series/national-survey-drug-use-and-health-nsduh-nid13517
 opioid.map <- ggplot(WV.opioid, aes(x = lon, y = lat, 
                             fill = cut_number(Diff.Opioid.Prescription.Rates, 
 5))) +
 geom_polygon(color = "gray10", size = 0.2) +
 coord_equal() +
 viridis::scale_fill_viridis(discrete = TRUE) +
 labs(title = "Test",
 fill = NULL) +
 theme_void() +
 theme(legend.position = "bottom",
 panel.background = element_rect(fill = NA, colour = "#cccccc"))
 opioid.map

I am not sure how to fix this plot and would appreciate help. I have also tried to use the rdgal package and add my data to a SpatialPoints object but have not been successful.

Could you please add the packages you are using to your code? — Jan Boyer, Mar 15 '18 at 15:23
Your code is not reproducible. There isn't any package called `rdgal`. Where did `opioid.dat` come from? — Tung, Mar 15 '18 at 16:32

score 0 · Answer 1 · answered Mar 16 '18 at 18:14

It's unclear what spatial resolution the opioid data is at, but I would recommend finding a shapefile that fits it, rather than using lat/lon county centroids. This type of transformation isn't really correct if you don't correct for the modifiable areal unit problem.

If the opioid data is by zip code, there are shapefiles at that resolution in the tigris package.

I think the main problem is in this transformation you are doing in the second code chunk. Once you are no longer doing that, it should be fairly simple to plot a chloropleth of the prescription rates.

Adding density to a spatial map in R

1 Answers1