Sampling Design
author: Derek Corcoran
Last update: 2015-04-30
First we load the spatial packages
library("raster", lib.loc="~/R/win-library/3.2")
library("rasterVis", lib.loc="~/R/win-library/3.2")
library("maps", lib.loc="~/R/win-library/3.2")
library("maptools", lib.loc="~/R/win-library/3.2")
library("rgdal", lib.loc="~/R/win-library/3.2")Then we read the needed rasters
PNF<- readGDAL("C:/Users/usuario/Bats_California/layers/PNF.asc")## C:/Users/usuario/Bats_California/layers/PNF.asc has GDAL driver AAIGrid
## and has 250 rows and 434 columnsPNF<-raster (PNF)
plot(PNF)
bc <- readGDAL("C:/Users/usuario/Bats_California/layers/burn_canopy.asc")## C:/Users/usuario/Bats_California/layers/burn_canopy.asc has GDAL driver AAIGrid
## and has 250 rows and 322 columnsbc<-raster (bc)
plot(bc)
bb <- readGDAL("C:/Users/usuario/Bats_California/layers/burn_basal.asc")## C:/Users/usuario/Bats_California/layers/burn_basal.asc has GDAL driver AAIGrid
## and has 250 rows and 322 columnsbb<-raster (bb)
plot(bb)
bs <- readGDAL("C:/Users/usuario/Bats_California/layers/burn_severity.asc")## C:/Users/usuario/Bats_California/layers/burn_severity.asc has GDAL driver AAIGrid
## and has 250 rows and 322 columnsbs<-raster (bs)
plot(bs)
topo <- readGDAL("C:/Users/usuario/Bats_California/layers/plumastopo.asc")## C:/Users/usuario/Bats_California/layers/plumastopo.asc has GDAL driver AAIGrid
## and has 103 rows and 151 columnstopo<-raster (topo)
plot(topo)
Vegetation_existing <- readGDAL("C:/Users/usuario/Bats_California/layers/Vegetation_existing.asc")## C:/Users/usuario/Bats_California/layers/Vegetation_existing.asc has GDAL driver AAIGrid
## and has 250 rows and 314 columnsVegetation_existing<-raster (Vegetation_existing)
plot(Vegetation_existing)
FireReturnIntervalDeparture <- readGDAL("C:/Users/usuario/Bats_California/layers/FireReturnIntervalDeparture.asc")## C:/Users/usuario/Bats_California/layers/FireReturnIntervalDeparture.asc has GDAL driver AAIGrid
## and has 250 rows and 329 columnsFireReturnIntervalDeparture<-raster (FireReturnIntervalDeparture)
plot(FireReturnIntervalDeparture)
TreatmentsStorrie <- readGDAL("C:/Users/usuario/Bats_California/layers/TreatmentsStorrie.asc")## C:/Users/usuario/Bats_California/layers/TreatmentsStorrie.asc has GDAL driver AAIGrid
## and has 271 rows and 250 columnsTreatmentsStorrie<-raster (TreatmentsStorrie)
plot(TreatmentsStorrie)
Change outlayers and extract NAs
In order to classify the raster we will get rid of unnecesary outlayers, and change NAs to 0
df.bb <- data.frame(id=c(NA,1,2,3,4,5,6,7,255), v=c(0,1,2,3,4,5,6,7,8))
bb1 <- subs(bb, df.bb,subswithNA=FALSE)
df.bs <- data.frame(id=c(NA,1,2,3,4,255), v=c(0,1,2,3,4,5))
bs1 <- subs(bs, df.bs,subswithNA=FALSE)
df.bc <- data.frame(id=c(NA,1,2,3,4,5,255), v=c(0,1,2,3,4,5,6))
bc1 <- subs(bc, df.bc,subswithNA=FALSE)Put all rasters in the same projection
roads.v <- spTransform(roads.v, CRS("+proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs"))
bb<-projectRaster(bb1, PNF)
bs<-projectRaster(bs1, PNF)
Vegetation_existing<-projectRaster(Vegetation_existing, PNF)
FireReturnIntervalDeparture<-projectRaster(FireReturnIntervalDeparture, PNF)
TreatmentsStorrie<-projectRaster(TreatmentsStorrie, PNF)Put them all in the same resolution and size
bc<-resample(bc, PNF)
bb<-resample(bb, PNF)
bs<-resample(bs, PNF)
Vegetation_existing<-resample(Vegetation_existing, PNF)
FireReturnIntervalDeparture<-resample(FireReturnIntervalDeparture, PNF)
TreatmentsStorrie<-resample(TreatmentsStorrie, PNF)
topo<-resample(topo,PNF)Prepare a distance from river/road raster
roads.v <- readOGR(dsn="C:/Users/usuario/Bats_California/layers",layer="Roads")## OGR data source with driver: ESRI Shapefile
## Source: "C:/Users/usuario/Bats_California/layers", layer: "Roads"
## with 4127 features
## It has 22 fields## Warning in readOGR(dsn = "C:/Users/usuario/Bats_California/layers", layer
## = "Roads"): Z-dimension discardedroads.v <- spTransform(roads.v, CRS("+proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs"))
plot(PNF)
lines(roads.v)
template <- PNF # this will be the template
template[] <- NA # assigns all values as NA
roads.r <- rasterize(roads.v, template, field=1)
summary(roads.r) # pixels crossed by a road have "1" ## layer
## Min. 1
## 1st Qu. 1
## Median 1
## 3rd Qu. 1
## Max. 1
## NA's 86439plot(roads.r, add=TRUE)
after rasterizing the roads we make the new raster
roaddist.r <- distance(roads.r)
class(roaddist.r)## [1] "RasterLayer"
## attr(,"package")
## [1] "raster"# Check:
plot(roaddist.r)
lines(roads.v)
#Check for correlation between rasters
burn <- stack(bc1, bs1, bb1)
AllLayers <-stack(topo,TreatmentsStorrie,FireReturnIntervalDeparture, Vegetation_existing, roaddist.r)
plot (AllLayers, colNA="black")
pairs(AllLayers)
AllLayers[is.na(AllLayers)] <- 0
plot(AllLayers)
Clasification example
Even though we shouldn’t classify using 3 layers of such high classification we will use the RasterBrick of the three burn classifications to exemplify how we will divide the area into areas of similar characteristics. Here we will ask R to use kmeans to sort the area into 3 types of habitat using the abovementioned rasterbrick:
## v.1 v.2 v.3
## [1,] 0 0 0
## [2,] 0 0 0
## [3,] 0 0 0
## [4,] 0 0 0
## [5,] 0 0 0
## [6,] 0 0 0
now with every layer
## v.1 v.2 v.3 band1.1 band1.2 band1.3 band1.4
## [1,] 0 0 0 0 0 0 0
## [2,] 0 0 0 0 0 0 0
## [3,] 0 0 0 0 0 0 0
## [4,] 0 0 0 0 0 0 0
## [5,] 0 0 0 0 0 0 0
## [6,] 0 0 0 0 0 0 0
## layer
## Min. 1
## 1st Qu. 1
## Median 4
## 3rd Qu. 5
## Max. 5
## NA's 0More info on how to do this clasification in https://geoscripting-wur.github.io/AdvancedRasterAnalysis/
Extract Random points from each habitat type
df.class.4 <- data.frame(id=c(1,2,3,4,5), v=c(NA,NA,NA,4,NA))
class4 <- subs(classes2, df.class.4,subswithNA=FALSE)
points4<-sampleRandom(class4,10, na.rm=TRUE, xy=TRUE)
df.class.3 <- data.frame(id=c(1,2,3,4,5), v=c(NA,NA,3,NA,NA))
class3 <- subs(classes2, df.class.3,subswithNA=FALSE)
points3<-sampleRandom(class3,10, na.rm=TRUE, xy=TRUE)
plot(classes2)
points (points4, col= "black")
points (points3, col= "red")