Skip to content

Commit

Permalink
Updated methodology for the definition of fire danger classes (thresh…
Browse files Browse the repository at this point in the history
…olds)
  • Loading branch information
cvitolo committed Mar 10, 2017
1 parent f45a8d2 commit 0da8069
Show file tree
Hide file tree
Showing 11 changed files with 377 additions and 129 deletions.
7 changes: 4 additions & 3 deletions DESCRIPTION
Original file line number Diff line number Diff line change
Expand Up @@ -2,20 +2,22 @@ Package: caliver
Type: Package
Title: CALIbration and VERification of gridded model outputs
Version: 0.1
Date: 2017-03-10
Authors@R: c(person("Claudia", "Vitolo", email = "[email protected]", role = c("aut", "cre")),
person("Mirko", "D'Andrea", role = c("aut", "ctb")),
person("Francesca", "Di Giuseppe", role = c("ctb")))
person("Francesca", "Di Giuseppe", role = c("aut", "ctb")))
Maintainer: Claudia Vitolo <[email protected]>
URL: https://github.com/ecmwf/caliver
BugReports: https://github.com/ecmwf/caliver/issues
Description: Utility functions for the post-processing, calibration and
validation of gridded model outputs. Initial test cases include the outputs of
the following forest fire models: GEFF and RISICO.
Depends:
R (>= 3.1)
R (>= 3.1), maps
Imports:
rgdal,
ncdf4,
ggplot2,
raster,
sp,
grDevices,
Expand All @@ -27,7 +29,6 @@ Imports:
stringr,
lubridate,
rhdf5,
maps,
maptools,
grid
Suggests:
Expand Down
21 changes: 19 additions & 2 deletions NAMESPACE
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
# Generated by roxygen2: do not edit by hand

export(RiskThresholds)
export(DangerLevels)
export(ThresholdsByRegion)
export(countryMaskCrop)
export(decompressGZ)
Expand All @@ -12,13 +12,26 @@ export(getGriddedCDF)
export(getPercRiskIndex)
export(getSeason)
export(mergetime)
export(plotPDF)
export(plotPercentiles)
export(readRisicoBinary)
export(regionalBBOX)
export(regionalMask)
export(shiftMap)
import(maps)
import(ncdf4)
import(rgdal)
importFrom(ggplot2,aes)
importFrom(ggplot2,element_text)
importFrom(ggplot2,geom_density)
importFrom(ggplot2,geom_segment)
importFrom(ggplot2,geom_text)
importFrom(ggplot2,ggplot)
importFrom(ggplot2,scale_colour_manual)
importFrom(ggplot2,theme)
importFrom(ggplot2,theme_bw)
importFrom(ggplot2,xlab)
importFrom(ggplot2,ylab)
importFrom(grDevices,heat.colors)
importFrom(graphics,plot)
importFrom(grid,grid.layout)
Expand All @@ -31,7 +44,6 @@ importFrom(httr,http_error)
importFrom(httr,write_disk)
importFrom(latticeExtra,layer)
importFrom(lubridate,yday)
importFrom(maps,map)
importFrom(maptools,map2SpatialPolygons)
importFrom(raster,crop)
importFrom(raster,extent)
Expand All @@ -47,6 +59,11 @@ importFrom(rhdf5,h5read)
importFrom(rworldmap,getMap)
importFrom(sp,CRS)
importFrom(sp,sp.lines)
importFrom(stats,density)
importFrom(stats,median)
importFrom(stats,na.omit)
importFrom(stats,quantile)
importFrom(stats,sd)
importFrom(stats,uniroot)
importFrom(stringr,str_pad)
importFrom(utils,download.file)
178 changes: 178 additions & 0 deletions R/DangerLevels.R
Original file line number Diff line number Diff line change
@@ -0,0 +1,178 @@
#' Define Risk Threshold
#'
#' @description This function calculates the danger levels (VeryLow-Low-Moderate-High-VeryHigh-Extreme) for a given country.
#'
#' @param fireIndex fire index to calculate the thresholds for (default is fwi = fire weather index)
#' @param countryName string describing the country name.
#' @param baseDir this is the directory where the reanalysis data are saved.
#' @param dataDates dates of the reanalysis data to use.
#' @param fireSeasonIndex vector of indices (same length as dataDates) with TRUE if the date falls in the fire season, FALSE otherwise.
#' @param returnExtremeValues logical value, by default set to FALSE. If this is set to TRUE, it returns the yearly extreme values (useful to calculate trends).
#'
#' @return A numeric vector listing the thresholds.
#'
#' @export
#'
#' @examples
#' \dontrun{
#'
#' baseDir <- "."
#' countryName <- "Spain"
#'
#' # Define period for Reanalysis
#' dataDates <- seq.Date(from = as.Date(strptime(paste("1980", 1),
#' format="%Y %j")),
#' to = as.Date(strptime(paste("2016", 366),
#' format="%Y %j")),
#' by = "day")
#'
#' # Create an index of fire season dates
#' seasons <- getFireSeason(dataDates,
#' FSS = NULL, FSE = NULL, emisphere = "north")
#' fireSeasonIndex <- which(seasons == TRUE)
#'
#' # Calculate the thresholds
#' thresholdFWI <- DangerLevels(fireIndex = "fwi", countryName = "Spain",
#' baseDir, dataDates, fireSeasonIndex)
#' }
#'

DangerLevels <- function(fireIndex = "fwi",
countryName = "Spain",
baseDir = getwd(),
dataDates,
fireSeasonIndex,
returnExtremeValues = FALSE){

filename <- file.path(baseDir, paste0(fireIndex, "_rotated.nc"))
if (!file.exists(filename)){

message("Merging reanalysis data")

# Merge (rotated) files into a single NetCDF
mergetime(dirs = paste0(baseDir, "/rotated/", fireIndex),
varname = fireIndex,
recursive = FALSE,
outFile = filename,
outDir = baseDir)

}

filename <- file.path(baseDir, paste0(fireIndex, "_rotated_masked.grd"))
if (!file.exists(filename)){

# Load single NetCDF into a single RasterBrick
IDX <- raster::brick(file.path(baseDir, paste0(fireIndex, "_rotated.nc")))

message("Masking reanalysis data using JRC's fuel map")

# Mask using JRC fuel_map
IDX_masked <- fuelmodelMask(IDX)
raster::writeRaster(IDX_masked,
filename=filename,
bandorder='BIL', overwrite=TRUE, progress = 'text')
rm(IDX)

}else{

IDX_masked <- raster::brick(filename)

}

filename <- file.path(baseDir,
paste0(fireIndex, "_rotated_masked_",
countryName, ".grd"))
if (!file.exists(filename)){

message("Masking and cropping reanalysis data over given country")

# Mask & Crop over the country of interest
IDXcountry <- countryMaskCrop(IDX_masked, countryName)
raster::writeRaster(IDXcountry,
filename = filename,
bandorder='BIL', overwrite=TRUE, progress = 'text')

}else{

IDXcountry <- raster::brick(filename)

}

rm(IDX_masked)

# Identify & remove non-prone areas
# (do not apply this, IF there is not much difference!)
# proneAreas <- raster::calc(x = IDXspain, fun = sum, progress = "text")
# proneAreas[proneAreas == 0] <- NA
# IDXspainProne <- raster::mask(IDXspain, proneAreas) # no changes!
# w <- array(EUROIDXsummerProne)
# w <- w[!is.na(w)]
# plot(density(w), main = "IDX (removed non-prone areas) - summer - Europe")
# round(quantile(w, c(0.50, 0.70, 0.80, 0.90, 0.99)), 1)

# Identify & remove stationary areas
# (do not apply this, IF there is not much variability!)
# s1 <- raster::stack()
# for (myYear in 1980:2016){
# print(myYear)
# idx <- which(years[fireSeasonIndex] == myYear)
# IDXIdx <- raster::subset(IDXspainProne, idx)
# s1 <- raster::stack(s1, raster::calc(x = IDXIdx, fun = sd))
# }
# statAreas <- raster::calc(x = s1, fun = sd)
# hist(as.vector(statAreas))

message("Calculating thresholds of danger levels")

# Calculate extreme yearly danger
years <- lubridate::year(dataDates[fireSeasonIndex])

# Number of days per year = extreme yearly danger (assumption)
ndays <- 4

# Calculate percentile related to the above assumption
extremePercentile <- floor(x = (1-ndays/365)*100)/100
extremeValues <- c()

for (FireYear in unique(lubridate::year(dataDates))){
print(FireYear)
yearIDX <- which(years == FireYear)
sub1 <- raster::subset(IDXcountry, fireSeasonIndex)
sub2 <- raster::subset(sub1, yearIDX)
IDXyear <- quantile(na.omit(as.vector(sub2)), extremePercentile)
extremeValues <- c(extremeValues, as.numeric(IDXyear))
}

# Transform FWI threshold into Intensity (I)
# see formula 31 and 32 in http://cfs.nrcan.gc.ca/pubwarehouse/pdfs/19927.pdf
f <- function (Icomponent, extremeDanger = median(extremeValues)){

log(0.289*Icomponent) - 0.980*(log(extremeDanger))^1.546

}

# Inspect f: curve(f, from = 0, to = 1000000); abline(h = 0, lty = 3)
Icomponent <- uniroot(f = f, interval = c(0, 100000))$root
a <- Icomponent^(1/5)

# We want to get 5 danger classes
nClasses <- 5
thresholds <- c()
for (i in 1:nClasses){
# Transform back into FWI
thresholds[i] <- round(exp(1.013*(log(0.289*a^i))^0.647), 0)
# If threshold is NA, return 0!
if (is.na(thresholds[i])) thresholds[i] <- 0
}

if (returnExtremeValues == TRUE){

return(list("thresholds" = thresholds, "extremeValues" = extremeValues))

}else{

return(thresholds)

}

}
66 changes: 0 additions & 66 deletions R/RiskThresholds.R

This file was deleted.

9 changes: 6 additions & 3 deletions R/caliver-package.R
Original file line number Diff line number Diff line change
Expand Up @@ -5,23 +5,26 @@
#' @name caliver
#' @docType package
#'
#' @import maps
#' @import rgdal
#' @import ncdf4
#' @importFrom ggplot2 ggplot aes element_text geom_density geom_segment
#' @importFrom ggplot2 geom_text scale_colour_manual theme theme_bw
#' @importFrom ggplot2 xlab ylab
#' @importFrom raster raster stack mask crop rotate extent extract t plot
#' @importFrom grDevices heat.colors
#' @importFrom rworldmap getMap
#' @importFrom sp sp.lines CRS
#' @importFrom rasterVis levelplot
#' @importFrom latticeExtra layer
#' @importFrom graphics plot
#' @importFrom stats quantile
#' @importFrom httr GET authenticate write_disk http_error
#' @importFrom stringr str_pad
#' @importFrom lubridate yday
#' @importFrom rhdf5 h5read
#' @importFrom utils download.file
#' @importFrom maps map
#' @importFrom maptools map2SpatialPolygons
#' @importFrom grid grid.newpage pushViewport viewport grid.layout
#'
#' @importFrom stats median na.omit sd uniroot density quantile
#'
NULL
14 changes: 12 additions & 2 deletions R/getFireSeason.R
Original file line number Diff line number Diff line change
Expand Up @@ -36,13 +36,23 @@ getFireSeason <- function(DATES, FSS = NULL, FSE = NULL, emisphere = "north"){
}

if (emisphere == "north"){
season <- ifelse(d >= FSS & d <= FSE, 1, 0)
season <- ifelse(d >= FSS & d <= FSE, TRUE, FALSE)
}else{
FirstJan <- as.Date("2012-01-01", format = "%Y-%m-%d")
LastDec <- as.Date("2012-12-31", format = "%Y-%m-%d")
season <- ifelse(d >= FirstJan & d <= FSE | d >= FSS & d <= LastDec, 1, 0)
season <- ifelse(d >= FirstJan & d <= FSE | d >= FSS & d <= LastDec,
TRUE, FALSE)
}

# if (zone == "tropics"){
# season <- ifelse(d >= FSS & d <= FSE, TRUE, FALSE)
# }else{
# FirstJan <- as.Date("2012-01-01", format = "%Y-%m-%d")
# LastDec <- as.Date("2012-12-31", format = "%Y-%m-%d")
# season <- ifelse(d >= FirstJan & d <= FSE | d >= FSS & d <= LastDec,
# TRUE, FALSE)
# }

return(season)

}
Loading

0 comments on commit 0da8069

Please sign in to comment.