Loading raster, reprojection, extract covariates

OFRA Spatial

#EDA OFRA Spatial
#load neccessary packages
library(maptools)
library(sp)
library(rgdal)


#Importing data;
#maptools: Tools for Reading and Handling Spatial Objects
#rgdal has issues with IDs when importing spatial objects

#loading the crop mask and setting projection
#Getting crop mask shape file into R and set projection
kenya_crop_mask=readShapeSpatial("C:/Users/machariam/Desktop/Spatial_layers/Kenya/Agric_Land/Kenya_Agric_land.shp", proj4string=CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"))
plot(kenya_crop_mask)
names(kenya_crop_mask)
class(kenya_crop_mask)
#Reproject vector file to LAEA
kenya_crop_cover=spTransform(kenya_crop_mask, CRS("+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"))

#loading the kenya map
kenya_county=readShapeSpatial("C:/Users/machariam/Desktop/Spatial_layers/Kenya/County/County.shp", proj4string=CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"))
kenya_admin_map=spTransform(kenya_county, CRS("+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"))
plot(kenya_admin_map)
class(kenya_admin_map)

#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"))
plot(Kenya_legacy)
class(Kenya_legacy)

#Random selection of points on the polygon
kcp=spsample (kenya_crop_cover, n=200, type='random')
plot(kcp)
class(kcp)

#Overlying the 3 spatial layers
plot(kenya_admin_map)# County map
plot(kenya_crop_cover, col='red', add=TRUE)# crop distribution layer
plot(kcp, col='orange',pch=3,cex=0.8, add=TRUE)#Random sample from crop distribution layer
plot(Kenya_legacy, col='blue',cex=0.8,add=TRUE)#legacy data results

#Testing the plots separately
par(mfrow=c(2,2))
plot(Kenya_legacy, col='blue',pch=3,cex=0.5)#legacy data results
plot(kenya_admin_map, col="light green")
plot(kenya_crop_cover, col='red')# crop distribution layer
plot(kcp, col='orange',pch=3,cex=0.5)

#overlay kenya admin on crop trials
plot(kenya_admin_map, col="light green")
plot(Kenya_crop_data, col='blue',pch=3,cex=0.5, add=TRUE)

#Load spatial layer

library(raster)
library(rgdal)
library(RCurl)
url = "ftp://disc2.nascom.nasa.gov/data/TRMM/Gridded/3B43_V7/2012/"
filenames = getURL(url, ftp.use.epsv = FALSE, ftplistonly=TRUE, crlf=TRUE)
filenames = paste(url, strsplit(filenames, "\r*\n") [[1]], sep="")
filenames = filenames[grep(filenames, pattern="*precipitation.accum")] 
# this line allows to download only files with a certain pattern, 
#e.g. only certain months or days. "*precipitation.accum" 
#means monthly accumulated rainfall here.
filenames # list all files
mapply(download.file, filenames, basename("ftp://disc2.nascom.nasa.gov/data/TRMM/Gridded/3B43_V7/2012")) # download files

##################################################################################################
###########EDA continuation#######################################################################
#load neccessary packages
library(maptools)
library(sp)
library(rgdal)


#Importing data;
#maptools: Tools for Reading and Handling Spatial Objects
#rgdal has issues with IDs when importing spatial objects

#loading the crop mask and setting projection
#Getting crop mask shape file into R and set projection
kenya_crop_mask=readShapeSpatial("C:/Users/machariam/Desktop/Spatial_layers/Kenya/Agric_Land/Kenya_Agric_land.shp", proj4string=CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"))
plot(kenya_crop_mask)
names(kenya_crop_mask)
class(kenya_crop_mask)
#Reproject vector file to LAEA
kenya_crop_cover=spTransform(kenya_crop_mask, CRS("+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"))

#loading the kenya map
kenya_county=readShapeSpatial("C:/Users/machariam/Desktop/Spatial_layers/Kenya/County/County.shp", proj4string=CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"))
kenya_admin_map=spTransform(kenya_county, CRS("+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"))
plot(kenya_admin_map)
class(kenya_admin_map)

#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"))
plot(Kenya_legacy)
class(Kenya_legacy)

#Random selection of points on the polygon
kcp=spsample (kenya_crop_cover, n=200, type='random')
plot(kcp)
class(kcp)

#Overlying the 3 spatial layers
plot(kenya_admin_map)# County map
plot(kenya_crop_cover, col='red', add=TRUE)# crop distribution layer
plot(kcp, col='orange',pch=3,cex=0.8, add=TRUE)#Random sample from crop distribution layer
plot(Kenya_legacy, col='blue',cex=0.8,add=TRUE)#legacy data results

#Testing the plots separately
par(mfrow=c(2,2))
plot(Kenya_legacy, col='blue',pch=3,cex=0.5)#legacy data results
plot(kenya_admin_map, col="light green")
plot(kenya_crop_cover, col='red')# crop distribution layer
plot(kcp, col='orange',pch=3,cex=0.5)

#overlay kenya admin on crop trials
plot(kenya_admin_map, col="light green")
plot(Kenya_crop_data, col='blue',pch=3,cex=0.5, add=TRUE)
########################################################################
##############Extract climate covariates#################################
# 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 referenc 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)
#The rgdal package reads numerous spatial data formats using the standalone, 
#open-source Geospatial Data Abstraction Library (GDAL). This is also the recommended package for 
#performing map projections and datum transformations, using the external PROJ.4 library.
library(raster)
library(maptools)
#The maptools package provides functions for reading, writing, and manipulating spatial data, 
#especially shapefiles.
#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)
Kenya_crop_data
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")

#https://help.nceas.ucsb.edu/r:spatial