Part_2.Rmd

# **Habitat Suitability and Distribution Models**
### with Applications in R
\
**by A. Guisan (1), W. Thuiller (2), N.E. Zimmermann (3) **,\
\
with contribution by V. Di Cola, D. Georges and A. Psomas\
\
_(1) University of Lausanne, Switzerland_\
_(2) CNRS, Université Grenoble Alpes, France_\
_(3) Swiss Federal Research Institute WSL, Switzerland_\


#### Cambridge University Press

http://www.cambridge.org/gb/academic/subjects/life-sciences/quantitative-biology-biostatistics-and-mathematical-modellin/habitat-suitability-and-distribution-models-applications-r

*Citation:* 
@book{
  title={Habitat Suitability and Distribution Models: With Applications in R},
  author={Guisan, A. and Thuiller, W. and Zimmermann, N.E.},
  isbn={9780521758369},
  series={Ecology, Biodiversity and Conservation},
  year={2017},
  publisher={Cambridge University Press}
}

*If you use any of these figures and code examples in a presentation or lecture, somewhere in your set of slides we would really appreciate if you please add the paragraph: "Some of the figures in this presentation are taken from "Habitat Suitability and Distribution Models: with applications in R"  (CUP, 2017) with permission from the authors: A. Guisan, W. Thuiller and N.E. Zimmerman " 
If you wish to use any of these figures in a publication, you must get permission from CUP, and each figure must be accompanied by a similar acknowledgement.*



# Part II "Data Acquisition, Sampling Design, and Spatial Scales"
# Chapter 6: Environmental predictors: issues of processing and selection
## Performing simple GIS analyses in R
### Introduction

```{r, message=FALSE,warning=FALSE}
library(sp)
library(rgdal)
library(raster)
library(dismo)
library(maptools)
```


Set Working Directory
```{r setwd}
setwd("~/data")
```



### Loading the data and initial exploration
```{r, message=FALSE,warning=FALSE}
bio3    <- raster("raster/bioclim/current/grd/bio3")
bio7    <- raster("raster/bioclim/current/grd/bio7")
bio11   <- raster("raster/bioclim/current/grd/bio11")
bio12   <- raster("raster/bioclim/current/grd/bio12")
bbox(bio7)
ncol(bio7); nrow(bio7) ; res(bio7)
elev <- raster("raster/topo/GTOPO30.tif")
projection(elev) <- "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"
bbox(elev); ncol(elev); nrow(elev); res(elev)
```




### Resampling, spatial alignment and indices
```{r, message=FALSE,warning=FALSE,fig.height=8,fig.width=8}
# to create the new resampled file
elev1 <- resample(elev,bio7,method="ngb")
```

```{r variables 6.1, message=FALSE,warning=FALSE,fig.height=8,fig.width=6}
par(mfrow=c(3,1)) 
plot(elev1,col=rev(topo.colors(50)),main="Elevation")
plot(bio3,col=heat.colors(100),main="Bio.3")
plot(bio11,col=rainbow(100),main="Bio.11")
par(mfrow=c(1,1))
```

```{r cor_var 6.2, message=FALSE,warning=FALSE,fig.height=10,fig.width=10}
par(mfrow=c(2,2))
plot(bio3,bio12,xlab="bio3",ylab="bio12",col="gray55")
plot(bio3,bio7,xlab="bio3",ylab="bio7",col="gray55")
plot(bio3,bio11,xlab="bio3",ylab="bio11",col="gray55")
plot(bio3,elev1,xlab="bio3",ylab="elevation",col="gray55")
par(mfrow=c(1,1))
```


### Working with contours and lines
```{r globe_elev 6.3, message=FALSE,warning=FALSE,fig.height=10,fig.width=10}
plot(elev1, main="Elevation Contours")
contour(elev, nlevels=7, levels=c(0, 500, 1000, 1500, 2000, 3000, 4000, 5000), add=TRUE, labels="", lwd=.2)
```


```{r SA_elev 6.4, message=FALSE,warning=FALSE}
elev_sa <- crop(elev, extent(-85,-30,-60,15))
elev_na <- crop(elev, extent(-188,-50,15,90))
plot(elev_sa, main="Elevation Contours South America ")
contour(elev_sa, nlevels=7, levels=c(0, 500, 1000, 1500, 2000, 3000, 4000, 5000), add=TRUE, labels="", lwd=.2)
```


```{r}
iso<-rasterToContour(elev, nlevels=7, levels=c(0, 500, 1000, 1500, 2000, 3000, 4000, 5000))
writeOGR(iso, dsn="vector/globe", layer="isolines", "ESRI Shapefile",check_exists=T, overwrite_layer=T)
iso1<- shapefile("vector/globe/isolines.shp")
```

```{r NA_elev, message=FALSE,warning=FALSE}
plot(elev_na,col=terrain.colors(40),main="Elevation Contours") 
lines(iso1, lwd=0.2)
```


### Raster analyses of type "Global" 
```{r 6.5, message=FALSE,warning=FALSE,fig.height=10,fig.width=8}
lat<-raster("raster/other/latitude.tif")
lon<-raster("raster/other/longitude.tif")
tcold<-4.138e+02 + (-3.624e-02 * elev1) + (-8.216e+00 * abs(lat)) + 
  (-1.794e+00 * abs(lon)) + (7.122e-03 * lon^2)
diff_obs_model_temp <-bio11 - tcold
par(mfrow=c(2,1)) 
plot(tcold, col=rev(rainbow(100))[20:100],main="Modelled mean 
     temperature of the coldest quarter")
contour(elev, nlevels=7, levels=c(0, 500, 1000, 1500, 2000, 3000, 4000, 
                                  5000), add=T, labels="", lwd=.3)
plot(diff_obs_model_temp, col=rev(rainbow(100))[20:100],main= 
       "Difference between modelled and observed temperatures")
contour(elev, nlevels=7, levels=c(0, 500, 1000, 1500, 2000, 3000, 4000, 
                                  5000), add=T, labels="", lwd=.3)
par(mfrow=c(1,1))

```


### Raster analyses of type focal and terrain analyses 
```{r topExp 6.6, message=FALSE,warning=FALSE}
tmp<-crop(elev, extent(-85,-30,-60,15))
elev_sa10<-aggregate(tmp, fact=10, fun=mean, expand=TRUE, na.rm=TRUE)
w <- matrix(rep(1,225), nr=15,nc=15)
TopEx <- elev_sa10 - focal(elev_sa10, w=w, fun=mean,na.rm=TRUE)
plot(TopEx,col=gray.colors(200),main="Topographic Exposure over South America")
contour(elev_sa10, nlevels=7, levels=c(0, 500, 1000, 1500, 2000, 3000, 4000, 5000), add=T, labels="", lwd=.1)
```

```{r hillshade 6.7, message=FALSE,warning=FALSE}
#slope <- terrain(elev, opt="slope")
#aspect <- terrain(elev, opt="aspect")
#hillshade <- hillShade(slope, aspect, 30, 315)
#writeRaster(hillshade, "raster/hillshade.tif", overwrite=T)
hillshade<-raster("raster/topo/hillshade.tif")
plot_extent<-extent(-124,-66,24,50)
hillsh_na<-crop(hillshade, extent(-188,-50,15,90))
plot(elev_na,col=terrain.colors(100),alpha=.5,add=T,ext=plot_extent)
dem.c<-colorRampPalette(c("aquamarine", "lightgoldenrodyellow","lightgoldenrod", "yellow", "burlywood", "burlywood4", "plum", "seashell"))
cols<-dem.c(100)
cols<-paste(cols,"99",sep="",collate="")
#cols<-paste(gray.colors(100),"99",sep="",collate="")
plot(hillsh_na, col=grey(0:100/100), legend=FALSE, axes=F, ext=plot_extent)
plot(elev_na, col=cols,add=T, ext=plot_extent)
```

### Stacking grids to a grid stack
```{r world_stk 6.8, message=FALSE,warning=FALSE,fig.height=9,fig.width=9}
world.stk <- stack(elev1,bio3,bio7,bio11,bio12)
summary(world.stk)
plot(world.stk[[2:5]], col=rainbow(100,start=.0,end=.8))
world.stk3 <- aggregate(world.stk,fact=2, method="mean")
```

### Stacking grids to a grid stack
```{r, message=FALSE,warning=FALSE,fig.height=5,fig.width=9}
prec_yearly_usa <- raster("raster/prism/prec_30yr_normal_annual.asc")
tave_yearly_usa <- raster("raster/prism/tave_30yr_normal_annual.asc") 
extent(prec_yearly_usa)
projection(prec_yearly_usa)
prec_yearly_usa_wgs <- projectRaster(prec_yearly_usa, res=0.025, crs="+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0",method="bilinear")
tave_yearly_usa_wgs <- projectRaster(tave_yearly_usa, res=0.025,crs="+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0",method="bilinear")
```

### Importing and overlaying species data
```{r, message=FALSE,warning=FALSE}
pinus_edulis <- read.table("tabular/species/pinus_edulis_occ.csv", sep=",",header=TRUE)
class(pinus_edulis)                             
coordinates(pinus_edulis) <- c("lon", "lat")
head(coordinates(pinus_edulis))
projection(pinus_edulis) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

library(dismo)
pinus_edulis  <- gbif('pinus', 'edulis', download=T, geo=T, sp=T, removeZeros=T)
names(pinus_edulis)[1]<-"dwnld.date"

write.table(data.frame(pinus_edulis@coords, pinus_edulis@data), 
            "tabular/species/p_edulis.txt", sep="," ,row.names=FALSE)
write.csv(data.frame(pinus_edulis@coords, pinus_edulis@data), 
          "tabular/species/p_edulis.csv")
```

```{r pinus 6.9, fig.height=5,fig.width=9}
pts.clim<-extract(world.stk, pinus_edulis, method="bilinear")
pin_edu.clim<-data.frame(cbind(coordinates(pinus_edulis), pts.clim, pinus_edulis@data))
coordinates(pin_edu.clim)<-c("lon","lat")
map.ext<-extent(-120,-100,30,44)
plot(hillsh_na, col=grey(0:100/100), legend=FALSE, axes=F,ext=map.ext)
plot(elev_na,col=cols,add=T,ext=map.ext)
plot(pinus_edulis, pch=16, cex=.5, add=T)
```


### Generating a uniform spatial data structure for modelling and analysis
```{r 6.10, message=FALSE,warning=FALSE}
pts.cal<-read.table("tabular/species/cal.txt")
pts.eva<-read.table("tabular/species/eva.txt")
plot(pts.cal[which(pts.cal$VulpesVulpes==0),1:2],pch=15,cex=.3, col="grey50",xlab="Longitude",ylab="Latitude")
points(pts.eva[which(pts.eva$VulpesVulpes==0),1:2],pch=15,cex=.3, col="grey85")
points(pts.cal[which(pts.cal$VulpesVulpes==1),1:2],pch=16,cex=.4, col="firebrick3")
points(pts.eva[which(pts.eva$VulpesVulpes==1),1:2],pch=16,cex=.4, col="seagreen3")
```

```{r 6.11, message=FALSE,warning=FALSE}
pts.cal.ovl<-cbind(pts.cal[,8],extract(world.stk[[2:5]],pts.cal[,1:2]))
pts.eva.ovl<-cbind(pts.eva[,8],extract(world.stk[[2:5]],pts.eva[,1:2]))
pts.cal.ovl<-data.frame(na.omit(pts.cal.ovl))
pts.eva.ovl<-data.frame(na.omit(pts.eva.ovl))
names(pts.cal.ovl)[1]<-"Vulpes.vulpes"
names(pts.eva.ovl)[1]<-"Vulpes.vulpes"
vulpes.full <- glm(Vulpes.vulpes~bio3+I(bio3^2)+bio7+I(bio7^2)+bio11+ I(bio11^2)+bio12+I(bio12^2), family="binomial", data=pts.cal.ovl)
vulpes.step <- step(vulpes.full, direction="both", trace=F)
library(ecospat)
ecospat.adj.D2.glm(vulpes.full)
ecospat.adj.D2.glm(vulpes.step)
summary(vulpes.step)
vulpes.map<-predict(world.stk,vulpes.step, type="response")
plot(vulpes.map, col=rev(heat.colors(10)), main="Predicted distribution: Vulpes vulpes")
points(pts.cal[which(pts.cal$VulpesVulpes==1),1:2], pch=15, cex=.25)
```

## RS-based Predictors
### Importing, resampling and grid stacking
```{r 6.12, message=FALSE,warning=FALSE,fig.height=6,fig.width=10}
Cantons <-shapefile("vector/swiss/Swiss_Cantons.shp")
Zurich<-Cantons[Cantons$NAME=="ZUERICH",]
band1_blue<-raster("raster/landsat/L7_194027_2001_08_24_B10.TIF")
band2_green<-raster("raster/landsat/L7_194027_2001_08_24_B20.TIF")
band3_red<-raster("raster/landsat/L7_194027_2001_08_24_B30.TIF")
band4_nir<-raster("raster/landsat/L7_194027_2001_08_24_B40.TIF")
band5_swir1<-raster("raster/landsat/L7_194027_2001_08_24_B50.TIF")
band7_swir2<-raster("raster/landsat/L72194027_2001_08_24_B70.TIF")
band1_blue_crop<-crop(band1_blue,extent(Zurich))
band2_green_crop<-crop(band2_green,extent(Zurich))
band3_red_crop<-crop(band3_red,extent(Zurich))
band4_nir_crop<-crop(band4_nir,extent(Zurich))
band5_swir1_crop<-crop(band5_swir1,extent(Zurich))
band7_swir2_crop<-crop(band7_swir2,extent(Zurich))
L7_010824<-brick(band1_blue_crop, band2_green_crop, band3_red_crop,band4_nir_crop, band5_swir1_crop,band7_swir2_crop) 
tmp<-stack(band1_blue,band2_green,band3_red,band4_nir,band5_swir1,band7_swir2)
L7_010824<-crop(tmp, extent(Zurich))
names(L7_010824)<-c("band1_blue","band2_green","band3_red","band4_nir","band5_swir1","band7_swir2")
names(L7_010824)
par(mfrow=c(1,3))
plotRGB(L7_010824,3,2,1,stretch="lin")
plotRGB(L7_010824,4,3,2,stretch="lin")
plotRGB(L7_010824,6,4,3,stretch="lin")
par(mfrow=c(1,1))
```

### Data Processing for ecological analyses
```{r, message=FALSE,warning=FALSE,fig.height=10,fig.width=10}
NDVI <- (L7_010824$band4_nir - L7_010824$band3_red)/(L7_010824$band4_nir + L7_010824$band3_red)
NDWI <- (L7_010824$band4_nir - L7_010824$band5_swir1)/(L7_010824$band4_nir + L7_010824$band5_swir1)
SR <- L7_010824$band4_nir / L7_010824$band3_red
L7_010824_sp <- as(L7_010824, "SpatialGridDataFrame")
L7_010824_sp1 <- L7_010824_sp[1]
L7_010824_sp2 <- L7_010824_sp[2]
L7_010824_sp3 <- L7_010824_sp[3]
L7_010824_sp4 <- L7_010824_sp[4]
L7_010824_sp5 <- L7_010824_sp[5]
L7_010824_sp7 <- L7_010824_sp[6]
library(landsat)
L7_010824_refl_sp1 <- radiocorr(L7_010824_sp1, Grescale=0.76282,Brescale=-1.52, sunelev= 48.29, edist=ESdist("2011-08-24"), Esun=1957, method="apparentreflectance")
L7_010824_refl_sp2 <- radiocorr(L7_010824_sp2, Grescale=1.44251,Brescale=-2.84, sunelev= 48.29, edist=ESdist("2011-08-24"), Esun=1826, method="apparentreflectance")
L7_010824_refl_sp3 <- radiocorr(L7_010824_sp3, Grescale=1.03988,Brescale=-1.17, sunelev= 48.29, edist=ESdist("2011-08-24"), Esun=1554, method="apparentreflectance")
L7_010824_refl_sp4 <- radiocorr(L7_010824_sp4, Grescale=0.87258,Brescale=-1.51, sunelev= 48.29, edist=ESdist("2011-08-24"), Esun=1036, method="apparentreflectance")
L7_010824_refl_sp5 <- radiocorr(L7_010824_sp5, Grescale=0.11988,Brescale=-0.37, sunelev= 48.29, edist=ESdist("2011-08-24"), Esun=215, method="apparentreflectance")
L7_010824_refl_sp7 <- radiocorr(L7_010824_sp7, Grescale=0.06529,Brescale=-0.15, sunelev= 48.29, edist=ESdist("2011-08-24"), Esun=80.67,method="apparentreflectance")
L7_010824_tc <- tasscap("L7_010824_refl_sp",  sat = 7)
L7_010824_Brightness <- raster(L7_010824_tc[[1]])
L7_010824_Greenness <- raster(L7_010824_tc[[2]])
L7_010824_Wetness   <- raster(L7_010824_tc[[3]])
L <- 0.5
SAVI <- ((raster(L7_010824_refl_sp4) -  raster(L7_010824_refl_sp3))/(raster(L7_010824_refl_sp4) +  raster(L7_010824_refl_sp3) + L) )* (1+L)

hill_250m_utm <- raster("raster/topo/hill_250m_utm.tif")

```
```{r 6.13, fig.height=10,fig.width=10}
par(mfcol=c(2,2))
plot(hill_250m_utm,col=grey(0:100/100), legend=FALSE, axes=F,ext=extent(SAVI), main ="NDVI")
plot(NDVI,col=rev(terrain.colors(20,alpha=0.6)),add=T)
plot(hill_250m_utm,col=grey(0:100/100), legend=FALSE, axes=F, ext=extent(SAVI), main ="SR")
plot(SR,col=rev(terrain.colors(20,alpha=0.6)),add=T)
plot(hill_250m_utm,col=grey(0:100/100), legend=FALSE, axes=F, ext=extent(SAVI), main ="NDWI")
plot(NDWI,col=rev(topo.colors(20,alpha=0.6)),add=T)
plot(hill_250m_utm,col=grey(0:100/100), legend=FALSE, axes=F, ext=extent(SAVI), main ="SAVI")
plot(SAVI,col=rev(terrain.colors(20,alpha=0.6)),add=T)
par(mfcol=c(1,1))
```

```{r 6.14, message=FALSE,warning=FALSE,fig.height=5,fig.width=10}
ygb.c <-colorRampPalette(c("yellow","#7FFF7F","forestgreen","deepskyblue4","#00007F"))
par(mfcol=c(1,3))
plot(hill_250m_utm,col=grey(0:100/100), legend=FALSE, axes=F, ext=extent(SAVI), main ="Brightness")
plot(L7_010824_Brightness,col=rev(paste(ygb.c(20),"B3", sep="")),add=T)
plot(hill_250m_utm,col=grey(0:100/100), legend=FALSE, axes=F, ext=extent(SAVI), main ="Greenness")
plot(L7_010824_Greenness,col=rev(terrain.colors(20,alpha=0.6)),add=T)
plot(hill_250m_utm,col=grey(0:100/100), legend=FALSE, axes=F, ext=extent(SAVI), main ="Wetness")
plot(L7_010824_Wetness,col=rev(topo.colors(20,alpha=0.6)),add=T)
par(mfcol=c(1,1))
```

## Properties and selection of variables
### Correlation, collinearity and variance inflation
```{r 6.15, message=FALSE,warning=FALSE,fig.height=8,fig.width=8}
data<-read.csv("tabular/bioclim/current/bioclim_table.csv", header = TRUE, sep = ",")
ecospat.cor.plot(data[,4:8])
```
```{r}
library(usdm)
vif(data[,4:8])
vif(data[,c(4,6:8)])
vifstep(world.stk[[2:5]])
vif(world.stk[[2:5]])
vifcor(data[,4:8], th=.7)
```

### Variable Pre-selection
```{r, message=FALSE,warning=FALSE}
cor(data$bio12,data$bio7)
cor(log(data$bio12+.001),data$bio7)
```

# Chapter 7:	Species data - issues of acquisition and design
## Spatial Autocorrelation and Pseudo-replicates
```{r 7.1, message=FALSE,warning=FALSE,fig.height=4,fig.width=10}
library(ape)
xy <- pts.cal[,1:2]
dists <- as.matrix(dist(xy))
dists.inv <- 1/dists
diag(dists.inv) <- 0
Moran.I(vulpes.step$residuals, dists.inv)
library(ncf)
rsd<-vulpes.step$residuals
rnd<-sample(1:length(rsd),500,replace=T)
spat.cor<-correlog(xy[rnd,1],xy[rnd,2],rsd[rnd],increment=2,resamp=10)
par(mfcol=c(1,2))
plot(spat.cor$mean,spat.cor$corr,ylim=c(-.05,.05),xlim=c(0,50), pch=16,col="firebrick3", ylab="Correlation", xlab="Distance class", main="Spatial Correlogram", font.lab=2)
lines(spat.cor$mean,spat.cor$corr,col="firebrick3")
abline(0,0,col="grey50",lty=3)
plot(xy[order(rsd),], pch=15, col=rev(heat.colors(5900)), cex=.3, main="Residuals from Vulpes vulpes GLM ", xlab="Latitude", ylab="Longitude", font.lab=2)
par(mfcol=c(1,1))
```

## Sample Size, Prevalence and Sample accuracy
```{r, message=FALSE,warning=FALSE,fig.height=10,fig.width=10}
library(ecospat)
library(PresenceAbsence)
vulvul.pa<-cbind(pts.cal.ovl,runif(dim(pts.cal.ovl)[1],1,100))
names(vulvul.pa)[6]<-"srt"
vulvul.pa<-vulvul.pa[order(vulvul.pa$srt),]
vulvul.p<-cbind(vulvul.pa[which(vulvul.pa[,1]==1),])
vulvul.a<-cbind(vulvul.pa[which(vulvul.pa[,1]==0),])

yb1<-c(10,20,40,60,80,100,125,150,200,250,300,400,600,800,1000,1500,length(vulvul.p$srt))
pr.qual<-data.frame(matrix(data=NA,nrow=length(yb1),ncol=12))
names(pr.qual)=c("adjD2.f","adjD2.s","AUC.f","AUC.s","AUC.x","Kappa.f","Kappa.s","Kappa.x","TSS.f","TSS.s","TSS.x","Prev")
pr.qual[12]<-yb1/(yb1+length(vulvul.a$srt))

for (i in 1:length(yb1)){
  paok<-rbind(vulvul.p[1:yb1[i],],vulvul.a)
  rownames(paok)<-1:dim(paok)[1]
  # Full
  paok1f<-glm(Vulpes.vulpes~bio3+I(bio3^2)+bio11+I(bio11^2)+bio12+
                I(bio12^2),family="binomial",data=paok)
  paok0<-predict(paok1f,paok,type="response")
  pr.qual[i,1]<-ecospat.adj.D2.glm(paok1f)
  tmp1 <- data.frame(1:length(paok0),paok[,1],paok0)
  names(tmp1) <- c("ID","Observed","Predicted")
  pr.qual[i,3]<-auc(tmp1)$AUC
  pr.qual[i,6]<-ecospat.max.kappa(paok0,paok[,1])[[2]][1,2] 
  pr.qual[i,9]<-ecospat.max.tss(paok0,paok[,1])[[2]][1,2] 
  
  # Step
  paok1s<-step(paok1f,direction="both",trace=F)
  paok0<-predict(paok1s,paok,type="response")
  pr.qual[i,2]<-ecospat.adj.D2.glm(paok1s)
  tmp1 <- data.frame(1:length(paok0),paok[,1],paok0)
  names(tmp1) <- c("ID","Observed","Predicted")
  pr.qual[i,4]<-auc(tmp1)$AUC
  pr.qual[i,7]<-ecospat.max.kappa(paok0,paok[,1])[[2]][1,2] 
  pr.qual[i,10]<-ecospat.max.tss(paok0,paok[,1])[[2]][1,2] 
  
  # Xval
  paok1x<-ecospat.cv.glm(paok1s)
  tmp1 <- data.frame(1:length(paok0),paok[,1],paok1x$predictions)
  names(tmp1) <- c("ID","Observed","Predicted")
  pr.qual[i,5]<-auc(tmp1)$AUC
  pr.qual[i,8]<-ecospat.max.kappa(paok1x$predictions,paok[,1])[[2]][1,2] 
  pr.qual[i,11]<-ecospat.max.tss(paok1x$predictions,paok[,1])[[2]][1,2] 
}
```
```{r 7.2, fig.height=10,fig.width=10}
plot(pr.qual$Prev,pr.qual$Kappa.f,ty="l",lwd=5,col="#00FF00B4",ylim=c(0,1.0),xlim=c(0,.5),xlab="Prevalence", ylab="Model Quality",main="Prevalence Effects")
points(pr.qual$Prev,pr.qual$Kappa.s,ty="l",lwd=5,col="#00CD00B4",lty=3)
points(pr.qual$Prev,pr.qual$Kappa.x,ty="l",lwd=5,col="#008B00B4",lty=2)
points(pr.qual$Prev,pr.qual$TSS.f,ty="l",lwd=5,col="#ADD8E6B4")
points(pr.qual$Prev,pr.qual$TSS.s,ty="l",lwd=5,col="#9FB6CDB4",lty=3)
points(pr.qual$Prev,pr.qual$TSS.x,ty="l",lwd=5,col="#0000FFB4",lty=2)
points(pr.qual$Prev,pr.qual$AUC.f,ty="l",lwd=5,col="#EE2C2CB4")
points(pr.qual$Prev,pr.qual$AUC.s,ty="l",lwd=5,col="#CD2626B4",lty=3)
points(pr.qual$Prev,pr.qual$AUC.x,ty="l",lwd=5,col="#8B1A1AB4",lty=2)
legend(.355,0.45,c("Kappa full","Kappa step","Kappa xval", "TSS full",
                   "TSS step", "TSS xval", "AUC full", "AUC step", "AUC xval"),
       lty=c(1,3,2,1,3,2,1,3,2),lwd=c(8),col=c("#00FF00B4","#00CD00B4",
                                               "#008B00B4","#ADD8E6B4","#9FB6CDB4","#0000FFB4","#EE2C2CB4",
                                               "#CD2626B4","#8B1A1AB4"))

```

## Sampling Design and Data Collection
### Preparing stratifications for spatial sampling design
```{r, message=FALSE,warning=FALSE,fig.height=6,fig.width=10}
library(ecospat)
library(classInt)
usa <- shapefile("vector/usa/USA_states.shp")
usa_contin <- usa[usa$STATE_NAME != "Alaska" & usa$STATE_NAME != "Hawaii", ]
empty_raster <- raster(bio3)
#usa_raster <- rasterize(usa_contin, empty_raster, field="DRAWSEQ")
usa_raster <- rasterize(usa_contin, empty_raster)
bio3.us <- crop(mask(bio3, usa_raster), extent(usa_contin))
bio12.us <- crop(mask(bio12, usa_raster), extent(usa_contin))
B3.rcl<-ecospat.rcls.grd(bio3.us,9) 
B12.rcl<-ecospat.rcls.grd(bio12.us,9)
B3B12.comb <- B3.rcl+ B12.rcl*10

cspan<-maxValue(B3B12.comb)-minValue(B3B12.comb)
yb<-rainbow(100)[round(runif(cspan,.5,100.5))]
```
```{r 7.3, fig.height=6,fig.width=10}
par(mfcol=c(1,2))
hist(B3B12.comb,breaks=100,col=heat.colors(cspan), main="Histogram values")
plot(B3B12.comb,col=yb,main="Stratified map", asp=1)
#click(B3B12.comb,n=5,type="p",xy=T)
par(mfcol=c(1,1))
```

### Spatial sampling design using built-in functions in *rgdal*
```{r 7.4, message=FALSE,warning=FALSE,fig.height=10,fig.width=10}
paok <- as(B3B12.comb, "SpatialPixelsDataFrame")
s.rand<- spsample(paok,n=100,type="random")
s.strt<- spsample(paok,n=100,type="stratified",cells=3)
s.regl<- spsample(paok,n=100,type="regular")
s.nona<-spsample(paok,n=100,type="nonaligned") 


par(mfcol=c(2,2))
plot(B3B12.comb,main="random",col=rev(terrain.colors(25)))  
points(s.rand, pch=3, cex=.5)
plot(B3B12.comb,main="stratified",col=rev(terrain.colors(25)))  
points(s.strt, pch=3, cex=.5)
plot(B3B12.comb,main="nonaligned",col=rev(terrain.colors(25)))  
points(s.nona, pch=3, cex=.5)
plot(B3B12.comb,main="regular",col=rev(terrain.colors(25)))  
points(s.regl, pch=3, cex=.5)
par(mfcol=c(1,1))
```



### Random, environmentally stratified sampling design
```{r 7.5, message=FALSE,warning=FALSE,fig.height=9,fig.width=5}
envstrat_equ<- ecospat.recstrat_regl(B3B12.comb,150)
envstrat_prp<- ecospat.recstrat_prop(B3B12.comb, 150)

par(mfcol=c(1,2))
plot(B3B12.comb,main="Proportional Sampling", col=rev(terrain.colors(25)))
points(envstrat_prp$x,envstrat_prp$y,pch=16,cex=.4,col=2)
plot(B3B12.comb,main="Equal Sampling",col=rev(terrain.colors(25)))
points(envstrat_equ$x,envstrat_equ$y,pch=16,cex=.4,col=2)
par(mfcol=c(1,1))
```

```{r 7.6, fig.height=9,fig.width=5}
par(mfrow=c(2,1))
barplot(table(envstrat_prp$class),col="firebrick", 
        main="Proportional point allocation")
barplot(table(envstrat_equ$class),col="slategray4",
        main="Equal point allocation")
par(mfrow=c(1,1))

```



### Sampling designs along linear features
```{r 7.7, message=FALSE,warning=FALSE,fig.height=4,fig.width=8}
dem_globe <- raster("raster/topo/GTOPO30.tif")
dem_usa <- crop(dem_globe, usa_contin)
dem_usa_10km <- aggregate(dem_usa, 10, fun=mean)

empty_raster <- raster(dem_usa_10km)
#usa_raster <- rasterize(usa_contin, empty_raster, field="DRAWSEQ")
usa_raster <- rasterize(usa_contin, empty_raster)
dem_usa_masked <- mask(dem_usa_10km, usa_raster)
iso_1000m<-rasterToContour(dem_usa_masked, nlevels=1, levels=c(1000))
plot(dem_usa_masked,col=dem.c(100),main="Elevation Contours at 1000m") 
lines(iso_1000m, lwd=1.3, col=2)
```

```{r 7.8, message=FALSE,warning=FALSE,fig.height=8,fig.width=6}
l.rand<-spsample(iso_1000m,n=150,type="random")
l.strt<-spsample(iso_1000m,n=150,type="stratified", cells=50)
l.regl<-spsample(iso_1000m,n=150,type="regular")
par(mfcol=c(3,1))
# plot(dem_usa_masked,col=dem.c(100),main="Random Sampling")
plot(dem_usa_masked,col= gray.colors(100),main="Random Sampling")
points(l.rand, pch=3, cex=.5)
# plot(dem_usa_masked,col=dem.c(100),main="Stratified Sampling")
plot(dem_usa_masked,col=gray.colors(100),main="Stratified Sampling")
points(l.strt, pch=3, cex=.5)
# plot(dem_usa_masked,col=dem.c(100),main="Regular Sampling")
plot(dem_usa_masked,col=gray.colors(100),main="Regular Sampling")
points(l.regl, pch=3, cex=.5)
par(mfcol=c(1,1))
```