Alright, let's do this thing

A - Prepare data/A1 Re-rasterize Rodents.R

@@ -0,0 +1,28 @@
+library(fasterize)
+library(sf)
+library(rgdal)
+
+setwd('C:/Users/cjcar/Dropbox/PlagueWNA/Layers')
+setwd('./PRISM climate layers/')
+
+pasty <- function(x) {raster(paste(paste(paste('./',x,sep=''),x,sep='/'),'bil',sep='.'))}
+
+tmin <- list.files(pattern="tmin")
+tmin <- stack(lapply(tmin,pasty))
+blank <- tmin[[1]]*0
+
+mammal <- readOGR('C:/Users/cjcar/Downloads/TERRESTRIAL_MAMMALS/TERRESTRIAL_MAMMALS.shp')
+mammal <- mammal[mammal$order_=='RODENTIA',]
+mammal <- st_as_sf(mammal)
+
+mammal.r <- fasterize(mammal, blank, fun='count')
+mammal.r <- mammal.r + blank
+
+mammal.r <- crop(mammal.r, extent(-130, -95, 10, 49))
+
+# COME BACK TO THIS LINE
+# IF YOU WANT TO RE-EXTENT THE MAP
+
+writeRaster(mammal.r, 
+            'C:/Users/cjcar/Dropbox/PlagueWNA/Layers/RodentRichness.tif',
+            overwrite=TRUE)

A - Prepare data/A2 Soil macronutrient kriging.R

@@ -0,0 +1,63 @@
+
+layerA <- read.csv('C:/Users/cjcar/Dropbox/PlagueWNA/Appendix_3a_Ahorizon_18Sept2013.csv')
+
+layerA$A_Ca <- as.character(layerA$A_Ca)
+layerA$A_Fe <- as.character(layerA$A_Fe)
+layerA$A_Na <- as.character(layerA$A_Na)
+
+layerA$A_Ca[layerA$A_Ca == 'N.S.'] <- NA
+layerA$A_Ca[layerA$A_Ca == 'INS'] <- 0
+layerA$A_Ca[layerA$A_Ca == '<0.01'] <- 0
+layerA$A_Ca <- as.numeric(as.character(layerA$A_Ca))
+
+layerA$A_Fe[layerA$A_Fe == 'N.S.'] <- NA
+layerA$A_Fe[layerA$A_Fe == 'INS'] <- 0
+layerA$A_Fe[layerA$A_Fe == '<0.01'] <- 0
+layerA$A_Fe <- as.numeric(as.character(layerA$A_Fe))
+
+layerA$A_Na[layerA$A_Na == 'N.S.'] <- NA
+layerA$A_Na[layerA$A_Na == 'INS'] <- 0
+layerA$A_Na[layerA$A_Na == '<0.01'] <- 0
+layerA$A_Na <- as.numeric(as.character(layerA$A_Na))
+
+library(sp)
+soilpts <- SpatialPointsDataFrame(layerA[,c("Longitude", "Latitude")], 
+                               layerA[,c('A_Fe', 'A_Ca', 'A_Na')])
+
+library(raster)
+blank <- raster('C:/Users/cjcar/Dropbox/PlagueWNA/Layers/RodentRichness.tif')*0
+blank <- as(blank, "SpatialPixelsDataFrame")
+crs(blank) <- NA
+
+library(automap)
+
+#Perform ordinary kriging and store results inside object of type "autoKrige" "list" 
+kriging_result = autoKrige(A_Fe~1, soilpts[complete.cases(soilpts$A_Fe),], blank)
+plot(kriging_result)
+
+kresult <- raster(kriging_result$krige_output)
+writeRaster(kresult, 'C:/Users/cjcar/Dropbox/PlagueWNA/kriged-iron.tif',
+            overwrite=TRUE)
+
+
+
+
+#Perform ordinary kriging and store results inside object of type "autoKrige" "list" 
+kriging_result = autoKrige(A_Ca~1, soilpts[complete.cases(soilpts$A_Ca),], blank)
+plot(kriging_result)
+
+kresult <- raster(kriging_result$krige_output)
+writeRaster(kresult, 'C:/Users/cjcar/Dropbox/PlagueWNA/kriged-calcium.tif',
+            overwrite=TRUE)
+
+
+
+
+#Perform ordinary kriging and store results inside object of type "autoKrige" "list" 
+kriging_result = autoKrige(A_Na~1, soilpts[complete.cases(soilpts$A_Na),], blank)
+plot(kriging_result)
+
+kresult <- raster(kriging_result$krige_output)
+writeRaster(kresult, 'C:/Users/cjcar/Dropbox/PlagueWNA/kriged-sodium.tif',
+            overwrite=TRUE)
+

A - Prepare data/A3 SoilGrids resampling.R

@@ -0,0 +1,34 @@
+
+setwd('C:/Users/cjcar/Dropbox/PlagueWNA/Layers/SoilGrids')
+
+clippy <- 0*raster('C:/Users/cjcar/Dropbox/PlagueWNA/Layers/RodentRichness.tif')
+
+clay <- raster('CLYPPT_M_sl1_250m.tif')
+clay.clip <- crop(clay, clippy)
+clay.sample <- resample(clay.clip, clippy)
+clay.sample <- sum(clay.sample,clippy,na.rm=FALSE)
+writeRaster(clay.sample,'clay-resample.tif',overwrite=TRUE)
+
+org <- raster('ORCDRC_M_sl1_250m.tif')
+org.clip <- crop(org, clippy)
+org.sample <- resample(org.clip, clippy)
+org.sample <- sum(org.sample,clippy,na.rm=FALSE)
+writeRaster(org.sample,'organic-resample.tif',overwrite=TRUE)
+
+snd <- raster('SNDPPT_M_sl1_250m.tif')
+snd.clip <- crop(snd, clippy)
+snd.sample <- resample(snd.clip, clippy)
+snd.sample <- sum(snd.sample,clippy,na.rm=FALSE)
+writeRaster(snd.sample,'sand-resample.tif',overwrite=TRUE)
+
+phi <- raster('PHIHOX_M_sl1_250m.tif')
+phi.clip <- crop(phi, clippy)
+phi.sample <- resample(phi.clip, clippy)
+phi.sample <- sum(phi.sample,clippy,na.rm=FALSE)
+writeRaster(phi.sample,'ph-resample.tif',overwrite=TRUE)
+
+chex <- raster('CECSOL_M_sl1_250m.tif')
+chex.clip <- crop(chex, clippy)
+chex.sample <- resample(chex.clip, clippy)
+chex.sample <- sum(chex.sample,clippy,na.rm=FALSE)
+writeRaster(chex.sample,'CEC-resample.tif',overwrite=TRUE)

B - Extract and generate model datasets/B1 Process and load environmental layers.R

@@ -0,0 +1,112 @@
+
+library(devtools)
+library(elevatr)
+library(raster)
+library(readr)
+library(rgeos)
+#library(rgdal) 
+library(fasterize)
+library(sf)
+
+#######################################################################
+#1. GENERATE EACH COMPONENT LAYER
+
+# Mammals
+
+setwd('C:/Users/cjcar/Dropbox/PlagueWNA/Layers/')
+rodents <- raster('RodentRichness.tif')
+rodents <- trim(rodents)
+
+blank <- rodents*0
+
+# Climate layers
+
+setwd('./PRISM climate layers/')
+
+pasty <- function(x) {raster(paste(paste(paste('./',x,sep=''),x,sep='/'),'bil',sep='.'))}
+
+tmin <- list.files(pattern="tmin")
+tmin <- stack(lapply(tmin,pasty))
+
+tmean <- list.files(pattern="tmean")
+tmean <- stack(lapply(tmean,pasty))
+
+tmax <- list.files(pattern="tmax")
+tmax <- stack(lapply(tmax,pasty))
+
+ppt <- list.files(pattern="ppt")
+ppt <- stack(lapply(ppt,pasty))
+
+tmin <- projectRaster(tmin, blank)
+tmin <- tmin+blank
+
+tmean <- projectRaster(tmean, blank)
+tmean <- tmean+blank
+
+tmax <- projectRaster(tmax, blank)
+tmax <- tmax+blank
+
+ppt <- projectRaster(ppt, blank)
+ppt <- ppt+blank
+
+# Some new code for generating anomalies
+
+tmean.norm <- (tmean-mean(tmean))/calc(tmean,var)
+ppt.norm <- (ppt-mean(ppt))/calc(ppt,var)
+
+# SoilGrids and other soil
+
+setwd('C:/Users/cjcar/Dropbox/PlagueWNA/Layers/SoilGridsCLayer')
+# ^ If it's just SoilGrids, it'll point to A Layer.
+
+soil.layers <- stack(raster('CEC-resample.tif'),
+                     raster('ph-resample.tif'),
+                     raster('sand-resample.tif'),
+                     raster('clay-resample.tif'),
+                     raster('organic-resample.tif'))
+soil.layers <- trim(soil.layers)
+
+iron <- raster('C:/Users/cjcar/Dropbox/PlagueWNA/kriged-iron.tif') + rodents*0
+calcium <- raster('C:/Users/cjcar/Dropbox/PlagueWNA/kriged-calcium.tif') + rodents*0
+sodium <- raster('C:/Users/cjcar/Dropbox/PlagueWNA/kriged-sodium.tif') + rodents*0
+
+soil.layers <- stack(soil.layers,iron)
+soil.layers <- stack(soil.layers,calcium)
+soil.layers <- stack(soil.layers,sodium)
+
+soil.layers@crs <- rodents@crs
+
+# Elevation
+
+library(elevatr)
+elev <- elevatr::get_elev_raster(locations = rodents, z=6)
+elev <- resample(elev,rodents)
+elev <- elev + rodents*0
+elev <- projectRaster(elev,soil.layers)
+
+#######################################################################
+#2. GENERATE AN OBJECT OF ANNUAL LAYERS
+
+rodents@crs <- soil.layers@crs
+
+stack.year <- function(yr) {
+
+  stack.fixed <- stack(rodents,soil.layers,elev)
+  year.guys <- stack(tmin[[yr-1949]],
+                     tmean[[yr-1949]],
+                     tmax[[yr-1949]],
+                     ppt[[yr-1949]],
+                     ppt.norm[[yr-1949]],
+                     tmean.norm[[yr-1949]])
+  pred <- stack(stack.fixed,year.guys)
+
+  names(pred) <- c('rodent','cec','ph','sand','clay','orgc','Fe','Ca','Na','elev',
+                   'tmin','tmean','tmax','ppt','ppt.n','tmp.n')
+  return(pred)
+
+}
+
+stackmaker <- lapply(c(1950:2017), function(x){print('one down')
+  stack.year(x)})
+
+write_rds(stackmaker, '~/Github/meridian100/Pipeline2020/stackmaker.RDS')

B - Extract and generate model datasets/B2 Extract NWDP data.R

@@ -0,0 +1,47 @@
+
+#library(gbm)
+library(dismo)
+library(dbarts)
+library(embarcadero)
+library(velox)
+
+library(tidyverse)
+nwdp <- read_csv('C:/Users/cjcar/Desktop/plague data download 1May2019.csv')
+
+nwdp$CollectionDate <- lubridate::dmy(nwdp$CollectionDate)
+nwdp$year <- lubridate::year(nwdp$CollectionDate)
+
+for(i in 2000:2017) {
+
+  nwdp %>% filter(year==i) -> sub 
+  if(nrow(sub)>0) {
+    env <- stackmaker[[i-1949]]
+
+    pres <- sub[sub$PlagueResults=='Positive',]
+    abs <- sub[sub$PlagueResults=='Negative',]
+
+    p.pts <- data.frame(cbind(1,raster::extract(env,pres[,c('Longitude','Latitude')])))
+    a.pts <- data.frame(cbind(0,raster::extract(env,abs[,c('Longitude','Latitude')])))
+    data <- rbind(p.pts,a.pts)
+    data$yr <- i
+
+    if(i==2000){big <- data} else {big <- rbind(big,data)}
+    print(i)
+  }
+}
+
+big <- data.frame(big)
+
+###############
+
+colnames(big) <- c('case','rodent','cec','ph','sand','clay','orgc','Fe','Ca','Na','elev',
+                   'tmin','tmean','tmax','ppt','ppt.n','tmp.n','yr')
+
+set.seed(69)
+big <- rbind(big %>% filter(case==0) %>% sample_n(nrow(big %>% filter(case==1))),
+             big %>% filter(case==1))
+
+big <- na.omit(big)
+
+write.csv(big, '~/Github/meridian100/Pipeline2020/NWDPdata.csv')
+

B - Extract and generate model datasets/B3 Extract human data.R

@@ -0,0 +1,57 @@
+library(sp)
+library(rgdal)
+library(raster)
+
+setwd('C:/Users/cjcar/Documents/Github/meridian100/Raw data')
+point.data <- read.csv('western USA-edit.csv')
+
+counties <- readOGR(layer='CountiesNoHole',dsn='.')
+counties[490,grepl('YEAR', colnames(counties@data))] <- 0
+counties[490,'TOTAL'] <- 0
+
+sp.point <- SpatialPointsDataFrame(point.data[,1:2], point.data[,3:59])
+sp.point@proj4string <- counties@proj4string
+
+j = 1 
+set.seed(10) # Added this - worth a double check, threw an error once without but stochastically?
+for (year in 1950:2005) {
+  for (co in 1:nrow(counties@data)) {
+    if(counties@data[co,year-1940]>0) {
+      s.sub <- spsample(counties[co,],n=counties@data[co,year-1940], "random")
+      s.sub <- SpatialPointsDataFrame(s.sub@coords,data=data.frame(year=rep(year,counties@data[co,year-1940])))
+      if (j == 1) { s <- s.sub} else {s <- bind(s, s.sub)}
+      j <- j+1
+    }
+  }}
+
+#plot(counties,col=(counties$TOTAL>0))
+#points(s,col='red',pch=16)
+
+county.raster <- fasterize::fasterize(sf::st_as_sf(readOGR(layer='CountiesNoHole',dsn='.')), rodents)
+
+for(i in 1950:2005) {
+
+  sub <- s[s$year==i,]
+  abs <- randomPoints(county.raster, n=7)
+  env <- stackmaker[[i-1949]]
+
+  p.pts <- data.frame(cbind(1,raster::extract(env,sub)))
+  a.pts <- data.frame(cbind(0,raster::extract(env,abs)))
+  data <- rbind(p.pts,a.pts)
+  data$yr <- i
+
+  if(i==1950){big <- data} else {big <- rbind(big,data)}
+  print(i)
+}
+
+big <- data.frame(big)
+
+colnames(big) <- c('case','rodent','cec','ph','sand','clay','orgc','Fe','Ca','Na','elev',
+                   'tmin','tmean','tmax','ppt','ppt.n','tmp.n','yr')
+
+#set.seed(69)
+#big <- rbind(big %>% filter(case==0) %>% sample_n(nrow(big %>% filter(case==1))),
+#             big %>% filter(case==1))
+# ^ this may have been flipped off at some point, strategically - maybe reverse if results get weird
+
+write.csv(big, 'C:/Users/cjcar/Documents/GitHub/meridian100/Pipeline2020/Humandata.csv')