Matching AfSIS spatial layers to OFRA legacy data

Extracting climate covariates

#Rgdal is what allows R to understand the structure of shapefiles by providing functions to read and 
#convert spatial data into easy to work with R dataframes. 
#Sp enables transformations and projections of the data and provides
#functions for working with the loaded spatial polygon objects.
#R to use .shp files.
# Some of these layers are in geographic coordinates;
#+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
#others in laea; +proj=laea +lat_0=5 +lon_0=20 +x_0=0 +y_0=0 +datum=WGS84
#+units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
#The two coordinate referencing systems have to be harmonised i.e all set as laea for example
#The purpose of the script is to load the the layer into R and then 
#extract the cell values within raster layers with reference to point data
#In this case for legacy field trial sites. The result is a polygonDataFrame
#with the information coming from the trial sites and combined with extracted 
#values, of which mean and standard deviation were calculated 
#this script forms part of script to construct environmental envelopes (regions of spatial similarity)
#based on trial locations
library(rgdal)
library(raster)
library(maptools)
#TRMM rainfall long term daily in mm LAEA 
#TRMM long term Average Precipitation Temporal Range 
TRMM_raster1=raster("C:/Users/machariam/Desktop/Spatial_layers/AfSIS_layers/TRMM3B42ALT.tif")
class(TRMM_raster1)
plot(TRMM_raster1)
TRMM_raster1
#TRMM rainfall long term daily SD LAEA
# Long term average number of rainy days temporal range
TRMM_raster2=raster("C:/Users/machariam/Desktop/Spatial_layers/AfSIS_layers/TRMM3B42ARDLT.tif")
class(TRMM_raster2)
plot(TRMM_raster2)
TRMM_raster2
#TRMM fourier LAEA
#TRMM long-term modified  fournier index temporal range
TRMM_raster3=raster("C:/Users/machariam/Desktop/Spatial_layers/AfSIS_layers/TRMM3B42MFILT.tif")
class(TRMM_raster3)
plot(TRMM_raster3)
TRMM_raster3
#Average annual precipitation in mm long term WorldCLIM longlat
#Long term average precipitation longlat
raster3=raster("C:/Users/machariam/Desktop/Spatial_layers/AfSIS_layers/BIO12ALT.tif")
class(raster3)
plot(raster3)
raster3
#Reprojecting raster3 from longlat to laea
newproj="+proj=laea +lat_0=5 +lon_0=20 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +TOWGS84=0,0,0"
raster_3 <- projectRaster(raster3, crs=newproj)
class(raster_3)
plot(raster_3)
raster_3
#Fourier rainfall WorldCLIM longlat
#BIOCLIM  12	LongTerm	Modified	Fournier	Index
rain_raster4=raster("C:/Users/machariam/Desktop/Spatial_layers/AfSIS_layers/BIO12MFILT.tif")
class(rain_raster4)
plot(rain_raster4)
rain_raster4
newproj="+proj=laea +lat_0=5 +lon_0=20 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +TOWGS84=0,0,0"
rain_raster_4 <- projectRaster(rain_raster4, crs=newproj)
rain_raster_4
#worldclim average temperature longlat
#BIOCLIM 1 longterm average temperature
raster0=raster("C:/Users/machariam/Desktop/Spatial_layers/AfSIS_layers/BIO1ALT.tif")
class(raster0)
plot(raster0)
raster0
newproj="+proj=laea +lat_0=5 +lon_0=20 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +TOWGS84=0,0,0"
raster_0 <- projectRaster(rain_raster4, crs=newproj)
raster_0

#modis average daily temperature laea
#LST day Long-Term average temporal range
raster1=raster("C:/Users/machariam/Desktop/Spatial_layers/AfSIS_layers/MY2LSTDALT200207_201412.tif")
class(raster1)
plot(raster1)
raster1
#modis average nightly temperature laea
#LST  Night	Long‐Term	Average	Temporal	Range	
raster2=raster("C:/Users/machariam/Desktop/Spatial_layers/AfSIS_layers/MY2LSTNALT200207_201412.tif")
class(raster2)
plot(raster2)
raster2
#Importing legacy data and reprojecting so as to merge the dataset with spatial covariates
#Importing legacy data and reprojecting
Kenya_crop_data=read.csv("C:/Users/machariam/Desktop/Spatial_layers/Kenya/kenya_legacy.csv",header=T)
coordinates(Kenya_crop_data)= ~X+Y
proj4string(Kenya_crop_data)=CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")
Kenya_legacy=spTransform(Kenya_crop_data,CRS("+proj=laea +lat_0=5 +lon_0=20 +x_0=0 +y_0=0 +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"), IDvar="id")
plot(Kenya_legacy)
class(Kenya_legacy)
Kenya_legacy
#Extracting covariates for the kenya legacy data
BIO12ALT=extract(raster_3,Kenya_legacy, small=TRUE, fun=mean)
BIO12ALT_sd=extract(raster_3,Kenya_legacy, small=TRUE, fun=sd)
BIO12MFI=extract(rain_raster_4,Kenya_legacy, small=TRUE, fun=mean)
BIO12MFI_sd=extract(rain_raster_4,Kenya_legacy, small=TRUE, fun=sd)

BI01ALT=extract(raster_0,Kenya_legacy, small=TRUE, fun=mean)
BI01ALT_sd=extract(raster_0,Kenya_legacy, small=TRUE, fun=sd)

MY2LSTDALT=extract(raster1,Kenya_legacy, small=TRUE, fun=mean)
MY2LSTDALT_sd=extract(raster1,Kenya_legacy, small=TRUE, fun=sd)

MY2LSTNALT=extract(raster2,Kenya_legacy, small=TRUE, fun=mean)
MY2LSTNALT_sd=extract(raster2,Kenya_legacy, small=TRUE, fun=sd)

TRMM3B42ALT=extract(TRMM_raster1,Kenya_legacy, small=TRUE, fun=mean)
TRMM3B42ARDLT=extract(TRMM_raster2,Kenya_legacy, small=TRUE, fun=mean)
TRMM3B42MFILT=extract(TRMM_raster3,Kenya_legacy, small=TRUE, fun=mean)

#make data frames for extracted values
Kenya_legacy_1=as.data.frame(Kenya_legacy)
Kenya_legacy_1
rain_CLIM=cbind(Kenya_legacy_1,BIO12ALT,BIO12ALT_sd, BIO12MFI,BIO12MFI_sd,
                BI01ALT,BI01ALT_sd,MY2LSTDALT, MY2LSTDALT_sd,MY2LSTNALT,MY2LSTNALT_sd,
                TRMM3B42ALT,TRMM3B42ARDLT,TRMM3B42MFILT)
write.csv(rain_CLIM, "/Users/machariam/Desktop/Spatial_layers/Kenya/rain_CLIM.csv")