Addition of two functions

DESCRIPTION

@@ -13,5 +13,13 @@ Depends:
     raster,
     proj4,
     sp,
-    tidyr
+    tidyr,
+    dplyr,
+    data.table
+Suggests:
+    ggplot2,
+    parallel,
+    sf,
+    viridis,
+    scales
 RoxygenNote: 6.0.1

R/rinmap.R

@@ -194,3 +194,103 @@ create_input_facility_all_but_one <- function(input_csv, path = basename(input_c
   }
 }
 
+#' Combine zipcode-linked inMAP output into single data object and link with spatial
+#' data for plotting.
+#' 
+#' \code{combine_inmap_output} uses spatial linkages of an inMAP outputs to combine output
+#' multiple runs. Employs the sf package to link with spatial data if it is installed.
+#' 
+#' @param path.out Directory that houses InMAP output csv files
+#' @param zcta_shapefile A ZCTA shapefile
+#' @param pattern A text or regex pattern common to all InMAP output you wish to be joined
+#' 
+#' @return A data frame of average PM2.5 values at the ZIP code level.
+combine_inmap_output <- function(path.out, pattern = NULL){
+
+	#list files for import, read in files
+	files = list.files(path.out,full.names=T)
+	names.f = gsub(paste(pattern,'_linked_zip.csv',sep='|'),'', list.files(path.out,full.names=F))
+	names(files) <- names.f
+	if (!is.null(pattern)) files = files[grep(pattern,files)]
+	if (length(files) == 0) stop(("No matching files in path.out!"))
+
+	im <- lapply(seq_along(files), function(x,f,n) {fin <- fread(f[x])[, V1 := NULL]
+							setnames(fin, 'PM25inMAP', n[x])
+							return(fin)}, 
+		     					files, names(files))
+
+	#reduce list to single data table
+	im <- Reduce(function(...) merge(..., all = TRUE, by = c("ZCTA","ZIP","PO_NAME","STATE","ZIP_TYPE")), im)
+
+	#convert a ZIP code from 3-digit to 5-digit format
+	im$ZIP <- formatC(im$ZIP, width = 5, format = "d", flag = "0")
+
+	return(im)
+}
+
+
+#' Using zipcode-linked inMAP output from \code{combine_inmap_output}, plot change in InMAP impacts
+#' at zip code level compared to base year.
+#' 
+#' \code{combine_inmap_output} uses spatial linkages of an inMAP outputs to combine output
+#' multiple runs. Requires the sf, ggplot2, and parallel packages.
+#' 
+#' @param read_inmap_d Directory that houses InMAP output csv files
+#' @param legend_lims Legend limits for zip code fill
+#' @param path.plot Output directory for plots. If it does not exist, it will be created
+#' @param cores Cores available to create plots. Defaults to one
+#' 
+#' @return A list of ggplot objects.
+plot_inmap <- function(read_inmap_d,
+   		       zcta_shapefile = "~/shared_space/ci3_nsaph/software/inmap/zcta/cb_2015_us_zcta510_500k.shp",
+		       legend_lims=c(-5,5),path.plot='InMAP_plots',cores=1){
+	#check if required packages are installed
+	if(F %in% (c('sf','parallel','ggplot2','viridis','scales') %in% (.packages()))) {stop("Required package missing! (need sf,parallel,ggplot2,viridis,scales)")}
+
+	#create directory if it does not exist
+	if (!file.exists(path.plot)) dir.create(path.plot)
+
+	#join with spatial zip data
+	zips <- st_read(zcta_shapefile)
+	zips$GEOID10 <- as.character(zips$GEOID10)
+	im_j <- left_join(im,zips,by=c('ZIP' = 'GEOID10'))
+
+	#extract names from combined data table/sf object
+	names.f <- names(read_inmap_d)[-grep(c('ZCTA|ZIP|PO_NAME|STATE|ZIP_TYPE|ZCTA5CE10|AFFGEOID10|ALAND10|AWATER10|geometry'),names(read_inmap_d))]
+
+	#read in USA and state shapes
+	cl <- makeCluster(cores)
+	clusterExport(cl,c( 'data.table','ggplot','aes','theme_bw','geom_sf','map_data',
+						'labs','geom_polygon','coord_sf','scale_color_viridis',
+						'scale_fill_viridis','theme','element_text','element_rect',
+						'unit','element_blank','ggsave','setnames','squish'))
+
+	#create plotting object
+	ggplotter <- function(x,im_j,n,ll){
+		x1 <- data.table(im_j)[,c('ZIP',n[x],'geometry'),with=F]
+		setnames(x1,n[x],'PM')
+		usa.states <- map_data("state")
+
+		gg <- ggplot(data = x1, aes(fill = -PM,color= -PM)) + 
+					theme_bw() + labs(title=paste('InMAP exposure change - ',n[x],sep='')) +
+					geom_sf(size=0.05) +
+					geom_polygon(data= usa.states, aes(x = long, y = lat, group = group),
+								 fill = NA, colour = "grey50",size=.25) + 
+					coord_sf(xlim = c(-123,-69),ylim=c(24,50),datum=NA) + 
+					scale_color_viridis(discrete=F,option='D',limits= ll,oob=squish,direction=-1) +
+					scale_fill_viridis(discrete=F,option='D',limits= ll,oob=squish,direction=-1) + 
+					theme(legend.position = c(.25,.15),	axis.text=element_blank(), 
+						axis.title.x = element_blank(),axis.title.y = element_blank(),legend.title=element_blank(),
+						axis.title=element_text(size=24),legend.text=element_text(size = 14),
+						legend.background=element_rect(fill='transparent'),strip.text=element_text(size=16),
+						legend.key.size = unit(.05,'npc'),plot.title = element_text(size = 16,hjust = 0.5),
+						legend.direction = 'horizontal',strip.background=element_rect(fill='white'))
+		invisible(ggsave(file.path('.',path.plot,paste('plot_',n[x],'.png',sep='')), 	
+						gg,width=13.5,height=7.79,unit='in'))
+		return(gg)
+	}
+	out <- parLapply(cl,seq_along(names.f),ggplotter, read_inmap_d,names.f,ll=legend_lims) #
+	stopCluster(cl)
+	return(out)
+}
+