Settings

library(tidyverse)
library(multidplyr)
library(zoo)
library(slider)
library(pROC)

plot_theme <- theme_bw() 

region_name <- "Southern 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_csv

Southern Rockies

Preprocessing

  1. Get MTBS database
  • Extract fire polys by location of interest (Ecoregion)
  • Create centroids for each poly (Pole of inaccessability in QGIS)
  1. Get LANDFIRE EVT (Expected Vegetation Type data)
  • Use statistics tool in QGIS to get mean, median, and mode EVT by pixel count. Only mode is used going further.
  1. Create master csv file with points of interest. Needs columns ‘Event_ID’, ‘Incid_Type’, ‘Ig_Date’, ‘Lon’, ‘Lat, and ’LF20_EVT_5’ (Landfire 2020 EVT type). Further analysis can look at finer resolution LANDFIRE vegetation types offered.
  2. Get Historical WB data for each centroid from Mike
  3. 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: 290 Columns: 46
## ── 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  (27): 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 153
non-forest 86

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 79 303
nonforest_q 7 301

Compute Percentiles

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")

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
    slice(window + 1:n()) %>%
    mutate(across(vars,
                  percent_rank)) %>%
    collect()

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()) %>%
    mutate(across(vars,
                  percent_rank)) %>%
    collect()

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.004210840 0.004102275 0.003901830 0.004273405 0.017141572 0.017138879 0.025857621 0.016519833 0.004859043 0.004698204
##  [12] 0.004967783 0.017141653 0.036637941 0.025904217 0.056576727 0.004693648 0.004692941 0.019880844 0.023289764 0.025862543 0.018157906
##  [23] 0.004811514 0.017149606 0.035760432 0.025897104 0.052838976 0.019894276 0.022851123 0.025857680 0.017968712 0.024078136 0.041995795
##  [34] 0.022966519 0.027991148 0.017307693 0.027213390 0.005324944 0.005517806 0.023948541 0.048922257 0.025912573 0.056843853 0.005430329
##  [45] 0.017671483 0.036644339 0.025927587 0.056576758 0.023767706 0.048616766 0.025933750 0.056807574 0.053975735 0.050131940 0.071726369
##  [56] 0.041193947 0.058709551 0.058124733 0.005347027 0.024308142 0.047195718 0.025902857 0.053190178 0.019894705 0.023825802 0.026794941
##  [67] 0.018877239 0.025835084 0.042180009 0.023057883 0.032030029 0.023668930 0.028380189 0.024115752 0.046890788 0.025921474 0.053146292
##  [78] 0.052418254 0.049766449 0.067069226 0.040033405 0.054501513 0.053916794 0.025684806 0.042162751 0.023038015 0.031619309 0.023199783
##  [89] 0.028202868 0.042415461 0.024112772 0.043946342 0.028235801 0.005964561 0.024413692 0.048941163 0.025935702 0.056843897 0.024025587
## [100] 0.048972379 0.026829323 0.057009619 0.056971134 0.055486206 0.073381065 0.050838489 0.062696682 0.058360560 0.024228335 0.048633943
## [111] 0.025956771 0.056807624 0.053980743 0.050252164 0.071726887 0.041251629 0.058716931 0.058172424 0.057113815 0.055284129 0.073229756
## [122] 0.050489030 0.062537215 0.058337282 0.064228326 0.075050684 0.071731464 0.060367647 0.024395260 0.047256962 0.026820464 0.053372103
## [133] 0.055035120 0.054618150 0.068343621 0.048999143 0.058311844 0.054243225 0.026008864 0.042274929 0.023201137 0.032817106 0.024145578
## [144] 0.029311403 0.042487326 0.026278423 0.044698083 0.033816017 0.055157978 0.054439274 0.068226596 0.048653102 0.058141950 0.054211692
## [155] 0.061977314 0.069687642 0.067147674 0.055604247 0.042503279 0.026114523 0.044702084 0.033320974 0.045891467 0.024497779 0.048991702
## [166] 0.026859899 0.057009622 0.056973561 0.055619817 0.073381568 0.050949128 0.062733585 0.058407799 0.057174979 0.055641173 0.073753247
## [177] 0.050961017 0.062754816 0.058434193 0.064492044 0.075051111 0.073446411 0.064316185 0.057116147 0.055419577 0.073230194 0.050593934
## [188] 0.062572470 0.058385085 0.064409771 0.075058989 0.071731569 0.060384838 0.064550474 0.075053760 0.073291848 0.064222340 0.075056163
## [199] 0.055209555 0.054753965 0.068634321 0.049138716 0.058387720 0.054340621 0.062169093 0.069687686 0.068549724 0.059336768 0.042545967
## [210] 0.026412132 0.044702300 0.034469038 0.045919209 0.062214476 0.069693032 0.068427176 0.059219205 0.069787043 0.045915922 0.057176808
## [221] 0.055770658 0.073754197 0.051076283 0.062792353 0.058477793 0.064677832 0.075059219 0.073446976 0.064316655 0.064704380 0.075418966
## [232] 0.073810971 0.064325755 0.075057447 0.064735808 0.075062695 0.073292433 0.064223211 0.075067396 0.075057901 0.062342975 0.069971145
## [243] 0.068828342 0.059362935 0.069789923 0.045923878 0.069787233 0.064884659 0.075425613 0.073811129 0.064326121 0.075068488 0.075423107
## [254] 0.075069354 0.070065170 0.075431913
## 
## $IJ
##                 I            J       Total
## T     0.010711561 -0.006500721 0.004210840
## SOIL  0.001645553  0.002456721 0.004102275
## GDD   0.007884330 -0.003982500 0.003901830
## AWSSM 0.001606291  0.002667114 0.004273405
## AET   0.013946197  0.003195376 0.017141572
## D     0.011081196  0.006057683 0.017138879
## RD    0.010961011  0.014896611 0.025857621
## VPD   0.017595776 -0.001075942 0.016519833
## 
## $I.perc
##       ind.exp.var
## T       14.200304
## SOIL     2.181508
## GDD     10.452247
## AWSSM    2.129458
## AET     18.488457
## D       14.690328
## RD      14.530999
## VPD     23.326699
## 
## $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.0002649891 0.0283302686 0.0001220891 0.0282875421 0.0411654229 0.0121652223 0.0326482691 0.0053336440 0.0538501377
##  [11] 0.0199139291 0.0538836429 0.0419257381 0.0470049142 0.0346279483 0.0261097834 0.0520726544 0.0284291480 0.0533511303 0.0311207253
##  [21] 0.0473664733 0.0300472128 0.0521072976 0.0422434652 0.0436711529 0.0341929409 0.0233034004 0.0534388751 0.0311106057 0.0472498960
##  [31] 0.0300314755 0.0500711585 0.0570800686 0.0489359500 0.0334009658 0.0163131317 0.0340963320 0.0668533513 0.0538931852 0.0606912054
##  [41] 0.0596266443 0.0711501325 0.0666498602 0.0668042812 0.0611358809 0.0651323405 0.0533929190 0.0461542849 0.0608284180 0.0595974557
##  [51] 0.0710127808 0.0665661250 0.0610709286 0.0570831048 0.0568622879 0.0512440759 0.0471602062 0.0347211182 0.0521184510 0.0595170024
##  [61] 0.0569310096 0.0694125562 0.0633414468 0.0545322707 0.0311208625 0.0489506838 0.0300523183 0.0533519855 0.0643320168 0.0536472807
##  [71] 0.0636132102 0.0312792600 0.0656039613 0.0596537368 0.0569082601 0.0692781710 0.0632701808 0.0590038280 0.0570931151 0.0551408218
##  [81] 0.0485997598 0.0437316707 0.0342083875 0.0534389032 0.0643392356 0.0537188674 0.0634038728 0.0312638974 0.0654536698 0.0586017826
##  [91] 0.0500846203 0.0572116975 0.0506513911 0.0668883984 0.0752706861 0.0741080965 0.0845274173 0.0805903989 0.0615365748 0.0596525478
## [101] 0.0716479887 0.0667985632 0.0656461704 0.0738381905 0.0719294459 0.0712824528 0.0671154980 0.0712003086 0.0753433222 0.0740253805
## [111] 0.0843293131 0.0804354611 0.0793438851 0.0757685820 0.0764308676 0.0693431100 0.0654404740 0.0536335732 0.0656944977 0.0738090150
## [121] 0.0719617526 0.0711374467 0.0670245430 0.0710641336 0.0641455128 0.0623247101 0.0583140469 0.0537295033 0.0604135878 0.0569437352
## [131] 0.0699049417 0.0634454732 0.0635985828 0.0724918373 0.0692369984 0.0697882441 0.0638841025 0.0694387917 0.0545383594 0.0643493475
## [141] 0.0547259680 0.0643494525 0.0312804049 0.0663699183 0.0694414088 0.0545188484 0.0695896414 0.0660606367 0.0636543584 0.0724691305
## [151] 0.0692811515 0.0696404245 0.0638067072 0.0693033206 0.0626035740 0.0600211643 0.0575592617 0.0523412873 0.0694262175 0.0546294216
## [161] 0.0695644060 0.0658910255 0.0617640417 0.0755936386 0.0745091872 0.0858083994 0.0812734106 0.0815812091 0.0882281994 0.0872525615
## [171] 0.0845281985 0.0807322702 0.0846774023 0.0659306910 0.0738679368 0.0720016768 0.0717711944 0.0672451913 0.0716900717 0.0738662546
## [181] 0.0723834759 0.0741757789 0.0715832178 0.0815870716 0.0881593948 0.0872345982 0.0843294415 0.0805675220 0.0844809940 0.0823884755
## [191] 0.0810044849 0.0774491270 0.0711377005 0.0738352449 0.0724248632 0.0741556586 0.0714336404 0.0642388251 0.0639696565 0.0725097107
## [201] 0.0693683921 0.0702703572 0.0639714768 0.0699402652 0.0726556914 0.0697561232 0.0725332182 0.0698559767 0.0694450612 0.0556502524
## [211] 0.0696250675 0.0667766245 0.0700539741 0.0726286530 0.0698123179 0.0725138381 0.0697063436 0.0629954229 0.0700271600 0.0815921628
## [221] 0.0885142410 0.0872678632 0.0858084230 0.0813873428 0.0859398644 0.0882769422 0.0873758987 0.0888975036 0.0847804784 0.0738973726
## [231] 0.0724703388 0.0741884185 0.0720376841 0.0745711141 0.0882128243 0.0873576430 0.0888406269 0.0845724431 0.0823901294 0.0745515494
## [241] 0.0726768792 0.0699074311 0.0725456401 0.0703197263 0.0729614688 0.0700781673 0.0729426592 0.0885598251 0.0873874226 0.0891054039
## [251] 0.0860088823 0.0889987827 0.0745779571 0.0889382282 0.0729679358 0.0891916056
## 
## $IJ
##                 I             J        Total
## T     0.014064952 -0.0137999626 0.0002649891
## SOIL  0.011416712  0.0169135569 0.0283302686
## GDD   0.012495297 -0.0123732083 0.0001220891
## AWSSM 0.011371858  0.0169156842 0.0282875421
## AET   0.015756334  0.0254090884 0.0411654229
## D     0.006087209  0.0060780134 0.0121652223
## RD    0.013559411  0.0190888578 0.0326482691
## VPD   0.004439832  0.0008938118 0.0053336440
## 
## $I.perc
##       ind.exp.var
## T       15.769367
## SOIL    12.800209
## GDD     14.009499
## AWSSM   12.749920
## AET     17.665715
## D        6.824867
## RD      15.202564
## VPD      4.977859
## 
## $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

CWD Forest

ecdf_fn <- ecdf(ignitions_forest$D)
summary(ecdf_fn)
## Empirical CDF:     151 unique values with summary
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.7483  0.8913  0.8199  0.9778  1.0000
estimated <- 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$dec))

##plot(lm)

summary(lm)
## 
## Call:
## lm(formula = log(proportions) ~ quants, data = estimated)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.42501 -0.06183 -0.02935  0.05987  0.25933 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -4.64736    0.04396 -105.72   <2e-16 ***
## quants       4.44971    0.05166   86.14   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1022 on 97 degrees of freedom
## Multiple R-squared:  0.9871, Adjusted R-squared:  0.987 
## F-statistic:  7420 on 1 and 97 DF,  p-value: < 2.2e-16
a <- 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) = 0.0095869 e^{4.4497083 x}\end{align*}\)

VPD Forest

ecdf_fn <- ecdf(ignitions_forest$VPD)
summary(ecdf_fn)
## Empirical CDF:     152 unique values with summary
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.01587 0.74339 0.90225 0.81542 0.97276 1.00000
estimated <- 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$dec))

##plot(lm)

summary(lm)
## 
## Call:
## lm(formula = log(proportions) ~ quants, data = estimated)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.47105 -0.07338 -0.02498  0.06839  0.26338 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -4.55048    0.04428 -102.76   <2e-16 ***
## quants       4.36308    0.05233   83.38   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1056 on 97 degrees of freedom
## Multiple R-squared:  0.9862, Adjusted R-squared:  0.9861 
## F-statistic:  6952 on 1 and 97 DF,  p-value: < 2.2e-16
a <- 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) = 0.0105621 e^{4.3630793 x}\end{align*}\)

CWD Non-Forest

ecdf_fn <- ecdf(ignitions_nonforest$D)
summary(ecdf_fn)
## Empirical CDF:     86 unique values with summary
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.7746  0.9014  0.8193  0.9772  1.0000
estimated <- 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$dec))

##plot(lm)

summary(lm)
## 
## Call:
## lm(formula = log(proportions) ~ quants, data = estimated)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.40749 -0.08412 -0.03126  0.10090  0.28049 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -4.42528    0.05212  -84.90   <2e-16 ***
## quants       4.18927    0.06113   68.53   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1281 on 97 degrees of freedom
## Multiple R-squared:  0.9798, Adjusted R-squared:  0.9796 
## F-statistic:  4697 on 1 and 97 DF,  p-value: < 2.2e-16
a <- 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) = 0.0119709 e^{4.1892655 x}\end{align*}\)

VPD Non-Forest

ecdf_fn <- ecdf(ignitions_nonforest$VPD)
summary(ecdf_fn)
## Empirical CDF:     84 unique values with summary
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.04694 0.70921 0.88820 0.80687 0.98529 1.00000
estimated <- 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$dec))

##plot(lm)

summary(lm)
## 
## Call:
## lm(formula = log(proportions) ~ quants, data = estimated)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.75860 -0.09779 -0.01807  0.10064  0.34912 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -4.27903    0.06204  -68.97   <2e-16 ***
## quants       4.05557    0.07336   55.28   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.158 on 97 degrees of freedom
## Multiple R-squared:  0.9692, Adjusted R-squared:  0.9689 
## F-statistic:  3056 on 1 and 97 DF,  p-value: < 2.2e-16
a <- 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.0138562 e^{4.055571 x}\end{align*}\)