Forest diversity with GAM

#! Modelying forest diversity using GAM
#! Paulo Eduardo Cardoso
require(RColorBrewer)
require(graphics)
require(mgcv)
library(maptools)
.trPaths <- paste(paste(Sys.getenv('APPDATA'), '\\Tinn-R\\tmp\\', sep=''), 
  c('', 'search.txt', 'objects.txt', 'file.r', 'selection.r', 'block.r', 
'lines.r'), sep='') 

limite_pt <- readShapeSpatial("D:/Sig/Administrativo/Portugal/LimitePTSimplificadoIgeoE2010.shp") #! read shapefile with 2 polys: country limits and buffer
setwd("D:/Strix/Flora_Digital/")
a<-read.table("Riqueza_Especies_Forestais.csv", sep=",", header=T)
names(a)[5:6]<-c("x", "y")

#! Construir o modelo GAM para os dados do Nº de espécies em função das coordenadas x e y
g<-with(a, gam(NumSp~s(x,y, k=150), gamma=1))
#plot(g, se=F, asp=1)

#! == Criar uma malha fina para estimar os valores ajustados com o GAM
#! Método expand.grid
resolucao<-1000 #! 2500 pelo Zé Pedro
x1<-seq(min(a$x),max(a$x), by=resolucao)
y1<-seq(min(a$y), max(a$y), by=resolucao)
valores<-expand.grid(x=x1, y=y1)
valores<-as.data.frame(unclass(valores))
ovalores<-valores
coordinates(ovalores)<- ~x+y #! SPointDF
valores$pos <- overlay(ovalores,limite_pt) #! Identify in/out points
#valores <- valores[valores$pos==2,][1:2] #! Subset


#! Método espacial - pontos aleatórios
so<-sample.Spatial(limite_pt,100000,"random") #! Random points inside a polygon
opt<-overlay(so,limite_pt)
pto_pt<- so[is.na(opt)]
pto_pt <- as.data.frame(pto_pt)
pto_pt[1:10,]
names(pto_pt)<-c("x","y")
x
plot(pto_pt,asp=1)

#! == Estimativas com base no GAM para os pontos da malha fina

pred<-predict(g, valores[,1:2], type="response", se=F)

valores$classes<-unclass(cut(pred, breaks=seq(min(pred)-.0001, max(pred), length=120)))
valores$pred <- pred
valores$pred[valores$pos==1] <- NA
#! == Gráficos
#par(cex.lab=0.6)
#with(a, plot(x,y, col=rainbow(131)[NumSp+1], pch=19, asp=1,cex=.7))
plot(limite_pt[2,], col="NA", border="NA",asp=1) # set up background polygon and paint hole
with(valores[valores$pos==2,][1:2], plot(x,y,xlim=c(70000,370000),ylim=c(0,600000),col=rainbow(120)[valores$classes[valores$pos==2]], asp=1, pch=19,cex.axis=0.7),add=T)
with(valores,contour(x1, y1, ldw=2,labcex=0.5,nlevels=8,matrix(unclass(valores$pred), nrow=length(x1), ncol=length(y1)), add=T,col="grey10"))  #   length(x1)
plot(limite_pt[2,], col="NA", border="black",asp=1,add=T) # set up background polygon and paint hole

#*lengthwin.graph()
filled.contour(x1, y1, matrix(unclass(pred), nrow=length(x1), ncol=length(y1)),asp=1)
#plot(valores, col=rainbow(131)[classes], asp=1, pch=19, cex=.7)
#contour(x1, y1, matrix(unclass(pred), nrow=200, ncol=200), add=T)
#win.graph()

#! With Image
plot(limite_pt[2,], col="NA", border="NA",asp=1) # set up background polygon and paint hole
with(valores,image(x1, y1, asp=1,matrix(unclass(valores$pred), nrow=length(x1), ncol=length(y1)),col= brewer.pal(9,"RdYlGn"),add=T))
with(valores,contour(x1, y1, ldw=2,labcex=0.5,nlevels=8,matrix(unclass(valores$pred), nrow=length(x1), ncol=length(y1)), add=T,col="grey10"))
plot(limite_pt[2,], col="NA", border="black",add=T,asp=1) # set up background polygon and paint hole

with(a, points(x,y, col=rainbow(131)[NumSp+1], pch=19, asp=1,cex=.7))