Tsetse-fly-Ensemble-Modelling-Kenya-data_processing.R

# ------------------------------------------------------------------------------
# Script: Tsetse Habitat Suitability Modelling - Kenya
# Purpose: Processes raster datasets for modelling purposes
# Context: Supports tsetse habitat suitability modelling in Kenya
# Author: Raphael Mongare
# R Version: 4.3.2
# ------------------------------------------------------------------------------

# ----------------------------- START ------------------------------------------
## Install needed packages
install.packages('rasterVis')
install.packages('gdalUtils', dependencies = TRUE)
install.packages('rgdal')
install.packages('ggplot2')
install.packages('terra')
install.packages('sp')
install.packages('sf')
install.packages('stars')
install.packages("devtools")
devtools:::install_github("gearslaboratory/gdalUtils")

#-------------------------------------------------------------------------------
## load needed librairies
library(raster)
library(rasterVis)
library(gdalUtils)
library(rgdal)
library(ggplot2)
library(sp)
library(terra)
library(stars)
library(sf)

#-------------------------------------------------------------------------------
## set working directory
setwd('C:/Raw_Data_Annual')
list.files()

## Define the projection
sr = "+proj=utm +zone=37 +south +datum=WGS84 +units=m +no_defs"

## Read the shp tsetse traps####
tsetse_traps<- readOGR('Tsetse-occuerence.shp') ## convert csv file into a shp
tsetse_Prj <- spTransform(tsetse_traps, CRS(sr))

## Read the kenya shapefile (download Kenya shapefiles from gadm: https://gadm.org/maps/KEN.html)
ke <- readOGR('Kenya.shp')

## Project Kenya shapefile

kenya = spTransform(ke, CRS(sr))
plot(kenya)

plot(tsetse_Prj, add = TRUE)

## Create a mask layer to use in the projection and resampling the predictors
masklayer = raster("mask.tif")
ML <- projectRaster(masklayer, crs = crs(sr))
plot(ML)

#-------------------------------------------------------------------------------
## Scale the variables that require scaling
## working with the stacked version is easier. Read all the single rasters and then create a stack
skintemp_2017_2020 <- stack('skintemp_combat.tif') - 273.15
soiltemp_2017_2020 <- stack('soiltemp_combat.tif') - 273.15
airtemp_2017_2020 <- stack('airtemp_combat.tif') - 273.15
lst_day_2017_2020 <- stack('lst_day_combat.tif') * 0.02 - 273.15
lst_night_2017_2020 <- stack('lst_night_combat.tif') * 0.02 - 273.15
lst_med_2017_2020 <- stack('lst_med_combat.tif') * 0.02 - 273.15
totprec_2017_2020 <- stack('totprec_combat.tif') * 1000
ndvi_2017_2020 <- stack('ndvi_combat.tif')
sbsmoist_2017_2020 <- stack('sbsmoist_combat.tif')
ssmoist_2017_2020 <- stack('ssmoist_combat.tif')
surfRunoff_2017_2020 <- stack('surfRunoff_combat.tif')
hist(lst_night_2017_2020)

#-------------------------------------------------------------------------------
## Stack all of them together
env_17_20 <- stack(skintemp_2017_2020,soiltemp_2017_2020,airtemp_2017_2020,
                             lst_day_2017_2020,lst_night_2017_2020,lst_med_2017_2020,
                             totprec_2017_2020,ndvi_2017_2020,sbsmoist_2017_2020,
                             ssmoist_2017_2020,surfRunoff_2017_2020)

env_17_20_prj <- projectRaster(env_17_20, crs = crs(sr))

#-------------------------------------------------------------------------------

## process other variable with different extent from the rest of the variables
## mask and crop the files

## Elevation
elev <- raster('Elevation.tif')
elev_prj = projectRaster(elev, crs = crs(sr))
masked <- mask(elev_prj, kenya)
elev_crop <- crop(masked, extent(kenya))

## Sand content
sand <- raster('sand_combat.tif')
sand_prj <- projectRaster(sand, crs = crs(sr))
masked <- mask(sand_prj, kenya)
sand_crop <- crop(masked, extent(kenya))

## Topographic wetness index (TWI: https://doi.org/10.1016/j.scitotenv.2020.143785)
twi <- raster('NewTWI.tif')
twi_prj <- projectRaster(twi, crs = crs(sr))
masked <- mask(twi_prj, kenya)
twi_crop <- crop(masked, extent(kenya))

## Human population density
pop <- raster('Kenya_Human_Pop_2020_1km.tif')
pop_prj <- projectRaster(pop, crs = crs(sr))
masked <- mask(pop_prj, kenya)
pop_crop <- crop(masked, extent(kenya))

## Cattle densisty
cattle <- raster('Cattle_Density_2020.tif')
cattle_prj <- projectRaster(cattle, crs = crs(sr))
masked <- mask(cattle_prj, kenya)
cattle_crop <- crop(masked, extent(kenya))

## Distance to protected areas (PA) (Proximity to PA)
protected_areas <- raster('EucDist_ProtectedAreas.tif')
protected_areas_prj <- projectRaster(protected_areas, crs = crs(sr))
masked <- mask(protected_areas_prj, kenya)
protected_areas_crop <- crop(masked, extent(kenya))

## Buffalo density
buffalo <- raster('Buffalo.tif')
buffalo_prj <- projectRaster(buffalo, crs = crs(sr))
masked <- mask(buffalo_prj, kenya)
buffalo_crop <- crop(masked, extent(kenya))

## Sheep density
sheep <- raster('Sheep.tif')
sheep_prj <- projectRaster(sheep, crs = crs(sr))
masked <- mask(sheep_prj, kenya)
sheep_crop <- crop(masked, extent(kenya))

## Goat density
goat <- raster('Goat.tif')
goat_prj <- projectRaster(goat, crs = crs(sr))
masked <- mask(goat_prj, kenya)
goat_crop <- crop(masked, extent(kenya))

#-------------------------------------------------------------------------------
### Crop to Kenya counties only
##First mask the files then apply the crop function
masked <- mask(env_17_20_prj, kenya)
env_17_20_crop <- crop(masked, extent(kenya)) ##Crop the masked files
plot(env_17_20_crop)

#-------------------------------------------------------------------------------
## Resample the cropped files and create a dummy raster for reference
extent(protected_areas_crop)## choose one raster file that you want to use as the 
                            ## base for resampling all the layers to
ncol(protected_areas_crop)
nrow(protected_areas_crop)
protected_areas_crop

e = extent(-67084.8, 822915.2, 9479734, 10566734 ) # use the coordinates of the base layer to set the extent
s = raster(e, nrows=1087, ncols=890, crs=protected_areas_crop@crs) #generate a raster with extent and number of columns of the base raster

#-------------------------------------------------------------------------------
#Resample the variables
env_17_20.resample<-resample(env_17_20_crop,s,method="ngb")
plot(env_17_20.resample)

elev.resample<-resample(elev_crop,s,method="ngb")
sand.resample<-resample(sand_crop,s,method="ngb")
twi.resample<-resample(twi_crop,s,method="ngb")
human_pop.resample<-resample(pop_crop,s,method="ngb")
cattle.resample<-resample(cattle_crop,s,method="ngb")
protected_areas.resample<-resample(protected_areas_crop,s,method="ngb")
buffalo.resample<-resample(buffalo_crop,s,method="ngb")
sheep.resample<-resample(sheep_crop,s,method="ngb")
goat.resample<-resample(goat_crop,s,method="ngb")

#-------------------------------------------------------------------------------
##Stack all variables
annual_var_stack <- stack(env_17_20.resample,elev.resample,sand.resample,
                      twi.resample,human_pop.resample,cattle.resample,
                      protected_areas.resample,buffalo.resample,sheep.resample,
                      goat.resample)
names(annual_var_stack)

#-------------------------------------------------------------------------------
## Extract raster values to points
annual_rasValue <- extract(annual_var_stack, tsetse_Prj)

### combine the points with the raster values and save as a .csv
annual_rasValues <- cbind(tsetse_Prj,annual_rasValue)

## Save the extracted csv file for the subsequent steps

write.table(annual_rasValues,file='Annual-Tsetse-predictors.csv', append=FALSE, sep= ",", row.names = FALSE, col.names=TRUE)

## Save the processed rasters for the subsequent steps
writeRaster(annual_var_stack, names(annual_var_stack), bylayer=TRUE, format='GTiff',overwrite=TRUE)

#-------------------------------------------------------------------------------
##      Data cleaning

## Read the csv with all the predictors from the directory
pred.data <- read.csv("Annual-Tsetse-predictors.csv")
View(pred.data)


## Remove the rows with NA from the dataset
tseClean <- na.omit(pred.data)
View(tseClean)
summary(tseClean)
print(tseClean)

## Remove the Buffalo column using dplyr package

library("dplyr")

print(tseClean[ , !names(tseClean) %in% 
      c('Buffalo')])
summary(tseClean)
str(tseClean)

#-------------------------------------------------------------------------------
## Save the cleaned data for further processing
write.table(tseClean,file='Tsetse-predictors-Cleaned.csv', append=FALSE, sep= ",", row.names = FALSE, col.names=TRUE)

## Plot two columns from the dataframe to check for completeness
dat <- select(tseClean,elevation, maxNDVI17)
View(dat)
plot.default(dat)

#-------------------------------------- END -----------------------------------------------------