FireModels.Rmd

---
title: "Lab 7 Estimating Parameters for Modules from Local Data"
author: "Steve Cumming"
date: "February 23 2017"
output: html_document
---



In this lab we integrate a landscape fire mode, data input from real landscapes as seen in Lab 1, with new modules that estimate the fire model parameters from the empirical data, with some tuning to the characteristics of the actual landscape.  

The fire model is simple extension of that we saw in Lab 5. It is based on the same percolation spread process as before, implemented using SpaDES::spread(). Fire is treated as a three stage stochastic process of ignition, escape from the cell of origin, and subsequent spread. The model was first described by Cumming et al. (1998) and more accessibly in Armstrong and Cumming (2003).

There is considerable evidence that fire suppression by initial attack directly effects the size distribution of fires by limiting the proportion of ognitions that exceed a few ha in size (Cumming 2005; Arienti et al. 2006. Typical simulation runs use cell sizes appropriate to this phenomena: < 10ha or as high as 25ha (Leroux et al. 2007). We simulate this by breaking out the first stage of the spread process as separate call to spread. Fires start from the ignition loci returned by scfmIgnition. The three modules are

* scfmIgnition
* scfmEscape
* scfmSpread

   
More details are needed here as to why we have three modules rather than use the two-stage capabilities build into spread(). 
   
   
The following code illustrates the agemodule integrated with the three stabe burn process. The three m
   
```{r scfmModule}
library(SpaDES)
library(magrittr)
savedir<-file.path("outputs")
inputDir <- file.path("inputs") 
outputDir <- file.path("outputs")

timeunit<-"year"
times <- list(start = 0, end = 50)
mapDim <- 200
defaultInterval <- 1.0
defaultPlotInterval <- 1.0
defaultInitialSaveTime <- NA #don't be saving nuffink
parameters <- list(
  .globals = list("neighbours"=8), # globals(sim)$neighbours
  .progress = list(type = "text", interval = 1),
  ageModule = list(
    initialAge=100, 
    maxAge=200, 
    returnInterval = defaultInterval, 
    startTime = times$start,
    .plotInitialTime = times$start,
    .plotInterval = defaultPlotInterval,
    .saveInitialTime = defaultInitialSaveTime, 
    .saveInterval=defaultInterval),
   scfmIgnition = list(
    pIgnition=0.0001, 
    returnInterval = defaultInterval, 
    startTime = times$start,
    .plotInitialTime = NA,
    .plotInterval = defaultPlotInterval,
    .saveInitialTime = defaultInitialSaveTime, 
    .saveInterval = defaultInterval),
  scfmEscape = list(
    p0=0.05, 
    returnInterval = defaultInterval, 
    startTime = times$start,
    .plotInitialTime = NA,
    .plotInterval = defaultPlotInterval,
    .saveInitialTime = defaultInitialSaveTime, 
    .saveInterval = defaultInterval),
  scfmSpread = list(
    pSpread=0.235, 
    returnInterval = defaultInterval, 
    startTime = times$start,
    .plotInitialTime = times$start,
    .plotInterval = defaultPlotInterval,
    .saveInitialTime = defaultInitialSaveTime, 
    .saveInterval = defaultInterval)
)

modules <- list("scfmIgnition", "scfmLandcoverInit", "scfmRegime", "ageModule", "scfmSpread", "scfmEscape") 

objects <- list(mapDim = mapDim) #note that these definitions are critical

paths <- list(
    cachePath = file.path(outputDir, "cache"),
    modulePath = file.path("modules"),
    inputPath = inputDir,
    outputPath = outputDir
)

setPaths(cachePath = paths$cachePath, 
         modulePath = paths$modulePath, 
         inputPath = paths$inputPath, 
         outputPath = paths$outputPath)
getPaths()

mySim <- simInit(times = times, params = parameters, modules = modules,
                 objects = objects, paths = paths)

set.seed(2343)
outSim <- spades(mySim, progress = FALSE)
#setwd(savedir)
```




```{r}
library(tools)
saveLazy <- function(..., file) {
  objs <- list(...)
  objs[sapply(objs, is.character)] <- sapply(objs[sapply(objs, is.character)], get, simplify=FALSE)
  names(objs) <- objs
  list2env(x = objs) %>%
    tools:::makeLazyLoadDB(., file)
  }



```




```{r model-inputs, eval=FALSE}
# download data (these are large files!)

baseDir <- "inputs"

if (TRUE){
  require("rgdal")
  require("memoise")
  require("magrittr")
  require("raster")

  tmp<-file.path(baseDir,"spread_firesize_curves.csv")
  
  if (file.exists(tmp)){
    spreadCalibrationTable<-read.table(tmp,header=T,sep=",")
  } else {
    stop("No Calibration Table")
  }
  
  if(!exists("readOGRmem")) readOGRmem <- memoise(readOGR)
  cacheLoc<-file.path("/tmp/cacheSpaDES")
  
  if (!file.exists(file.path(baseDir, "LandCoverOfCanada2005_V1_4", "LCC2005_V1_4a.tif"))) {
    url <- "ftp://ftp.ccrs.nrcan.gc.ca/ad/NLCCLandCover/LandcoverCanada2005_250m/LandCoverOfCanada2005_V1_4.zip"
    download.file(url = url,
                  destfile = file.path(baseDir, "LandCoverOfCanada2005_V1_4.zip"))
    unzip(zipfile = file.path(baseDir, "LandCoverOfCanada2005_V1_4.zip"),
          exdir = file.path(baseDir, "LandCoverOfCanada2005_V1_4"))
  }
  
  if (!file.exists(file.path(baseDir, "age", "age.tif"))) {
    dir.create(file.path(baseDir, "age"))
    url <- "ftp://ftp.daac.ornl.gov/data/nacp/NA_TreeAge//data/can_age04_1km.tif"
    download.file(url = url, destfile = file.path(baseDir, "age", "age.tif"))
  }
}
# load data
#savedir<-setwd(baseDir)
ageMapInit <- raster(file.path(baseDir, "age", "age.tif"))
lcc05 <- raster(file.path(baseDir, "LandCoverOfCanada2005_V1_4", "LCC2005_V1_4a.tif"))
lcc05CRS <- crs(lcc05)
firePointsInput<-readOGRmem(dsn="CanadianNationalFireDatabase/NFDB_point/",layer="NFDB_point_20141222")
#setwd(savedir)

#this is where you would read in your shapefile if you had one.
# Random polygon
areaKm2 <- 3e4
minX <- -1072250.2
maxX <- minX + sqrt(areaKm2*1e6)
minY <- 7438877-1.6e5
maxY <- minY + sqrt(areaKm2*1e6)
meanY <- mean(c(minY, maxY))

# Add random noise to polygon
#set.seed(5567788)
xAdd <- round(runif(1, -5e5, 1.5e6))
yAdd <- round(runif(1, 1e5, 5e5)) - xAdd/2
nPoints <- 20
betaPar <- 0.6
X = c(
  jitter(sort(rbeta(nPoints, betaPar, betaPar)*(maxX-minX)+minX)),
  jitter(sort(rbeta(nPoints, betaPar, betaPar)*(maxX-minX)+minX, decreasing = TRUE))
)
Y = c(
  jitter(sort(rbeta(nPoints/2, betaPar, betaPar)*(maxY-meanY)+meanY)),
  jitter(sort(rbeta(nPoints, betaPar, betaPar)*(maxY-minY)+minY, decreasing = TRUE)),
  jitter(sort(rbeta(nPoints/2, betaPar, betaPar)*(meanY-minY)+minY))
)

inputMapPolygon <- cbind(X+xAdd, Y+yAdd) %>% 
  Polygon %>%
  list %>%
  Polygons("s1") %>%
  list %>%
  SpatialPolygons(1L)
crs(inputMapPolygon) <- lcc05CRS


Plot(lcc05,new=TRUE)
Plot(inputMapPolygon,addTo="lcc05")
Plot(inputMapPolygon)
#downloadModule("forestAge",".") #and for other missing modules that have been written allready

```

This is how you might access instead some of the ready made study areas.

```{r eval=FALSE}
library(SpaDES)
library(spatial.tools)
library(magrittr) #implements pipes as %>%
library(archivist)
library(memoise)
library(raster)


if(!exists("readOGRmem")) readOGRmem <- memoise(readOGR)

    cacheLoc<-file.path("/tmp/cacheSpaDES")
    #archivist::createEmptyRepo(cacheLoc)
    foo<-raster(file.path(baseDir,"lccSquares", "grid15_lcc05_ascii.txt")) %>% setMinMax
    studyArea<-bbox_to_SpatialPolygons(bbox(foo))
    crs(studyArea)<-crs(foo)
    savedir<-setwd("~/Dropbox/SpaDES/Data")
    vegInput<-raster("LandCoverOfCanada2005_V1_4/LCC2005_V1_4a.tif")
    ageMapInit<-raster("age/age.tif")
    firePointsInput<-readOGRmem(dsn="CanadianNationalFireDatabase/NFDB_point/",layer="NFDB_point_20141222")
    setwd(savedir)


#readOGR<-memoise(readOGR)
#Steve: try to implement with cache for greater persistence

#Plot(foo,new=TRUE)
#Plot(studyArea,addTo="foo")
#downloadModule("forestAge",".") #and for oter missing modules that have been written allready

#read the full lcc05 with raster and extract the colours, use code now in scfmCrop

```



```{r fireRegimeSpecs echo=TRUE}
library(SpaDES)
library(spatial.tools)
library(magrittr) #implements pipes as %>%
library(archivist)
library(memoise)
library(raster)

cacheLoc<-file.path("/tmp/cacheSpaDES")
#there is supposed to be someway to get simInit to read this, too, but the method are indocumented and I can not get them to work...s'matter of specifying the names arguments of readOGR, and also in what environment or package to look for the memoised version we just created here.
if(!exists("readOGRmem")) readOGRmem <- memoise(readOGR)

dsnPath<-normalizePath(file.path(baseDir,"CanadianNationalFireDatabase","NFDB_point","")) #somebody don't like ~stevec. ahem, readOGR,
firePointsInput<-readOGRmem(dsn=dsnPath,layer="NFDB_point_20141222")

inputs <- data.frame(file=c(file.path("age", "age.tif"), 
                            file.path("LandCoverOfCanada2005_V1_4", "LCC2005_V1_4a.tif"),
                            file.path("spread_firesize_curves.csv")),
                     fun=c("raster","raster","read.csv"),
                     package=c("raster","raster","utils"),
                     objectName=c("ageMapInit","vegInput","spreadCalibrationTable"))

#dsnPath readOGR, rgdal,firePointsInput
#inputArgs(firePointsInput) <- list(dsn=dsnPath,layer="NFDB_point_20141222")

times <- list(start=0, end=10, timeunit="year")

plotInit <- NA #times$start+1

parameters <- list(.globals=list("neighbours"=8),
                   .progress=list(type=NA, interval=1), #was type="graphical"
                   scfmCrop=list(useCache=TRUE),
                   scfmLandcoverInit=list(useCache=TRUE),
                   scfmRegime=list(fireEpoch=c(1961,2010),fireCause=c("L"))
              )


modules <- list("scfmCrop","scfmLandcoverInit", "scfmRegime", "scfmDriver")
            

if(Sys.getenv("ComputerName")=="W-VIC-A105200") {
  paths <- list(modulePath=file.path("~/Documents/GitHub/BEACONs/BEACONs"),
             inputPath="M:/data",
             outputPath="c:/temp/SpaDES")
} else {
  paths <- list(modulePath=file.path("~/Dropbox/Courses/7043H16/Lab/scfmModules"),
             inputPath=baseDir,
             outputPath="/tmp/SpaDES")
}

objects <- list(#vegInput = lcc05, 
                #ageMapInit = ageMapInit,
                firePointsInput = firePointsInput,
                studyArea = inputMapPolygon,
                cacheLoc = cacheLoc
                )
                #spreadCalibrationTable = spreadCalibrationTable)

savedir<-setwd("~/Dropbox/SpaDES/Data")
try(rm(mySim))
mySim <- simInit(times=times, params=parameters, modules=modules,paths=paths, objects=objects,inputs=inputs) 
#tmpSim<- spades(mySim,debug=TRUE)
```



```{r Data2Models echo=TRUE}
library(SpaDES)
library(igraph)
library(spatial.tools)
library(magrittr) #implements pipes as %>%
library(archivist)
library(memoise)
library(raster)
savedir<-setwd("/Users/stevec/Dropbox/Courses/7043H17/Lab")
baseDir<-"/Users/stevec/Dropbox/SpaDES/Data"
cacheLoc<-file.path("/tmp/cacheSpaDES")
#there is supposed to be someway to get simInit to read this, too, but the method are undocumented and I can not get them to work...s'matter of specifying the names arguments of readOGR, and also in what environment or package to look for the memoised version we just created here.
if(!exists("readOGRmem")) readOGRmem <- memoise(readOGR)

dsnPath<-normalizePath(file.path(baseDir,"CanadianNationalFireDatabase","NFDB_point","")) #somebody don't like ~stevec. ahem, readOGR,
firePointsInput<-readOGRmem(dsn=dsnPath,layer="NFDB_point_20141222")

inputs <- data.frame(file=c(file.path("age", "age.tif"), 
                            file.path("LandCoverOfCanada2005_V1_4", "LCC2005_V1_4a.tif"),
                            file.path("spread_firesize_curves.csv"),
                            file.path("FiresN1000MinFiresSize2NoLakes.csv")),
                     fun=c("raster","raster","read.csv", "read.csv"),
                     package=c("raster","raster","utils", "utils"),
                     objectName=c("ageMapInit","vegMapInit","spreadCalibrationTable","cTable2"),
                     arguments=c(list(NULL),list(NULL),list(NULL),list(header=TRUE)))

#dsnPath readOGR, rgdal,firePointsInput
#inputArgs(firePointsInput) <- list(dsn=dsnPath,layer="NFDB_point_20141222")

times <- list(start=0, end=30, timeunit="year")

plotInit <- NA #times$start+1
defaultInterval<-1
defaultInitialSaveTime<-NA

parameters <- list(.globals=list("neighbours"=8),
                   .progress=list(type=NA, interval=1), #was type="graphical"
                   scfmCrop=list(useCache=TRUE),
                   scfmLandcoverInit=list(useCache=TRUE),
                   scfmRegime=list(fireEpoch=c(1981,2010),fireCause=c("L")),
   ageModule = list(
    initialAge=100, 
    maxAge=200, 
    returnInterval = defaultInterval, 
    startTime = times$start,
    .plotInitialTime = times$start,
    .plotInterval = defaultInterval,
    .saveInitialTime = defaultInitialSaveTime, 
    .saveInterval=defaultInterval),
   scfmIgnition = list(
    pIgnition=0.0002, 
    returnInterval = defaultInterval, 
    startTime = times$start,
    .plotInitialTime = NA,
    .plotInterval = defaultInterval,
    .saveInitialTime = defaultInitialSaveTime, 
    .saveInterval = defaultInterval),
  scfmEscape = list(
    p0=0.1, 
    returnInterval = defaultInterval, 
    startTime = times$start,
    .plotInitialTime = NA,
    .plotInterval = defaultInterval,
    .saveInitialTime = defaultInitialSaveTime, 
    .saveInterval = defaultInterval),
  scfmSpread = list(
    pSpread=0.235, 
    returnInterval = defaultInterval, 
    startTime = times$start,
    .plotInitialTime = times$start,
    .plotInterval = defaultInterval,
    .saveInitialTime = defaultInitialSaveTime, 
    .saveInterval = defaultInterval)
              )

#modules<-list("scfmParent")
modules <- list("scfmCrop","scfmLandcoverInit", "scfmRegime", "scfmDriver","ageModule", "scfmIgnition","scfmEscape", "scfmSpread") # "caribou"
            
paths <- list(modulePath=file.path("scfmModules"),
             inputPath=baseDir,
             outputPath="/tmp/SpaDES")

foo<-raster(file.path(baseDir,"lccSquares", "grid17_lcc05_ascii.txt")) %>% setMinMax
inputMapPolygon<-bbox_to_SpatialPolygons(bbox(foo))
objects <- list(firePointsInput = firePointsInput,
                studyArea = inputMapPolygon,
                cacheLoc = cacheLoc
                )

try(rm(mySim))
mySim <- simInit(times=times, params=parameters, modules=modules,paths=paths, objects=objects,inputs=inputs) 
#dev()
#tmpSim<- spades(myim,debug=TRUE)
```