This R Markdown Document Will Serve as a Demonstration the code Used for my Term Project “An Attempt to Estimate 2021 West Nile Incident Rates in Michigan Peninsula to Validate Results from the Disaggregation R Package”

load packages

#note that INLA is not on cran. INLA needs to be downloaded as instructed here: https://www.r-inla.org/download-install
library(SpatialEpi)
library(dplyr)
library(disaggregation)
library(ggplot2)
library(sf)
library(terra)
library(tidyverse)
library(raster)
library(INLA)

prepare and plot the West Nile cases for 2021 by Michigan county

#read in 2021 Michigan Counties by Human Incident Rate
micounty2021 <- st_read("MIShapefile2021.shp")
## Reading layer `MIShapefile2021' from data source 
##   `D:\anderson\OneDrive - W.E. Upjohn Institute\Desktop\PSU\Geog588\Final Project\FinalProject\MIShapefile2021.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 83 features and 27 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -90.41814 ymin: 41.69612 xmax: -82.41347 ymax: 48.2388
## Geodetic CRS:  NAD83
#remove 0s and missing values from the 2021 Michigan County data
micounty2021_clean <- micounty2021[micounty2021$ReportedHu != 0, ]

#plot & save 2021 Michigan Counties cleaned data by Human Incident Rate 
ggplot() +
  geom_sf(data = micounty2021, fill = NA) +
  geom_sf(data = micounty2021_clean, aes(fill = ReportedHu)) +
  labs(fill = "Cases") +
  theme_classic() +
  ggtitle("2021 Reported Human Case by County") +
  labs(caption = "Source: CDC 2023")
Fig 1. 2021 Michigan West Nile Cases by County.

Fig 1. 2021 Michigan West Nile Cases by County. 

Cases2021 <- ggplot() +
                geom_sf(data = micounty2021, fill = NA) +
                geom_sf(data = micounty2021_clean, aes(fill = ReportedHu)) +
                labs(fill = "Cases") +
                theme_classic() +
                ggtitle("2021 Reported Human Case by County") +
                labs(caption = "Source: CDC 2023")

ggsave("Cases2021.png", plot = Cases2021, dpi = 300)

prepare and plot the West Nile cases for 2020 by Michigan county

#read in Michigan Counties by Human Incident Rate
micounty2020 <- st_read("MIShapefile2020.shp")
## Reading layer `MIShapefile2020' from data source 
##   `D:\anderson\OneDrive - W.E. Upjohn Institute\Desktop\PSU\Geog588\Final Project\FinalProject\MIShapefile2020.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 83 features and 26 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -90.41814 ymin: 41.69612 xmax: -82.41347 ymax: 48.2388
## Geodetic CRS:  NAD83
#remove 0s and missing values from the Michigan County data
micounty2020_clean <- micounty2020[micounty2020$Reported_h != 0, ]

#plot Michigan Counties cleaned data by Human Incident Rate 
ggplot() +
  geom_sf(data = micounty2020, fill = NA) +
  geom_sf(data = micounty2020_clean, aes(fill = Reported_h)) +
  scale_fill_continuous(breaks = c(5, 10)) +
  labs(fill = "Cases") +
  theme_classic() +
  ggtitle("2020 Reported Human Case by County") +
  labs(caption = "Source: CDC 2023")
Fig 2. 2020 Michigan West Nile Cases by County.

Fig 2. 2020 Michigan West Nile Cases by County. 

Cases2020 <- ggplot() +
  geom_sf(data = micounty2020, fill = NA) +
  geom_sf(data = micounty2020_clean, aes(fill = Reported_h)) +
  scale_fill_continuous(breaks = c(5, 10)) +
  labs(fill = "Cases") +
  theme_classic() +
  ggtitle("2020 Reported Human Case by County") +
  labs(caption = "Source: CDC 2023")

ggsave("Cases2020.png", plot = Cases2020, dpi = 300)

#buffer the ploygon geometry to account for small inconsistencies
micounty2020_clean <- sf::st_buffer(micounty2020_clean, dist = 0)

load in LST, EVI, and POP rasters, adjusting LST and EVI rasters

#load in raster files
rasterEVI <- raster("rasters/EVI.tif")
rasterLST <- raster("rasters/LST.tif")
rasterPOP <- raster("usa_ppp_2020_constrained.tif")

#find extents of EVI and LST rasters
extEVI <- extent(rasterEVI)
extLST <- extent(rasterLST)
print(extEVI)
## class      : Extent 
## xmin       : -8417414 
## xmax       : -5443494 
## ymin       : 4636402 
## ymax       : 5363803
print(extLST)
## class      : Extent 
## xmin       : -8427013 
## xmax       : -5443484 
## ymin       : 4636402 
## ymax       : 5364730
#get max extent of both EVI and LST rasters
max_extent <- raster::union(extent(rasterEVI), extent(rasterLST))

#set EVI and LST raster extents to be the same 
rasterEVI_max <- setExtent(rasterEVI, max_extent)
rasterLST_max <- setExtent(rasterLST, max_extent)

#set EVI and LST rasters to same extent
EVIresample <- projectRaster(rasterEVI_max, rasterEVI_max, method = 'ngb')
LSTresample <- projectRaster(rasterLST_max, rasterEVI_max, method = 'ngb')

#set target crs and project each resampled raster to be in the same extent as the MI County data
target_crs <- "+proj=longlat +datum=NAD83 +no_defs"
EVIresample_NAD83 <- projectRaster(EVIresample, crs = target_crs)
LSTresample_NAD83 <- projectRaster(LSTresample, crs = target_crs)

plot and save EVI raster

plot(EVIresample_NAD83,
     col = terrain.colors(255),
     main = "EVI Raster")

mtext("Source: NASA; MODIS enhanced vegetation index (EVI)", side = 1, line = 4)
Fig 3. EVI Raster.

Fig 3. EVI Raster. 

#turn on png
png("EVI_plot.png", width = 800, height = 600)

plot(EVIresample_NAD83,
     col = terrain.colors(255),
     main = "EVI Raster")

mtext("Source: NASA; MODIS enhanced vegetation index (EVI)", side = 1, line = 4)

#turn off png
dev.off()
## png 
##   2

plot and save LST raster

plot(LSTresample_NAD83,
     col = terrain.colors(255),
     main = "LST Raster")

mtext("Source: NASA; MODIS Land Surface Temperature and Emissivity (LST)", side = 1, line = 4)
Fig 4. LST Raster.

Fig 4. LST Raster. 

png("LST_plot.png", width = 800, height = 600)

plot(LSTresample_NAD83,
     col = terrain.colors(255),
     main = "LST Raster")

mtext("Source: NASA; MODIS Land Surface Temperature and Emissivity (LST)", side = 1, line = 4)

dev.off()
## png 
##   2

creating the cov_stack

#create cov_stack of EVI and LST rasters
output_folder <- "rastersResampled"
output_EVI <- "output_EVI.tif"
output_EVI_filepath <- file.path(output_folder, output_EVI)
writeRaster(EVIresample_NAD83, filename = output_EVI_filepath, format = "GTiff", overwrite=TRUE)
output_LST <- "output_LST.tif"
output_LST_filepath <- file.path(output_folder, output_LST)
writeRaster(LSTresample_NAD83, filename = output_LST_filepath, format = "GTiff", overwrite=TRUE)
cov_stack <- getCovariateRasters('rastersResample', '.tif$', micounty2020)

plot(cov_stack,
     col = terrain.colors(255),
     main = c("EVI", "LST"))

mtext("Source: NASA; MODIS enhanced vegetation index (EVI) & MODIS Land Surface Temperature and Emissivity (LST)", side = 1, line = 4)
Fig 5. cov-stack of EVI and LST Raster.

Fig 5. cov-stack of EVI and LST Raster. 

#plot and save the cov_stack
png("cov_stack_plot.png", width = 800, height = 600)

plot(cov_stack,
     col = terrain.colors(255),
     main = c("EVI", "LST"))

mtext("Source: NASA; MODIS enhanced vegetation index (EVI) & MODIS Land Surface Temperature and Emissivity (LST)", side = 1, line = 4)

dev.off()
## png 
##   2

this code clips and saves the rasterPOP to the EVIresample_NAD83, it is excluded from running here because it takes ~20 min

#load in the cipped rasterPOP file
rasterPOP_clip <- raster("rastersResampled/output_POP.tif")
#set POP raster to same extent as other rasters
POPresample <- projectRaster(rasterPOP_clip, rasterEVI_max, method = 'ngb')
#set target crs and project each resampled raster to be in the same extent as the MI County data
POPresample_NAD83 <- projectRaster(POPresample, crs = target_crs)
#make POP raster a SpatRaster class
SpatRaster_POP <- as(POPresample_NAD83, "SpatRaster")
#get extent of cov_stack
rasterExt <- ext(cov_stack)
#set new SpatRaster_POP to same extent as cov_stack
SpatRaster_POP <- crop(SpatRaster_POP, rasterExt)

plot of SpatRaster_POP

#plot and save SpatRaster_POP
plot(SpatRaster_POP,
     col = terrain.colors(255),
     main = "Michigan Population")

mtext("Source: WorldPop, United States Population 2020", side = 1, line = 4)
Fig 6. WorldPop Raster.

Fig 6. WorldPop Raster. 

png("SpatRaster_POP.png", width = 800, height = 600)

plot(SpatRaster_POP,
     col = terrain.colors(255),
     main = "Michigan Population")

mtext("Source: WorldPop, United States Population 2020", side = 1, line = 4)

dev.off()
## png 
##   2

using the prepare_data function, with larger mesh edges

#prep data for the model
data_for_model_1 <- prepare_data(polygon_shapefile = micounty2020_clean,
                               covariate_rasters = cov_stack,
                               aggregation_raster = SpatRaster_POP,
                               response_var = 'Reported_h',
                               id_var = 'GEOID',
                               mesh.args = list(cut = 0.01,
                                                offset = c(0.1, 0.5),
                                                max.edge = c(0.2, 0.5),
                                                resolution = 250),
                               na.action = TRUE)

plotting the prepare_data function, with larger mesh edges

#plot and save data_for_model

plot(data_for_model_1)
Fig 7. Data for Model, with larger mesh edges.

Fig 7. Data for Model, with larger mesh edges. 

png("data_for_model_plot_1.png", width = 800, height = 600)

plot(data_for_model_1)

dev.off()
## png 
##   2

using the prepare_data function, with smaller mesh edges

#prep data for the model
data_for_model_2 <- prepare_data(polygon_shapefile = micounty2020_clean,
                               covariate_rasters = cov_stack,
                               aggregation_raster = SpatRaster_POP,
                               response_var = 'Reported_h',
                               id_var = 'GEOID',
                               mesh.args = list(cut = 0.01,
                                                offset = c(0.1, 0.5),
                                                max.edge = c(0.1, 0.2),
                                                resolution = 250),
                               na.action = TRUE)

plotting the prepare_data function, with smaller mesh edges

#plot and save data_for_model

plot(data_for_model_2)
Fig 8. Data for Model, with smaller mesh edges.

Fig 8. Data for Model, with smaller mesh edges. 

png("data_for_model_plot_2.png", width = 800, height = 600)

plot(data_for_model_2)

dev.off()
## png 
##   2

using the prepare_data function, with no mesh arguments defined

#prep data for the model
data_for_model_3 <- prepare_data(polygon_shapefile = micounty2020_clean,
                               covariate_rasters = cov_stack,
                               aggregation_raster = SpatRaster_POP,
                               response_var = 'Reported_h',
                               id_var = 'GEOID',
                               na.action = TRUE)

plotting the prepare_data function, with no mesh argumnets defined

#plot and save data_for_model

plot(data_for_model_3)
Fig 9. Data for Model, with no mesh arguments defined.

Fig 9. Data for Model, with no mesh arguments defined. 

png("data_for_model_plot_3.png", width = 800, height = 600)

plot(data_for_model_3)

dev.off()
## png 
##   2

Unfortunately, this is where my analysis was halted by the error. From the sources I can find online, “the “NA/NaN gradient evaluation” is a generic nlminb error message, which can occur whenever something goes wrong in the objective and/or gradient function computation such that a bad value (NA or NaN) is returned, at which point nlminb gives up.” I’ve tried changing the raster inputs and removing null values from all the raster inputs and polygon inputs but that did not resolve the issue.

gradient evaluation error gradient evaluation error code