Rockies Fire Analysis
Stephen Huysman
2023-05-16
Settings
library(tidyverse)
library(multidplyr)
library(zoo)
library(slider)
library(pROC)
plot_theme <- theme_bw()
cores <- parallel::detectCores() - 4
cluster <- new_cluster(cores)
window <- 7 ## Window of days to calculate rolling values
vars <- c("RD", "VPD", "T", "SOIL", "AWSSM", "RAIN", "AET", "D", "GDD")
region_name <- "Middle Rockies" ## Set region name, should be same as L3 Ecoregion name
settings <- tribble(
~region, ~params,
"Southern Rockies", list(data_dir = "/home/steve/OneDrive/nothern-rockies-dryness",
out_dir = "/media/steve/storage/waterbalance/fire/out/ap23/",
master_csv = "SR_fire_points_extracted.csv"
),
"Middle Rockies", list(data_dir = "/home/steve/OneDrive/nothern-rockies-dryness",
out_dir = "/media/steve/storage/waterbalance/fire/out/",
master_csv = "MR_fire_points_extracted.csv"
),
)
settings <- settings %>% filter(region == region_name)
## Directory containing metadata file
data_dir <- settings$params[[1]]$data_dir
## Directory with csv files from merge-csv.R
out_dir <- settings$params[[1]]$out_dir
master_csv <- settings$params[[1]]$master_csvMiddle Rockies
Preprocessing
- Get MTBS database
- Extract fire polys by location of interest (Ecoregion)
- Create centroids for each poly (Pole of inaccessability in QGIS)
- Get LANDFIRE EVT (Existing Vegetation Type data)
- Use statistics tool in QGIS to get mean, median, and mode EVT by pixel count. Only mode is used going further.
- Create master csv file with points of interest. Needs columns ‘Event_ID’, ‘Incid_Type’, ‘Ig_Date’, ‘x’ (=Longitude), ‘y’ (=Latitiude), and ‘LF20_EVT_5’ (Landfire 2020 EVT type). Further analysis can look at finer resolution LANDFIRE vegetation types offered.
- Get Historical WB data for each centroid from Mike
- Get Historical gridmet climate data for Each centroid
- gridmet-historical.R - download gridmet data for each point
- pr, tmmx, tmmn, rmax, rmin, vpd
- creates one dir ‘gridmet-$var’ for each variable
- merge-csv.R
- join Mike’s data and historical gridmet data into one csv for each Event_ID
- use these csvs for input to this script
Master dataframe
Read in file containing meta data for fires of interest.
setwd(data_dir)
master <- read_csv(master_csv) %>%
filter(Incid_Type == "Wildfire") %>%
## filter(year(Ig_Date) <= 2014) %>%
rename("maj_veg_cl" = LF20_EVT_5,
"Acres" = BurnBndAc) %>%
## maj_class is the majority
## of pixels classification from Landfire 30m pixels in burned
## areas, this data was extracted in QGIS.
mutate(maj_veg_cl = case_match(maj_veg_cl,
c("Herb", "Shrub", "Sparse") ~ "non-forest",
"Tree" ~ "forest",
c("Agriculture", "Water") ~ "other"),
filename = paste0(Event_ID, ".csv"),
doy = yday(Ig_Date))## Rows: 466 Columns: 45
## ── Column specification ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────
## Delimiter: ","
## chr (18): Event_ID, irwinID, Incid_Name, Incid_Type, Map_Prog, Asmnt_Type, Pre_ID, Post_ID, Comment, LF20_EVT_2, LF20_EVT_4, LF20_EVT_...
## dbl (26): Map_ID, BurnBndAc, BurnBndLat, BurnBndLon, Perim_ID, dNBR_offst, dNBR_stdDv, NoData_T, IncGreen_T, Low_T, Mod_T, High_T, _EV...
## date (1): Ig_Date
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.Map
library(sp)
library(sf)
coords <- data.frame(
x = master$x,
y = master$y
)
spdf <- SpatialPointsDataFrame(coords = coords, data = master)
ecoregions <- st_transform(read_sf("/home/steve/OneDrive/nothern-rockies-dryness/data/us_eco_l3/us_eco_l3.shp"), crs = 4326)
ecoregion <- ecoregions %>% filter(US_L3NAME == region_name)
ggplot(ecoregion) +
geom_sf() +
geom_point(data = master, aes(x, y, color = maj_veg_cl))
summary_tbl <- master %>% group_by(maj_veg_cl) %>% summarise(n = n())
knitr::kable(summary_tbl)| maj_veg_cl | n |
|---|---|
| forest | 110 |
| non-forest | 307 |
Determine Fire Season
forest<-subset(master, maj_veg_cl=="forest")
forest_q<-quantile(forest$doy,c(0,1))#90% of data fall within this range 0.05-0.95# 0.27 to 0.75 or 0,1 to determine min max
nonforest<-subset(master, maj_veg_cl=="non-forest")
nonforest_q<-quantile(nonforest$doy,c(0,1))#90% of data fall within this range 0.05-0.95# 0.27 to 0.75 or 0,1 to determine min
#start and end dates of fire season by vegetation type derived from day of year in min_max table
doy_min_max<-as.data.frame(rbind(forest_q, nonforest_q))
forest_start<-doy_min_max[1,1]
nonforest_start<-doy_min_max[2,1]
forest_end<-doy_min_max[1,2]
nonforest_end<-doy_min_max[2,2]
knitr::kable(doy_min_max)| 0% | 100% | |
|---|---|---|
| forest_q | 93 | 289 |
| nonforest_q | 6 | 350 |
Compute Percentiles
cluster_copy(cluster, "master")
cluster_copy(cluster, "window")
cluster_copy(cluster, "vars")
cluster_copy(cluster, "forest_start")
cluster_copy(cluster, "nonforest_start")
cluster_copy(cluster, "forest_end")
cluster_copy(cluster, "nonforest_end")
cluster_library(cluster, 'tidyverse')
cluster_library(cluster, 'slider')
cluster_library(cluster, 'zoo')
files <- paste0(out_dir, master$filename)
cluster_assign_partition(cluster, files = files)
cluster_send(cluster, wbdata <- vroom::vroom(files, id = "PATH"))
wbdata <- party_df(cluster, "wbdata")
master <- party_df(cluster, "master")
wbdata <- wbdata %>%
mutate(Event_ID = str_extract(PATH, '(?<=/)[^/]*(?=\\.[^./]+$)')) %>% ## Extract file name without extension, use as Event_ID
group_by(Event_ID) %>%
mutate(RD = 100 - RH,
WHC = max(SOIL), ## Guestimate soil WHC from max soil moisture observed
AWSSM = WHC - SOIL) %>%
left_join(select(master,
Event_ID, Ig_Date, maj_veg_cl, Acres),
by = "Event_ID") %>%
mutate(fire = if_else(Date == Ig_Date, 1, 0))
start <- proc.time() # Start clock
wbdata_smoothed <- select(wbdata, c(Event_ID, Ig_Date, maj_veg_cl, Acres, fire, Date, vars)) %>%
## mutate_at(vars[1:5],
## ~ slide_dbl(.x, mean, .before = window, .after = 0, .complete = TRUE)) %>%
## mutate_at(vars[6:9],
## ~ slide_dbl(.x, sum, .before = window, .after = 0, .complete = TRUE)) %>%
mutate_at(vars[1:5],
~ rollapply(.x, window, mean, by = 1, partial = FALSE, fill = NA, align = "right")) %>%
mutate_at(vars[6:9],
~ rollapply(.x, window, sum, by = 1, partial = FALSE, fill = NA, align = "right"))
time_elapsed_parallel <- proc.time() - start # End clock
wbdata_percentiles_forest <- select(wbdata_smoothed, c(Event_ID, Ig_Date, maj_veg_cl, Acres, fire, Date, vars)) %>%
filter(maj_veg_cl == "forest") %>%
filter(Date <= Ig_Date,
yday(Date) >= forest_start, ## Filter to fire season
yday(Date) <= forest_end) %>% ## To get this to match results from original dryness_script2, we only filter beginning of season. Original script also didn't use doy, instead it filtered by absolute date, so would only chop off first pre-fireseason from first year in record.
## slice(window + 1):n()) %>%
collect() %>%
drop_na(vars) %>% ### Drop NAs from periods not included in rolling calcs
mutate(across(vars,
percent_rank))
wbdata_percentiles_nonforest <- select(wbdata_smoothed, c(Event_ID, Ig_Date, maj_veg_cl, Acres, fire, Date, vars)) %>%
filter(maj_veg_cl == "non-forest") %>%
filter(Date <= Ig_Date,
yday(Date) >= nonforest_start, ## Filter to fire season
yday(Date) <= nonforest_end) %>%
## slice((window + 1):n()) %>%
collect() %>%
drop_na(vars) %>% ### Drop NAs from periods not included in rolling calcs
mutate(across(vars,
percent_rank))
ignitions_forest <- wbdata_percentiles_forest %>% filter(fire == 1)
nonignitions_forest <- wbdata_percentiles_forest %>% filter(fire == 0)
ignitions_forest_long <- ignitions_forest %>% pivot_longer(cols = vars)
ignitions_nonforest <- wbdata_percentiles_nonforest %>% filter(fire == 1)
nonignitions_nonforest <- wbdata_percentiles_nonforest %>% filter(fire == 0)
ignitions_nonforest_long <- ignitions_nonforest %>% pivot_longer(cols = vars)
ignitions_long <- rbind(ignitions_forest_long, ignitions_nonforest_long)Corrplot
library(GGally)
ggcorr(wbdata_percentiles_forest, method = c("everything", "pearson")) ## Warning in ggcorr(wbdata_percentiles_forest, method = c("everything", "pearson")): data in column(s) 'Event_ID', 'Ig_Date', 'maj_veg_cl',
## 'Date' are not numeric and were ignored
ROC Analysis
demonstrate histogram of conditions associated with fire or no fire to see how much overlap there is. This helps interpret the ROC 
Determine contribution of various parameters to burned area 
## $gfs
## [1] 0.000000000 0.008428674 0.004943795 0.009037956 0.004943203 0.012320908 0.006039262 0.009198240 0.013260896 0.009231091 0.009292778
## [12] 0.009230613 0.013254085 0.008442127 0.011956989 0.020454097 0.009747693 0.004979038 0.019080426 0.006395096 0.010547810 0.013386536
## [23] 0.009747171 0.013549140 0.009112826 0.012353812 0.017489344 0.019078476 0.006394807 0.010547101 0.013386258 0.014017935 0.018668797
## [34] 0.016545246 0.010583923 0.015734132 0.014343938 0.010054835 0.009338704 0.020658692 0.010412629 0.012266820 0.020479946 0.010054198
## [45] 0.014027604 0.009329964 0.012673290 0.020545331 0.020655766 0.010411399 0.012266406 0.020479751 0.014155244 0.018756831 0.030039403
## [56] 0.012082734 0.022160951 0.020949226 0.009894121 0.020131587 0.011323724 0.012616083 0.017535476 0.019361126 0.006448104 0.010708686
## [67] 0.013611203 0.021613319 0.021675973 0.019184571 0.010816431 0.019649565 0.014413966 0.020129027 0.011322071 0.012615671 0.017535294
## [78] 0.014021254 0.018686697 0.024926939 0.012611187 0.019386597 0.017499481 0.021609197 0.021674241 0.019182810 0.010816020 0.019646350
## [89] 0.014413737 0.018708604 0.017629664 0.019033049 0.016341809 0.010259390 0.021222438 0.011436684 0.012967547 0.020572624 0.020943174
## [100] 0.010553592 0.012470552 0.021026615 0.021629980 0.025101425 0.035528834 0.013407063 0.024436435 0.020989151 0.021219023 0.011434941
## [111] 0.012967046 0.020572413 0.014773448 0.019338807 0.030047423 0.012850581 0.022272838 0.020999244 0.021625901 0.025097880 0.035523916
## [122] 0.013405719 0.024433321 0.020988904 0.019465942 0.030046671 0.030360725 0.022966225 0.020380675 0.011553898 0.012855068 0.017797168
## [133] 0.021649623 0.024267186 0.029480505 0.014147585 0.022105989 0.017548414 0.021947731 0.022095885 0.019482490 0.010980982 0.020160959
## [144] 0.014685344 0.025652919 0.026663444 0.022077916 0.019696490 0.021645372 0.024264181 0.029477713 0.014145895 0.022103265 0.017548228
## [155] 0.019013049 0.025049739 0.024947164 0.019477821 0.025647864 0.026657076 0.022075692 0.019693366 0.019415219 0.021627466 0.011706235
## [166] 0.013286805 0.021179235 0.022347037 0.025529738 0.035529209 0.014290651 0.024598403 0.021039928 0.021960405 0.025588364 0.036624750
## [177] 0.013654401 0.025236083 0.021552970 0.026014490 0.039554919 0.035954226 0.026015960 0.022342119 0.025525816 0.035524283 0.014288892
## [188] 0.024595085 0.021039669 0.019930636 0.030055610 0.030361450 0.023022032 0.026009597 0.039545358 0.035949213 0.026012075 0.030395448
## [199] 0.022003031 0.024670601 0.029902842 0.014488943 0.022615908 0.017805635 0.025728881 0.034294572 0.029494286 0.022588360 0.026205696
## [210] 0.027333279 0.022489038 0.020212982 0.027463063 0.025723992 0.034287804 0.029491473 0.022585424 0.025055463 0.027456809 0.022830305
## [221] 0.026149598 0.036642521 0.014687154 0.025497247 0.021650017 0.026573941 0.039568726 0.035957145 0.026090639 0.026557954 0.041007438
## [232] 0.037070688 0.026901219 0.040989714 0.026568355 0.039559058 0.035952156 0.026086612 0.030396610 0.040979071 0.026259504 0.035034083
## [243] 0.029924474 0.023079701 0.034514696 0.028172497 0.034507823 0.027284440 0.041072207 0.037074228 0.027048210 0.040989751 0.042548601
## [254] 0.040979113 0.035237029 0.042565817
##
## $IJ
## I J Total
## T 0.005323182 0.0031054921 0.008428674
## SOIL 0.002472641 0.0024711536 0.004943795
## GDD 0.002732331 0.0063056251 0.009037956
## AWSSM 0.002468716 0.0024744872 0.004943203
## AET 0.011733675 0.0005872333 0.012320908
## D 0.003478901 0.0025603615 0.006039262
## RD 0.003224324 0.0059739159 0.009198240
## VPD 0.011132047 0.0021288487 0.013260896
##
## $I.perc
## ind.exp.var
## T 12.505767
## SOIL 5.808983
## GDD 6.419074
## AWSSM 5.799761
## AET 27.565958
## D 8.172992
## RD 7.574914
## VPD 26.152551
##
## $params
## $params$full.model
## [1] "y ~ T + SOIL + GDD + AWSSM + AET + D + RD + VPD"
##
## $params$family
## [1] "gaussian"
##
## $params$link
## [1] "default"
##
## $params$gof
## [1] "Rsqu"
## $gfs
## [1] 0.0000000000 0.0091767168 0.0070542364 0.0091581213 0.0070522670 0.0002228995 0.0121621729 0.0003373826 0.0081576380 0.0130156440
## [11] 0.0092065236 0.0130156049 0.0101874452 0.0127901176 0.0091903916 0.0092485324 0.0130020149 0.0071654190 0.0123426963 0.0121665023
## [21] 0.0073312275 0.0107874688 0.0130019850 0.0101697661 0.0127860046 0.0091724422 0.0092332066 0.0123404428 0.0121665375 0.0073285374
## [31] 0.0107870527 0.0131092126 0.0008080025 0.0081914503 0.0134470190 0.0121650857 0.0091220143 0.0130390982 0.0130157032 0.0136000847
## [41] 0.0132444995 0.0138091654 0.0134948061 0.0130390160 0.0102049249 0.0128019360 0.0092174288 0.0092812116 0.0135996054 0.0132447291
## [51] 0.0138084337 0.0134947844 0.0131107808 0.0106468251 0.0102459130 0.0138448664 0.0137306622 0.0094382881 0.0130020464 0.0135892807
## [61] 0.0132361814 0.0138009029 0.0134746066 0.0124015588 0.0121933513 0.0075120510 0.0107985602 0.0140595731 0.0124077624 0.0128344276
## [71] 0.0134793313 0.0121765385 0.0128143310 0.0135888071 0.0132364165 0.0138001805 0.0134746037 0.0131104506 0.0106357773 0.0102260960
## [81] 0.0138444400 0.0137237894 0.0094299519 0.0140598140 0.0124050525 0.0128332035 0.0134793508 0.0121765999 0.0128129534 0.0139150154
## [91] 0.0135222256 0.0097382887 0.0135704964 0.0130394031 0.0136370264 0.0132617064 0.0138154157 0.0135107961 0.0136053313 0.0132516631
## [101] 0.0138288109 0.0134948142 0.0140940846 0.0140574337 0.0138815591 0.0143181791 0.0141334125 0.0138296337 0.0136364605 0.0132619109
## [111] 0.0138146639 0.0135107323 0.0131256132 0.0106504896 0.0102600309 0.0138449331 0.0137332959 0.0094603474 0.0140943041 0.0140564004
## [121] 0.0138811557 0.0143181058 0.0141338395 0.0138289916 0.0139373167 0.0139384111 0.0107003985 0.0140194723 0.0135943865 0.0132437619
## [131] 0.0138200144 0.0134746071 0.0140918029 0.0140486949 0.0138636122 0.0143160110 0.0141220330 0.0138195494 0.0140628535 0.0124919025
## [141] 0.0128603073 0.0134804547 0.0122000564 0.0128979481 0.0146723967 0.0142924084 0.0135276542 0.0137662535 0.0140920252 0.0140476697
## [151] 0.0138632241 0.0143159426 0.0141224685 0.0138189140 0.0139370323 0.0139333156 0.0106922053 0.0140179382 0.0146723242 0.0142928029
## [161] 0.0135258089 0.0137661572 0.0139512689 0.0136443761 0.0132673572 0.0138362872 0.0135109159 0.0141224025 0.0140720682 0.0139087581
## [171] 0.0143190933 0.0141388970 0.0138363790 0.0140966785 0.0140828370 0.0138849234 0.0143188407 0.0141612687 0.0138444009 0.0147448467
## [181] 0.0147603972 0.0140752924 0.0144585875 0.0141225675 0.0140709896 0.0139082760 0.0143190145 0.0141393106 0.0138357164 0.0139378810
## [191] 0.0139426534 0.0107037081 0.0140195597 0.0147447228 0.0147609066 0.0140743654 0.0144587162 0.0141300227 0.0140944692 0.0140736689
## [201] 0.0138667178 0.0143165870 0.0141513704 0.0138340120 0.0147432530 0.0147517192 0.0140645492 0.0144545427 0.0146727177 0.0143028398
## [211] 0.0135987032 0.0137694392 0.0146790988 0.0147431324 0.0147522373 0.0140636260 0.0144546766 0.0141283427 0.0146790797 0.0141238746
## [221] 0.0140997391 0.0139136109 0.0143198613 0.0141649890 0.0138521903 0.0147507696 0.0147732834 0.0140905809 0.0144595194 0.0147457472
## [231] 0.0147766347 0.0140954989 0.0144610506 0.0149077266 0.0147506252 0.0147737585 0.0140896092 0.0144596430 0.0141307246 0.0149079193
## [241] 0.0147441079 0.0147685943 0.0140846136 0.0144572277 0.0149031401 0.0146791161 0.0149033363 0.0147519986 0.0147873270 0.0141127634
## [251] 0.0144617987 0.0149138943 0.0149103339 0.0149140682 0.0149058523 0.0149160142
##
## $IJ
## I J Total
## T 0.0021394409 0.0070372759 0.0091767168
## SOIL 0.0020189542 0.0050352823 0.0070542364
## GDD 0.0021297129 0.0070284084 0.0091581213
## AWSSM 0.0020184337 0.0050338333 0.0070522670
## AET 0.0008361529 -0.0006132534 0.0002228995
## D 0.0034993839 0.0086627890 0.0121621729
## RD 0.0004550989 -0.0001177163 0.0003373826
## VPD 0.0018188367 0.0063388013 0.0081576380
##
## $I.perc
## ind.exp.var
## T 14.343248
## SOIL 13.535480
## GDD 14.278030
## AWSSM 13.531991
## AET 5.605740
## D 23.460583
## RD 3.051076
## VPD 12.193852
##
## $params
## $params$full.model
## [1] "y ~ T + SOIL + GDD + AWSSM + AET + D + RD + VPD"
##
## $params$family
## [1] "gaussian"
##
## $params$link
## [1] "default"
##
## $params$gof
## [1] "Rsqu"Plots of wb percentile vs. area burned show weak relationships 
ECDF
Fire Danger Model
First Order
CWD Forest
ecdf_fn <- ecdf(ignitions_forest$D)
summary(ecdf_fn)## Empirical CDF: 110 unique values with summary
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.2235 0.9043 0.9532 0.9098 0.9821 1.0000estimated <- data.frame(
dec = seq(from = 0.01, to = 0.99, by = 0.01),
quants = quantile(ecdf_fn, dec)
)
estimated$proportions <- ecdf_fn(estimated$quants)
lm <- lm(log(proportions) ~ quants, data = estimated)
estimated$predicted <- exp(predict(lm, data = estimated$quants))
##plot(lm)
summary(lm)##
## Call:
## lm(formula = log(proportions) ~ quants, data = estimated)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.5362 -0.2127 0.0214 0.1946 0.9200
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -7.9872 0.2178 -36.67 <2e-16 ***
## quants 7.6748 0.2362 32.49 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2562 on 97 degrees of freedom
## Multiple R-squared: 0.9158, Adjusted R-squared: 0.915
## F-statistic: 1056 on 1 and 97 DF, p-value: < 2.2e-16a <- exp(lm$coefficients[[1]])
b <- lm$coefficients[[2]]
plot(ecdf_fn)
points(y = estimated$proportions, x = estimated$quants)
points(y = estimated$predicted, x = estimated$quants, col = "red")
curve(a*exp(b*x), add = TRUE, col = "blue")
Estimated function predicting proportion of historic fires burning in forests at or below given percentile of climatic water deficit on day of ignition:
\(\begin{align*}f(x) = 3.3978819\times 10^{-4} e^{7.6748305 x}\end{align*}\)
VPD Forest
ecdf_fn <- ecdf(ignitions_forest$VPD)
summary(ecdf_fn)## Empirical CDF: 109 unique values with summary
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.1177 0.8893 0.9429 0.9076 0.9845 1.0000estimated <- data.frame(
dec = seq(from = 0.01, to = 0.99, by = 0.01),
quants = quantile(ecdf_fn, dec)
)
estimated$proportions <- ecdf_fn(estimated$quants)
lm <- lm(log(proportions) ~ quants, data = estimated)
estimated$predicted <- exp(predict(lm, data = estimated$quants))
##plot(lm)
summary(lm)##
## Call:
## lm(formula = log(proportions) ~ quants, data = estimated)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.48860 -0.15060 0.05087 0.14196 0.68756
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -8.5004 0.1972 -43.11 <2e-16 ***
## quants 8.2399 0.2141 38.48 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.219 on 97 degrees of freedom
## Multiple R-squared: 0.9385, Adjusted R-squared: 0.9379
## F-statistic: 1481 on 1 and 97 DF, p-value: < 2.2e-16a <- exp(lm$coefficients[[1]])
b <- lm$coefficients[[2]]
plot(ecdf_fn)
points(y = estimated$proportions, x = estimated$quants)
points(y = estimated$predicted, x = estimated$quants, col = "red")
curve(a*exp(b*x), add = TRUE, col = "blue")
Estimated function predicting proportion of historic fires burning in forests at or below given percentile of vapor pressure deficit on day of ignition:
\(\begin{align*}f(x) = 2.0338181\times 10^{-4} e^{8.239898 x}\end{align*}\)
CWD Non-Forest
ecdf_fn <- ecdf(ignitions_nonforest$D)
summary(ecdf_fn)## Empirical CDF: 300 unique values with summary
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.9108 0.9628 0.9100 0.9856 1.0000estimated <- data.frame(
dec = seq(from = 0.01, to = 0.99, by = 0.01),
quants = quantile(ecdf_fn, dec)
)
estimated$proportions <- ecdf_fn(estimated$quants)
lm <- lm(log(proportions) ~ quants, data = estimated)
estimated$predicted <- exp(predict(lm, data = estimated$quants))
##plot(lm)
summary(lm)##
## Call:
## lm(formula = log(proportions) ~ quants, data = estimated)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.45038 -0.24811 0.03398 0.21852 0.78489
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -7.3150 0.2020 -36.21 <2e-16 ***
## quants 6.9213 0.2184 31.69 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2702 on 97 degrees of freedom
## Multiple R-squared: 0.9119, Adjusted R-squared: 0.911
## F-statistic: 1004 on 1 and 97 DF, p-value: < 2.2e-16a <- exp(lm$coefficients[[1]])
b <- lm$coefficients[[2]]
plot(ecdf_fn)
points(y = estimated$proportions, x = estimated$quants)
points(y = estimated$predicted, x = estimated$quants, col = "red")
curve(a*exp(b*x), add = TRUE, col = "blue")
Estimated function predicting proportion of historic fires burning in nonforests at or below given percentile of climatic water deficit on day of ignition:
\(\begin{align*}f(x) = 6.6550577\times 10^{-4} e^{6.9212675 x}\end{align*}\)
VPD Non-Forest
ecdf_fn <- ecdf(ignitions_nonforest$VPD)
summary(ecdf_fn)## Empirical CDF: 304 unique values with summary
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.1948 0.8894 0.9602 0.9016 0.9865 1.0000estimated <- data.frame(
dec = seq(from = 0.01, to = 0.99, by = 0.01),
quants = quantile(ecdf_fn, dec)
)
estimated$proportions <- ecdf_fn(estimated$quants)
lm <- lm(log(proportions) ~ quants, data = estimated)
estimated$predicted <- exp(predict(lm, data = estimated$quants))
##plot(lm)
summary(lm)##
## Call:
## lm(formula = log(proportions) ~ quants, data = estimated)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.37214 -0.17621 -0.00755 0.19867 0.54198
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6.8199 0.1523 -44.78 <2e-16 ***
## quants 6.4541 0.1662 38.84 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2238 on 97 degrees of freedom
## Multiple R-squared: 0.9396, Adjusted R-squared: 0.9389
## F-statistic: 1508 on 1 and 97 DF, p-value: < 2.2e-16a <- exp(lm$coefficients[[1]])
b <- lm$coefficients[[2]]
plot(ecdf_fn)
points(y = estimated$proportions, x = estimated$quants)
points(y = estimated$predicted, x = estimated$quants, col = "red")
curve(a*exp(b*x), add = TRUE, col = "blue")
Estimated function predicting proportion of historic fires burning in nonforests at or below given percentile of vapor pressure deficit on day of ignition:
\(\begin{align*}f(x) = 0.0010918 e^{6.45408 x}\end{align*}\)
Second Order
Fits to ECDF curves may be improved by including a quadtratic term in regression
CWD Forest
ecdf_fn <- ecdf(ignitions_forest$D)
summary(ecdf_fn)## Empirical CDF: 110 unique values with summary
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.2235 0.9043 0.9532 0.9098 0.9821 1.0000estimated <- data.frame(
dec = seq(from = 0.01, to = 0.99, by = 0.01),
quants = quantile(ecdf_fn, dec)
)
estimated$proportions <- ecdf_fn(estimated$quants)
lm <- lm(log(proportions) ~ poly(quants, 2, raw = TRUE), data = estimated)
estimated$predicted <- exp(predict(lm, data = estimated$quants))
##plot(lm)
summary(lm)##
## Call:
## lm(formula = log(proportions) ~ poly(quants, 2, raw = TRUE),
## data = estimated)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.38889 -0.04976 0.01981 0.04779 0.28434
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.1356 0.2590 -4.384 2.98e-05 ***
## poly(quants, 2, raw = TRUE)1 -11.0575 0.6839 -16.168 < 2e-16 ***
## poly(quants, 2, raw = TRUE)2 12.1149 0.4393 27.576 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.08621 on 96 degrees of freedom
## Multiple R-squared: 0.9906, Adjusted R-squared: 0.9904
## F-statistic: 5040 on 2 and 96 DF, p-value: < 2.2e-16b0 <- lm$coefficients[[1]]
b1 <- lm$coefficients[[2]]
b2 <- lm$coefficients[[3]]
plot(ecdf_fn)
points(y = estimated$proportions, x = estimated$quants, col = "blue")
points(y = estimated$predicted, x = estimated$quants, col = "red")
curve(exp(b0 + b1*x + b2*(x^2)), add = TRUE, col = "red")
##curve(exp(predict(lm, newdata = data.frame(quants = x))), add = TRUE, col = "red")Estimated function predicting proportion of historic fires burning in forests at or below given percentile of CWD on day of ignition:
\(\begin{align*}f(x) = e^{-1.1356493 + -11.0574732x + 12.1148891x^2}\end{align*}\)
VPD Forest
ecdf_fn <- ecdf(ignitions_forest$VPD)
summary(ecdf_fn)## Empirical CDF: 109 unique values with summary
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.1177 0.8893 0.9429 0.9076 0.9845 1.0000estimated <- data.frame(
dec = seq(from = 0.01, to = 0.99, by = 0.01),
quants = quantile(ecdf_fn, dec)
)
estimated$proportions <- ecdf_fn(estimated$quants)
lm <- lm(log(proportions) ~ poly(quants, 2, raw = TRUE), data = estimated)
estimated$predicted <- exp(predict(lm, data = estimated$quants))
##plot(lm)
summary(lm)##
## Call:
## lm(formula = log(proportions) ~ poly(quants, 2, raw = TRUE),
## data = estimated)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.35282 -0.02766 0.00024 0.02699 0.28947
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.1035 0.2906 -3.797 0.000256 ***
## poly(quants, 2, raw = TRUE)1 -11.1678 0.7443 -15.005 < 2e-16 ***
## poly(quants, 2, raw = TRUE)2 12.2214 0.4663 26.211 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.07707 on 96 degrees of freedom
## Multiple R-squared: 0.9925, Adjusted R-squared: 0.9923
## F-statistic: 6319 on 2 and 96 DF, p-value: < 2.2e-16b0 <- lm$coefficients[[1]]
b1 <- lm$coefficients[[2]]
b2 <- lm$coefficients[[3]]
plot(ecdf_fn)
points(y = estimated$proportions, x = estimated$quants, col = "blue")
points(y = estimated$predicted, x = estimated$quants, col = "red")
curve(exp(b0 + b1*x + b2*(x^2)), add = TRUE, col = "red")
##curve(exp(predict(lm, newdata = data.frame(quants = x))), add = TRUE, col = "red")Estimated function predicting proportion of historic fires burning in forests at or below given percentile of VPD on day of ignition:
\(\begin{align*}f(x) = e^{-1.1035326 + -11.1678437x + 12.2213889x^2}\end{align*}\)
CWD Non-forest
ecdf_fn <- ecdf(ignitions_nonforest$D)
summary(ecdf_fn)## Empirical CDF: 300 unique values with summary
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.9108 0.9628 0.9100 0.9856 1.0000estimated <- data.frame(
dec = seq(from = 0.01, to = 0.99, by = 0.01),
quants = quantile(ecdf_fn, dec)
)
estimated$proportions <- ecdf_fn(estimated$quants)
lm <- lm(log(proportions) ~ poly(quants, 2, raw = TRUE), data = estimated)
estimated$predicted <- exp(predict(lm, data = estimated$quants))
##plot(lm)
summary(lm)##
## Call:
## lm(formula = log(proportions) ~ poly(quants, 2, raw = TRUE),
## data = estimated)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.53455 -0.12165 0.00516 0.11731 0.38825
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.634 0.383 -6.877 6.19e-10 ***
## poly(quants, 2, raw = TRUE)1 -6.535 1.051 -6.217 1.31e-08 ***
## poly(quants, 2, raw = TRUE)2 8.944 0.693 12.905 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1642 on 96 degrees of freedom
## Multiple R-squared: 0.9678, Adjusted R-squared: 0.9671
## F-statistic: 1443 on 2 and 96 DF, p-value: < 2.2e-16b0 <- lm$coefficients[[1]]
b1 <- lm$coefficients[[2]]
b2 <- lm$coefficients[[3]]
plot(ecdf_fn)
points(y = estimated$proportions, x = estimated$quants, col = "blue")
points(y = estimated$predicted, x = estimated$quants, col = "red")
curve(exp(b0 + b1*x + b2*(x^2)), add = TRUE, col = "red")
##curve(exp(predict(lm, newdata = data.frame(quants = x))), add = TRUE, col = "red")Estimated function predicting proportion of historic fires burning in nonforests at or below given percentile of CWD on day of ignition:
\(\begin{align*}f(x) = e^{-2.633803 + -6.5346555x + 8.9434708x^2}\end{align*}\)
VPD Non-forest
ecdf_fn <- ecdf(ignitions_nonforest$VPD)
summary(ecdf_fn)## Empirical CDF: 304 unique values with summary
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.1948 0.8894 0.9602 0.9016 0.9865 1.0000estimated <- data.frame(
dec = seq(from = 0.01, to = 0.99, by = 0.01),
quants = quantile(ecdf_fn, dec)
)
estimated$proportions <- ecdf_fn(estimated$quants)
lm <- lm(log(proportions) ~ poly(quants, 2, raw = TRUE), data = estimated)
estimated$predicted <- exp(predict(lm, data = estimated$quants))
##plot(lm)
summary(lm)##
## Call:
## lm(formula = log(proportions) ~ poly(quants, 2, raw = TRUE),
## data = estimated)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.46531 -0.10737 -0.03828 0.10860 0.33554
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.4958 0.2975 -11.750 < 2e-16 ***
## poly(quants, 2, raw = TRUE)1 -3.2051 0.8235 -3.892 0.000184 ***
## poly(quants, 2, raw = TRUE)2 6.4666 0.5467 11.828 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1435 on 96 degrees of freedom
## Multiple R-squared: 0.9754, Adjusted R-squared: 0.9749
## F-statistic: 1904 on 2 and 96 DF, p-value: < 2.2e-16b0 <- lm$coefficients[[1]]
b1 <- lm$coefficients[[2]]
b2 <- lm$coefficients[[3]]
plot(ecdf_fn)
points(y = estimated$proportions, x = estimated$quants, col = "blue")
points(y = estimated$predicted, x = estimated$quants, col = "red")
curve(exp(b0 + b1*x + b2*(x^2)), add = TRUE, col = "red")
##curve(exp(predict(lm, newdata = data.frame(quants = x))), add = TRUE, col = "red")Estimated function predicting proportion of historic fires burning in nonforests at or below given percentile of VPD on day of ignition:
\(\begin{align*}f(x) = e^{-3.4958205 + -3.2051301x + 6.4665794x^2}\end{align*}\)
Third Order
CWD Forest
ecdf_fn <- ecdf(ignitions_forest$D)
summary(ecdf_fn)## Empirical CDF: 110 unique values with summary
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.2235 0.9043 0.9532 0.9098 0.9821 1.0000estimated <- data.frame(
dec = seq(from = 0.01, to = 0.99, by = 0.01),
quants = quantile(ecdf_fn, dec)
)
estimated$proportions <- ecdf_fn(estimated$quants)
lm <- lm(log(proportions) ~ poly(quants, 3, raw = TRUE), data = estimated)
estimated$predicted <- exp(predict(lm, data = estimated$quants))
##plot(lm)
summary(lm)##
## Call:
## lm(formula = log(proportions) ~ poly(quants, 3, raw = TRUE),
## data = estimated)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.153256 -0.015155 0.000106 0.017564 0.104976
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -10.3754 0.4445 -23.34 <2e-16 ***
## poly(quants, 3, raw = TRUE)1 29.9760 1.9361 15.48 <2e-16 ***
## poly(quants, 3, raw = TRUE)2 -45.2506 2.6835 -16.86 <2e-16 ***
## poly(quants, 3, raw = TRUE)3 25.6546 1.1973 21.43 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03588 on 95 degrees of freedom
## Multiple R-squared: 0.9984, Adjusted R-squared: 0.9983
## F-statistic: 1.955e+04 on 3 and 95 DF, p-value: < 2.2e-16b0 <- lm$coefficients[[1]]
b1 <- lm$coefficients[[2]]
b2 <- lm$coefficients[[3]]
b3 <- lm$coefficients[[4]]
plot(ecdf_fn)
points(y = estimated$proportions, x = estimated$quants, col = "blue")
points(y = estimated$predicted, x = estimated$quants, col = "red")
curve(exp(b0 + b1*x + b2*(x^2) + b3*(x^3)), add = TRUE, col = "red")
##curve(exp(predict(lm, newdata = data.frame(quants = x))), add = TRUE, col = "red")Estimated function predicting proportion of historic fires burning in forests at or below given percentile of CWD on day of ignition:
\(\begin{align*}f(x) = e^{-10.3754488 + 29.9760209x + -45.2506376x^2 + 25.6546155x^3}\end{align*}\)
VPD Forest
ecdf_fn <- ecdf(ignitions_forest$VPD)
summary(ecdf_fn)## Empirical CDF: 109 unique values with summary
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.1177 0.8893 0.9429 0.9076 0.9845 1.0000estimated <- data.frame(
dec = seq(from = 0.01, to = 0.99, by = 0.01),
quants = quantile(ecdf_fn, dec)
)
estimated$proportions <- ecdf_fn(estimated$quants)
lm <- lm(log(proportions) ~ poly(quants, 3, raw = TRUE), data = estimated)
estimated$predicted <- exp(predict(lm, data = estimated$quants))
##plot(lm)
summary(lm)##
## Call:
## lm(formula = log(proportions) ~ poly(quants, 3, raw = TRUE),
## data = estimated)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.19886 -0.02833 -0.00305 0.03567 0.15165
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -9.737 1.174 -8.297 7.10e-13 ***
## poly(quants, 3, raw = TRUE)1 24.581 4.801 5.120 1.59e-06 ***
## poly(quants, 3, raw = TRUE)2 -35.252 6.337 -5.563 2.44e-07 ***
## poly(quants, 3, raw = TRUE)3 20.405 2.719 7.504 3.29e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.06139 on 95 degrees of freedom
## Multiple R-squared: 0.9953, Adjusted R-squared: 0.9951
## F-statistic: 6659 on 3 and 95 DF, p-value: < 2.2e-16b0 <- lm$coefficients[[1]]
b1 <- lm$coefficients[[2]]
b2 <- lm$coefficients[[3]]
b3 <- lm$coefficients[[4]]
plot(ecdf_fn)
points(y = estimated$proportions, x = estimated$quants, col = "blue")
points(y = estimated$predicted, x = estimated$quants, col = "red")
curve(exp(b0 + b1*x + b2*(x^2) + b3*(x^3)), add = TRUE, col = "red")
##curve(exp(predict(lm, newdata = data.frame(quants = x))), add = TRUE, col = "red")Estimated function predicting proportion of historic fires burning in forests at or below given percentile of VPD on day of ignition:
\(\begin{align*}f(x) = e^{-9.7370222 + 24.5810007x + -35.2519182x^2 + 20.4045319x^3}\end{align*}\)
CWD Non-forest
ecdf_fn <- ecdf(ignitions_nonforest$D)
summary(ecdf_fn)## Empirical CDF: 300 unique values with summary
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.9108 0.9628 0.9100 0.9856 1.0000estimated <- data.frame(
dec = seq(from = 0.01, to = 0.99, by = 0.01),
quants = quantile(ecdf_fn, dec)
)
estimated$proportions <- ecdf_fn(estimated$quants)
lm <- lm(log(proportions) ~ poly(quants, 3, raw = TRUE), data = estimated)
estimated$predicted <- exp(predict(lm, data = estimated$quants))
##plot(lm)
summary(lm)##
## Call:
## lm(formula = log(proportions) ~ poly(quants, 3, raw = TRUE),
## data = estimated)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.26542 -0.04641 0.01102 0.04011 0.18427
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -11.7509 0.5324 -22.07 <2e-16 ***
## poly(quants, 3, raw = TRUE)1 37.9079 2.4902 15.22 <2e-16 ***
## poly(quants, 3, raw = TRUE)2 -57.5114 3.6630 -15.70 <2e-16 ***
## poly(quants, 3, raw = TRUE)3 31.2572 1.7159 18.22 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.07788 on 95 degrees of freedom
## Multiple R-squared: 0.9928, Adjusted R-squared: 0.9926
## F-statistic: 4386 on 3 and 95 DF, p-value: < 2.2e-16b0 <- lm$coefficients[[1]]
b1 <- lm$coefficients[[2]]
b2 <- lm$coefficients[[3]]
b3 <- lm$coefficients[[4]]
plot(ecdf_fn)
points(y = estimated$proportions, x = estimated$quants, col = "blue")
points(y = estimated$predicted, x = estimated$quants, col = "red")
curve(exp(b0 + b1*x + b2*(x^2) + b3*(x^3)), add = TRUE, col = "red")
##curve(exp(predict(lm, newdata = data.frame(quants = x))), add = TRUE, col = "red")Estimated function predicting proportion of historic fires burning in nonforests at or below given percentile of CWD on day of ignition:
\(\begin{align*}f(x) = e^{-11.7509331 + 37.9078727x + -57.5114143x^2 + 31.2571648x^3}\end{align*}\)
VPD Non-forest
ecdf_fn <- ecdf(ignitions_nonforest$VPD)
summary(ecdf_fn)## Empirical CDF: 304 unique values with summary
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.1948 0.8894 0.9602 0.9016 0.9865 1.0000estimated <- data.frame(
dec = seq(from = 0.01, to = 0.99, by = 0.01),
quants = quantile(ecdf_fn, dec)
)
estimated$proportions <- ecdf_fn(estimated$quants)
lm <- lm(log(proportions) ~ poly(quants, 3, raw = TRUE), data = estimated)
estimated$predicted <- exp(predict(lm, data = estimated$quants))
##plot(lm)
summary(lm)##
## Call:
## lm(formula = log(proportions) ~ poly(quants, 3, raw = TRUE),
## data = estimated)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.24417 -0.05638 -0.01129 0.05733 0.19818
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -10.2383 0.4256 -24.06 <2e-16 ***
## poly(quants, 3, raw = TRUE)1 30.1834 2.0164 14.97 <2e-16 ***
## poly(quants, 3, raw = TRUE)2 -44.0238 2.9971 -14.69 <2e-16 ***
## poly(quants, 3, raw = TRUE)3 23.9467 1.4155 16.92 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.072 on 95 degrees of freedom
## Multiple R-squared: 0.9939, Adjusted R-squared: 0.9937
## F-statistic: 5135 on 3 and 95 DF, p-value: < 2.2e-16b0 <- lm$coefficients[[1]]
b1 <- lm$coefficients[[2]]
b2 <- lm$coefficients[[3]]
b3 <- lm$coefficients[[4]]
plot(ecdf_fn)
points(y = estimated$proportions, x = estimated$quants, col = "blue")
points(y = estimated$predicted, x = estimated$quants, col = "red")
curve(exp(b0 + b1*x + b2*(x^2) + b3*(x^3)), add = TRUE, col = "red")
##curve(exp(predict(lm, newdata = data.frame(quants = x))), add = TRUE, col = "red")Estimated function predicting proportion of historic fires burning in nonforests at or below given percentile of VPD on day of ignition:
\(\begin{align*}f(x) = e^{-10.2383034 + 30.1833941x + -44.0237879x^2 + 23.9466752x^3}\end{align*}\)