# `r spname`


```{r,echo=FALSE,warning=F,message=FALSE,results='hide'}

library(rgdal)
library(climatol)
library(sm)
library(vegan)
library(dismo)
library(rgeos)
library(mgcv)
library(randomForest)
fldir<-"/var/www/apps/Fichas/resources"
mexico<-readOGR("/var/www/apps/Fichas/resources/shapefiles","mexico")
meso<-readOGR("/var/www/apps/Fichas/resources/shapefiles","meso")
preds<-brick("/var/www/apps/Fichas/resources/predictors.tif")
ndvi<-brick("/var/www/apps/Fichas/resources/ndvimeso.tif")
soil<-brick("/var/www/apps/Fichas/resources/Soil.tif")
clim<-brick("/var/www/apps/Fichas/resources/Worldclim2min.tif")
back<-readGDAL("/var/www/apps/Fichas/resources/backdrop.png")
prednames<-c("X","Y","elev","anprec","mtemp","Trange","saet","msoil")
names(preds)<-prednames
preds1<-stack(preds[[c(5,6)]],soil[[c(1,3,5,7,9)]])
preds1<-brick(preds1)

 clean<-function(x){
    x<-as.character(x)
    x<-gsub("[^0-9]","",x)
    as.numeric(x)}
  



```



```{r,echo=FALSE,warning=F,message=FALSE}
spname<-"Quercus_segoviensis"
library(RODBC)
con<-odbcConnect("reddeam")
pts<-sqlQuery(con,"select st_astext(pts) y from models.spatial where binomial like 'Quercus_seg%'")
pts<-readWKT(pts$y)


```

```{r,echo=FALSE}
figpath <-paste("/var/www/apps/Fichas/results/figure/",spname,"/",sep="")
 opts_chunk$set(fig.path = figpath,echo=FALSE,warnings=FALSE)

```


```{remove_last r,warning=FALSE,echo=FALSE}

## Remove variables from last run so that if the model fails we simply get no output.
rm(vars,vars1,nd,sl,e,mod,gampreds,gampreds1,ndvimod,moddata,moddata1)
pts@proj4string<-CRS("+init=epsg:4326")
```

```{r range_map,echo=FALSE,warning=F,message=FALSE}


## Crop the predictor rasters to a 3 degree buffer of the convex hull
pts_proj<-spTransform(pts, CRS("+init=epsg:3408"))

chull<-convHull(pts_proj)@polygons


chull_area<-gArea(chull)/10^9
chull_ex<-gBuffer(chull,width=150000) #150 km buffer
chull<-spTransform(chull, CRS("+init=epsg:4326"))
chull_ex<-spTransform(chull_ex, CRS("+init=epsg:4326"))
preds<-crop(preds,chull_ex)
preds1<-crop(preds1,chull_ex)

cvhull<-gUnaryUnion(gBuffer(pts_proj,width=20000))
cvhull<-gBuffer(cvhull,width=-10000)
cvhull_area<-gArea(cvhull)/10^9
cvhull<-spTransform(cvhull, CRS("+init=epsg:4326"))


plot(chull_ex,main="")
image(back,r=1,g=2,b=3,add=T)
#contour(preds,3,add=T,nlevels=7,lwd=0.4)
image(back,r=1,g=2,b=3,add=T)
tit<-paste(spname,"\n Convex hull area = ",round(chull_area,1)," thousand km2 \n Concave buffer area = ",round(cvhull_area,1)," thousand km2 ",sep="")
title(main=tit)
plot(chull_ex,add=T,border=2,lwd=3)
plot(chull,add=T,border=2,lwd=1)
plot(cvhull,border="darkblue",lwd=2,add=T)

plot(pts,add=T,pch=21,bg=2,cex=0.3)
plot(meso,add=T)

```



```{r extract_vars,echo=FALSE,warning=F,message=FALSE}

# Extract the background points from both sets of predictors.

## The object pts has the presence coordinates. So after this vars constains the variables from the first set of predictors and vars1 from the second

backgrpts <- randomPoints(preds, 1000)
backgrpts<-data.frame(backgrpts)
coordinates(backgrpts)<-~x+y
backgr<-na.omit(data.frame(coordinates(backgrpts),extract(preds,backgrpts)))
vars<-na.omit(data.frame(coordinates(pts),extract(preds,pts)))

```



## Seasonal patterns

A standard Water and Leith climate diagram for the mean values of precipitation and temperature extracted from the species presence points.


```{r,echo=FALSE,warning=F,message=FALSE}
climate<-as.data.frame(extract(clim,pts))
climate<-apply(climate,2,mean)
diag<-rbind(climate[25:36],climate[1:12]/10,climate[13:24]/10,(climate[13:24]/10)-6)
diagwl(diag)

```

The Walter and Leith diagram assumes that the growing season occurs when rainfall is over 100mm. A more refined method is to extract the values from a bucket model that keeps track of input to the soil profile through precipitation and reductions in soil moisture through evaptranspiration over the course of the year. This can be compared to changes in NDVI at the collection points.


```{r seasonal,echo=FALSE,warning=F,message=FALSE}

nd<-as.data.frame(extract(ndvi,pts))
nd<-(nd-1000)/9000
nd<-data.frame(stack(nd))
nd$ind<-clean(nd$ind)
minndvi<-min(tapply(nd$values,nd$ind,mean))
maxndvi<-max(tapply(nd$values,nd$ind,mean))
ndvi_index<-(maxndvi-minndvi)/maxndvi
tit<-paste("NDVI \n Seasonality = ",round(ndvi_index*100,0),"%",sep="")
boxplot(nd$values~nd$ind,main=tit,ylim=c(0.2,1))

ndvimod<-gam(values~s(ind),data=nd)
ind<-1:12
lines(1:12,predict(ndvimod,newdata=data.frame(ind)),lwd=2,col="darkgreen")
sl<-as.data.frame(extract(soil,pts))

sl<-(sl-200)/300
sl<-data.frame(stack(sl))
sl$ind<-clean(sl$ind)
minsoil<-min(tapply(sl$values,sl$ind,function(x)mean(x,na.rm=T)))
maxsoil<-max(tapply(sl$values,sl$ind,function(x)mean(x,na.rm=T)))
soil_index<-(maxsoil-minsoil)/maxsoil


tit<-paste("Soil moisture \n Seasonality = ",round(soil_index,2),"%",sep="")
boxplot(sl$values~sl$ind,main=tit,ylim=c(0,1))

soilmod<-gam(values~s(ind),data=sl)
ind<-1:12
lines(1:12,predict(soilmod,newdata=data.frame(ind)),lwd=2,col=2)


```




## Comparison

In some cases NDVI will remain fairly constant, even when the balance model shows that soil water constant is lowered for part of the year. Providing SWC is above 50% of maximum levels the vegetation would not experience a great deal of hydric stress.

```{r compare,echo=FALSE,warning=F,message=FALSE}
plot(1:12,predict(soilmod,newdata=data.frame(ind)),lwd=2,col=2,type="l",ylim=c(0,1),ylab="Seasonality",main="Seasonality")
lines(1:12,predict(ndvimod,newdata=data.frame(ind)),lwd=2,col="darkgreen")
grid()

```

## Summary of the variable values at the collection points.


```{r summary,results='asis',echo=FALSE}
library(xtable)
print(xtable(t(apply(vars,2,summary))),"html")
```

## Niche space with relation to annual precipitation and mean annual temperature.

The following diagrams show kernel densities one two synthestic climate axes (total annual rainfall and mean annual temperature). If their are signs of multimodality this may indicate that the species has not fully explored its climate niche, or that there are disjunct populations with differing characteristics. The method will not show clear results for species with few collection points.

```{r kernel,echo=FALSE,warning=F,message=FALSE}
sm.density(x=vars[,c(5,4)],nbins=0,display="image",xlab="Annual temperature",ylab="Annual rainfall (mm)",hull=T)
sm.density(vars[,c(5,4)], display = "persp", col = "grey")

```


## Spatial clustering

The same analysis can be run to look at spatial clustering. The kernel densities are smoothed, so will only suggest multimodality if the points are very highly clustered.


```{r spatial_kernel,echo=FALSE,warning=F,message=FALSE}
a<-sm.density(vars[,c(1,2)],nbins = 0,display="image",xlab="Long",ylab="Lat",col.palette=topo.colors(12),ngrid=100)

plot(meso,add=T)
points(pts,pch=21,bg=2,cex=0.2)
```


## Finding spatial clusters

The significance of any spatial clusters can be checked using the silhouette width method. The width is calculated for values of k between 2 and 5. If any are higher than 0.52 the analysis will produce a diagram showing the clusters. 

```{r spatial_clusters,echo=FALSE,warning=F,message=FALSE}
txt<-""
## In order to use Euclidean distance for all variables need to mormalise the data (make the range from 0 to 1)
normalise<-function(x)(x-min(x))/(max(x)-min(x))
varsn<-data.frame(apply(vars,2,normalise))

library(fpc)
k<-2:5
f<-function(k)pam(vars[,c(1,2)], k) $ silinfo $ avg.width
kscores<-unlist(lapply(k,f))
plot(k,kscores,type="b")
if (max(kscores)>0.52){
   clustering<-TRUE
   bestk<-which.max(kscores)+1
   #plot(pam(vars[,c(1,2)], bestk))
   fit <- kmeans(vars[,c(1,2)],bestk)
   plot(pts, col = fit$cluster)
   plot(meso,add=T)
   axis(1)
   axis(2)
   box()
   txt<-paste("The collections form  ",bestk," spatial clusters",sep="")
   ch<-convHull(pts,n=bestk)
plot(ch@polygons,add=T)

}


if (max(kscores)<=0.52){
  clustering<-FALSE
  txt<-"No evidence of spatial clustering"
  }


txt
```

`r txt`

## Anosim

If there is evidence that the points fall into at least 2 groups, but fewer than 6 we can look at whether there is significant differences in variability between and within groups in the climatic conditions at the sites using Anosim. This is a sensitive test, as would be MANOVA, so there will often be significant differences. They should only be intepreted as important if R is much larger than 0.3.



```{r anosim,echo=FALSE,warning=F,message=FALSE}
if (clustering==TRUE) {

dis <- dist(varsn[,-c(1,2)],'euclidean')
a<-anosim(dis,fit$cluster,perm=100)
summary(a)
plot(a)
fit <- lm(vars$mtemp~ as.factor(fit$cluster))
anova(fit)
}

```

`r txt`



## Gam model using simple environmental variables

For comparison we can fit a model using mean temperature, temperature range and annual precipitation.

```{r,echo=FALSE,warning=F,message=FALSE}

moddata<-data.frame(cbind(pb = c(rep(1, nrow(vars)), rep(0, nrow(backgr))),rbind(vars, backgr)))
cutpoint<-median(vars$x)
East<-moddata[moddata$x<cutpoint,]
West<-moddata[moddata$x>=cutpoint,]
```


```{r,echo=FALSE,warning=F,message=FALSE}
GAM_form<-as.formula(paste("pb ~ s(", paste(names(preds)[c(4,5,6)], collapse=",k=4) + s("), ",k=4)", collapse=""))
```


```{r,echo=FALSE,warning=F,message=FALSE}

gammod<-gam(GAM_form,data=moddata)
par(mfcol=c(2,2))
plot(gammod)
eback<-evaluate(vars,backgr,gammod)
lim<-0.9
threshold<-eback@t[which.max(eback@TPR<lim)]

```

```{r echo=FALSE,warning=FALSE}
par(mfcol=c(2,2))
modpreds<-predict(preds,gammod)
plot(modpreds)
plot(pts,add=T,pch=21,cex=0.2,bg=2)
plot(meso,add=T)
axis(1)
axis(2)
box()
modpres<-modpreds>threshold
plot(modpres)
plot(pts,add=T,pch=21,cex=0.4,bg=2)
plot(meso,add=T)
axis(1)
axis(2)
box()
eback<-evaluate(vars,backgr,gammod)
plot(eback,'ROC')

```


```{r,echo=FALSE,warning=FALSE,message=FALSE}

par(mfcol=c(2,2))
plot(modpres)
points(East$X,East$Y,pch=21,bg=East$pb+1,cex=0.3,col=East$pb+1)
abline(v=cutpoint)
box()

testmod<-gam(GAM_form,data=East)
e<-evaluate(West[West$pb==1,-1],West[West$pb==0,-1],testmod)
plot(e,'ROC')
plot(modpres)
points(West$X,West$Y,pch=21,bg=West$pb+1,cex=0.3,col=West$pb+1)
abline(v=cutpoint)
box()
testmod<-gam(GAM_form,data=West)
e<-evaluate(East[East$pb==1,-1],East[East$pb==0,-1],testmod)
plot(e,'ROC')



```



## Bucket model GAM.

Now fit a gam using temperature range, mean temperature and the annual soil moisture dynamic as input.

```{r gam_bucket,echo=FALSE,warning=FALSE}
backgr<-na.omit(data.frame(coordinates(backgrpts),extract(preds1,backgrpts)))
vars<-na.omit(data.frame(coordinates(pts),extract(preds1,pts)))
var_names<-names(vars)[-c(1:2)]
moddata<-data.frame(cbind(pb = c(rep(1, nrow(vars)), rep(0, nrow(backgr))),rbind(vars, backgr)))
cutpoint<-median(vars$x)
East<-moddata[moddata$x<cutpoint,]
West<-moddata[moddata$x>=cutpoint,]
GAM_form<-as.formula(paste("pb ~ s(", paste(var_names, collapse=",k=4) + s("), ",k=4)", collapse=""))

gammod<-gam(GAM_form,data=moddata)
eback<-evaluate(vars,backgr,gammod)
lim<-0.9
threshold<-eback@t[which.max(eback@TPR<lim)]


```




```{r echo=FALSE,warning=FALSE}
par(mfcol=c(2,2))
modpreds<-predict(preds1,gammod)
plot(modpreds)
plot(pts,add=T,pch=21,cex=0.2,bg=2)
plot(meso,add=T)
axis(1)
axis(2)
box()
modpres<-modpreds>threshold
plot(modpres)
plot(pts,add=T,pch=21,cex=0.4,bg=2)
plot(meso,add=T)
axis(1)
axis(2)
box()
eback<-evaluate(vars,backgr,gammod)
plot(eback,'ROC')

```




```{r,echo=FALSE,warning=FALSE,message=FALSE}

par(mfcol=c(2,2))
plot(modpres)
points(East$x,East$y,pch=21,bg=East$pb+1,cex=0.3,col=East$pb+1)
abline(v=cutpoint)
plot(meso,add=T)
box()

testmod<-gam(GAM_form,data=East)
e<-evaluate(West[West$pb==1,-1],West[West$pb==0,-1],testmod)
plot(e,'ROC')
plot(modpres)
points(West$x,West$y,ch=21,bg=West$pb+1,cex=0.3,col=West$pb+1)
abline(v=cutpoint)
plot(meso,add=T)
box()
testmod<-gam(GAM_form,data=West)
e<-evaluate(East[East$pb==1,-1],East[East$pb==0,-1],testmod)
plot(e,'ROC')


```

## Random Forest






```{r,echo=FALSE,warning=FALSE}
backgr<-na.omit(data.frame(coordinates(backgrpts),extract(preds1,backgrpts)))
vars<-na.omit(data.frame(coordinates(pts),extract(preds1,pts)))
var_names<-names(vars)[-c(1:2)]
moddata<-data.frame(cbind(pb = c(rep(1, nrow(vars)), rep(0, nrow(backgr))),rbind(vars, backgr)))
cutpoint<-median(vars$x)
East<-moddata[moddata$x<cutpoint,]
West<-moddata[moddata$x>=cutpoint,]

mod_form<-as.formula(paste("pb ~ ", paste(var_names, collapse="+"), collapse=""))
mod<-randomForest(mod_form,data=moddata)

eback<-evaluate(vars,backgr,mod)
lim<-0.9
threshold<-eback@t[which.max(eback@TPR<lim)]


```




```{r echo=FALSE,warning=FALSE}
par(mfcol=c(2,2))
modpreds<-predict(preds1,mod)
plot(modpreds)
plot(pts,add=T,pch=21,cex=0.2,bg=2)
plot(meso,add=T)
axis(1)
axis(2)
box()
modpres<-modpreds>threshold
plot(modpres)
plot(pts,add=T,pch=21,cex=0.4,bg=2)
plot(meso,add=T)
axis(1)
axis(2)
box()
eback<-evaluate(vars,backgr,gammod)
plot(eback,'ROC')

```




```{r,echo=FALSE,warning=FALSE,message=FALSE}

par(mfcol=c(2,2))
plot(modpres)
points(East$x,East$y,pch=21,bg=East$pb+1,cex=0.3,col=East$pb+1)
abline(v=cutpoint)
plot(meso,add=T)
box()

testmod<-randomForest(mod_form,data=East)
e<-evaluate(West[West$pb==1,-1],West[West$pb==0,-1],testmod)
plot(e,'ROC')
plot(modpres)
plot(meso,add=T)
points(West$x,West$y,ch=21,bg=West$pb+1,cex=0.3,col=West$pb+1)
abline(v=cutpoint)
box()
testmod<-randomForest(mod_form,data=West)
e<-evaluate(East[East$pb==1,-1],East[East$pb==0,-1],testmod)
plot(e,'ROC')


```



