Regression kriging

Regression kriging

#Regression kriging for maize 50 kg N

library(rgdal)
library (raster)
library (maptools)
library(rgeos)
library(dismo)
library(sp)

#load tiffs of standardized variables and PCA
glist <- list.files(path="/Users/alexverlinden/grassdata/AF_grids/AF_geotiff_std", pattern='tif', full.names=T) #loads standaridzed tiffs
glist2 <- list.files(path="/Users/alexverlinden/grassdata/AF_PC_OFRA", pattern='tif', full.names=T) #loads 10 PCs

#stack grids
grids <- stack(glist)
grids2=stack(glist2)

#read legacy data
maize_50=readShapeSpatial("/Users/alexverlinden/grassdata/AfricaLaea/legacy_maize/maize_50N.shp")
proj4string(maize_50)=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(maize_50, axes = TRUE)
#dataframe
m=as.data.frame(maize_50)
x=m[,1:2]
points=extract(grids2,x)
pts=extract(grids,x)
reg=cbind(m,points)
regr=cbind(m,pts)
#glm on PCs
maize.lm=lm(reg$Y50_Y0~reg$AfGYGPC1+reg$AfGYGPC2+reg$AfGYGPC3+reg$AfGYGPC4+reg$AfGYGPC5+reg$AfGYGPC6+reg$AfGYGPC7+reg$AfGYGPC8+reg$AfGYGPC9+reg$AfGYGPC10)
maize.lm2=lm(regr$Y50_Y0~regr$BSAN+regr$BSAS+regr$BSAV+regr$CTI+regr$ELEV+regr$EVI+regr$LSTD+regr$LSTN+regr$REF1+regr$REF2+regr$REF3+regr$REF7+regr$RELI+regr$TMAP+regr$TMFI)