In statguy/STREM: Winter track counts

#library(devtools)
#install_github("statguy/STREM")
#install_github("ropengov/gisfin")
#install_github("statguy/Blur")
library(parallel)
library(doMC)
registerDoMC(cores=detectCores())
library(STREM)
source("~/git/STREM/setup/WTC-Boot.R")
library(ggplot2)
library(scales)
library(grid)
library(plyr)
library(gridExtra)
library(Blur)

context <- Context$new(resultDataDirectory=wd.data.results, processedDataDirectory=wd.data.processed, rawDataDirectory=wd.data.raw, scratchDirectory=wd.scratch, figuresDirectory=wd.figures)
boundaryDF <- ggplot2::fortify(FinlandWTCStudy(context=context, response="canis.lupus")$studyArea$boundary)
study <- FinlandWTCStudy$new(context=context, response="canis.lupus")
#study$studyArea$loadBoundary(thin=F, tolerance=0.001)
gray <- gray.colors(2, start=0.2, end=0.8)
labelPanel <- function(p, text, x, y) {
  g <- tableGrob(data.frame(a=text), rows=NULL, cols=NULL, gpar.corefill=gpar(fill="white", col=NA))
  p + annotation_custom(grob=g, xmin=x, ymin=y, xmax=x, ymax=y)
}

# Habitat
r <- raster::crop(study$studyArea$habitat, extent(3400000, 3410000, 7000000, 7010000))
#plot(r)
forest <- raster::calc(r, function(x) if (is.na(x) || x %in% c(0,45:255)) NA else x %in% 18:35)
#plot(forest)
coords <- raster::xyFromCell(r, cell=1:ncell(r))
x <- Blur::smoothContinuousSubsets(r=forest, coords=coords, kernel=Blur::ExponentialKernel$new(), scales=c(125, 1000), .parallel=T)
#plot(flip(rasterFromXYZ(x[,c("x","y","value")]), direction="y"), asp=1)
fileName <- study$context$getFileName(dir=study$context$processedDataDirectory, name="figure-b-1", region=study$studyArea$region)
save(x, file=fileName)
load(fileName)

# Tests:
if (F) {
r <- crop(forest, extent(forest, 10, 20, 10, 20))
x <- smoothContinuousSubsets(r=r, coords=xyFromCell(r, cell=1:ncell(r)), kernel=ExponentialKernel$new(), scales=125, .parallel=F)
p <- ggplot(x) + geom_raster(aes(x, y, fill=`125`)) + scale_fill_gradient(low=gray[1], high=gray[2]) + theme_raster() + labs(x=NULL, y=NULL) + scale_x_continuous(expand=c(0,0)) + scale_y_continuous(expand=c(0,0))
plot(r)
print(p)
xy <- xyFromCell(r, cell=1:ncell(r)); x=xy[1,1]; y=xy[1,2]
kernel <- ExponentialKernel$new()$constructKernel(r, 125)
smoothContinuousSubset(r=r, x=x, y=y, scale=125, kernel=kernel)
}

r1 <- as.data.frame(rasterToPoints(r))
colnames(r1) <- c("x","y","value")
r1$value <- factor(r1$value, levels=0:255)
r2 <- as.data.frame(rasterToPoints(forest))
colnames(r2) <- c("x","y","value")
r2$value <- factor(r2$value, levels=0:1)

p1 <- ggplot(r1) + geom_raster(aes(x, y, fill=value)) + scale_fill_grey() + theme_raster() + labs(x=NULL, y=NULL) + scale_x_continuous(expand=c(0,0)) + scale_y_continuous(expand=c(0,0))
p2 <- ggplot(r2) + geom_raster(aes(x, y, fill=value)) + scale_fill_grey() + theme_raster() + labs(x=NULL, y=NULL) + scale_x_continuous(expand=c(0,0)) + scale_y_continuous(expand=c(0,0))
p3 <- ggplot(x) + geom_raster(aes(x, y, fill=`125`)) + scale_fill_gradient(low=gray[1], high=gray[2]) + theme_raster() + labs(x=NULL, y=NULL) + scale_x_continuous(expand=c(0,0)) + scale_y_continuous(expand=c(0,0))
p4 <- ggplot(x) + geom_raster(aes(x, y, fill=`1000`)) + scale_fill_gradient(low=gray[1], high=gray[2]) + theme_raster() + labs(x=NULL, y=NULL) + scale_x_continuous(expand=c(0,0)) + scale_y_continuous(expand=c(0,0))

p1 <- labelPanel(p1, "a)", min(r1$x), max(r1$y))
p2 <- labelPanel(p2, "b)", min(r1$x), max(r1$y))
p3 <- labelPanel(p3, "c)", min(r1$x), max(r1$y))
p4 <- labelPanel(p4, "d)", min(r1$x), max(r1$y))
p <- arrangeGrob(p1, p2, p3, p4, ncol=2)
saveFigure(p, "figure-b-1a.png", width=17, height=17, dpi=200)


# Human population density

library(gisfin)
request <- gisfin::GeoStatFiWFSRequest$new()$getPopulation("vaestoruutu:vaki2005_1km")
client <- gisfin::GeoStatFiWFSClient$new(request)
population <- client$getLayer("vaestoruutu:vaki2005_1km")
population <- raster::stack(sp::SpatialPixelsDataFrame(coordinates(population), population@data, proj4string=population@proj4string))
population <- projectRaster(population$vaesto, crs="+init=epsg:2393")
r <- raster::crop(population$vaesto, extent(3400000, 3410000, 7000000, 7010000) * 10)
r[is.na(r)] <- 0
r1 <- as.data.frame(rasterToPoints(r))
colnames(r1) <- c("x","y","value")
coords <- raster::xyFromCell(r, cell=1:ncell(r))
x <- Blur::smoothContinuousSubsets(r=r, coords=coords, kernel=Blur::ExponentialKernel$new(), scales=c(2000), .parallel=T)

p5 <- ggplot(r1) + geom_raster(aes(x, y, fill=value)) + scale_fill_gradient(low=gray[1], high=gray[2]) + theme_raster() + labs(x=NULL, y=NULL) + scale_x_continuous(expand=c(0,0)) + scale_y_continuous(expand=c(0,0))
p6 <- ggplot(x) + geom_raster(aes(x, y, fill=`2000`)) + scale_fill_gradient(low=gray[1], high=gray[2]) + theme_raster() + labs(x=NULL, y=NULL) + scale_x_continuous(expand=c(0,0)) + scale_y_continuous(expand=c(0,0))


# Elevation

# Elevation model data: http://www.eea.europa.eu/data-and-maps/data/digital-elevation-model-of-europe
elevation <- raster(file.path(context$rawDataDirectory, "elevation1x1_new.tif"))
r <- raster::crop(elevation, extent(elevation, 200, 1500, 2600, 3500))
raster::projection(r) <- "+proj=laea +lat_0=52 +lon_0=20 +x_0=5071000 +y_0=3210000 +a=6378137 +b=6378137 +units=m +no_defs"
elevation <- raster::projectRaster(r, crs="+init=epsg:2393")
r <- raster::crop(elevation, extent(3400000, 3410000, 7000000, 7010000) * 10)
s <- raster(extent(3355003, 3455003, 6955246, 7055246), 100, 100, crs="+proj=laea +lat_0=52 +lon_0=20 +x_0=5071000 +y_0=3210000 +a=6378137 +b=6378137 +units=m +no_defs")
r <- raster::resample(r, s) # Changing projection does not work good enough. Have to resample to fix.
r1 <- as.data.frame(rasterToPoints(r))
colnames(r1) <- c("x","y","value")
coords <- raster::xyFromCell(r, cell=1:ncell(r))
x <- Blur::smoothContinuousSubsets(r=r, coords=coords, kernel=Blur::ExponentialKernel$new(), scales=c(2000), .parallel=T)

p7 <- ggplot(r1) + geom_raster(aes(x, y, fill=value)) + scale_fill_gradient(low=gray[1], high=gray[2]) + theme_raster() + labs(x=NULL, y=NULL) + scale_x_continuous(expand=c(0,0)) + scale_y_continuous(expand=c(0,0))
p8 <- ggplot(x) + geom_raster(aes(x, y, fill=`2000`)) + scale_fill_gradient(low=gray[1], high=gray[2]) + theme_raster() + labs(x=NULL, y=NULL) + scale_x_continuous(expand=c(0,0)) + scale_y_continuous(expand=c(0,0))


p5 <- labelPanel(p5, "e)", min(r1$x), max(r1$y))
p6 <- labelPanel(p6, "f)", min(r1$x), max(r1$y))
p7 <- labelPanel(p7, "g)", min(r1$x), max(r1$y))
p8 <- labelPanel(p8, "h)", min(r1$x), max(r1$y))
p <- arrangeGrob(p5, p6, p7, p8, ncol=2)
saveFigure(p, "figure-b-1b.png", width=17, height=17, dpi=200)

writeTable <- function(df, filename, context) {
  library(ReporteRs)
  doc <- docx()
  flextable <- FlexTable(df)
  addFlexTable(doc, flextable)
  dir <- normalizePath(context$figuresDirectory)
  writeDoc(doc, file=file.path(dir, filename))
}

models <- data.frame(Species=c("Canis Lupus"), Model=c("Intercept-only"), WAIC=c("1074.72"))
models <- rbind(models, data.frame(Species=c("Canis Lupus", "Canis Lupus"), Model=c("All PCs", "PC1 habitat 250 + PC2 all 64000"), WAIC=c("1074.78", "1072.57")))
models <- rbind(models, data.frame(Species=c("Canis Lupus", "Canis Lupus"), Model=c("All non-PCs 32000", "All non-PCs 16000"), WAIC=c("1061.45", "1067.00")))
models <- rbind(models, data.frame(Species=c("Lynx Lynx"), Model=c("Intercept-only"), WAIC=c("3922.51")))
models <- rbind(models, data.frame(Species=c("Lynx Lynx", "Lynx Lynx"), Model=c("All PCs", "PC3 all 8000"), WAIC=c("3930.95", "3914.81")))
models <- rbind(models, data.frame(Species=c("Lynx Lynx", "Lynx Lynx"), Model=c("All non-PCs 2000", "All non-PCs 8000"), WAIC=c("3917.34", "3917.35")))
writeTable(models, "table-b-1.docx", context)