Rockies Fire Analysis
Stephen Huysman
2023-05-15
Settings
library(tidyverse)
library(multidplyr)
library(zoo)
library(slider)
library(pROC)
plot_theme <- theme_bw()
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 (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.
- 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.
- 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
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.008427108 0.004936385 0.009037155 0.004935788 0.012310342 0.006035703 0.009214246 0.013261302 0.009225352 0.009295240
## [12] 0.009224873 0.013244905 0.008441034 0.011966691 0.020465597 0.009742468 0.004972326 0.019055532 0.006388953 0.010556672 0.013385038
## [23] 0.009741945 0.013540408 0.009113515 0.012363735 0.017499394 0.019053578 0.006388662 0.010555961 0.013384761 0.014005101 0.018671514
## [34] 0.016538550 0.010593574 0.015746903 0.014350564 0.010052489 0.009334237 0.020625296 0.010410051 0.012273349 0.020490646 0.010051852
## [45] 0.014022543 0.009333422 0.012686966 0.020558477 0.020622368 0.010408815 0.012272934 0.020490452 0.014140549 0.018759996 0.030042884
## [56] 0.012094393 0.022181921 0.020959069 0.009890331 0.020100485 0.011323698 0.012622831 0.017544607 0.019337679 0.006442870 0.010719812
## [67] 0.013611564 0.021585140 0.021663999 0.019162178 0.010824285 0.019669647 0.014419079 0.020097922 0.011322038 0.012622419 0.017544426
## [78] 0.014008122 0.018689573 0.024932752 0.012624578 0.019405112 0.017509401 0.021581013 0.021662266 0.019160417 0.010823872 0.019666416
## [89] 0.014418852 0.018710233 0.017632952 0.019034458 0.016357990 0.010259073 0.021191041 0.011438559 0.012977680 0.020584909 0.020911054
## [100] 0.010552478 0.012479507 0.021040210 0.021600649 0.025081151 0.035512468 0.013417362 0.024461963 0.020997961 0.021187625 0.011436809
## [111] 0.012977179 0.020584701 0.014762385 0.019346007 0.030051487 0.012867038 0.022296095 0.021010353 0.021596564 0.025077602 0.035507554
## [122] 0.013416011 0.024458837 0.020997716 0.019463247 0.030050911 0.030363309 0.022986246 0.020350969 0.011555750 0.012864369 0.017808488
## [133] 0.021623118 0.024248963 0.029473875 0.014160294 0.022129803 0.017557373 0.021920853 0.022086691 0.019461820 0.010991111 0.020184353
## [144] 0.014692997 0.025642391 0.026666162 0.022063583 0.019717133 0.021618858 0.024245953 0.029471080 0.014158596 0.022127067 0.017557188
## [155] 0.019011742 0.025057724 0.024952916 0.019496130 0.025637325 0.026659776 0.022061357 0.019713995 0.019421640 0.021597862 0.011710518
## [166] 0.013300368 0.021195109 0.022320590 0.025511639 0.035512930 0.014305646 0.024626476 0.021049958 0.021932508 0.025570903 0.036610679
## [177] 0.013667411 0.025265581 0.021564236 0.025998950 0.039557019 0.035936740 0.026042107 0.022315664 0.025507714 0.035508007 0.014303877
## [188] 0.024623144 0.021049702 0.019931628 0.030060523 0.030364219 0.023043705 0.025994047 0.039547443 0.035931733 0.026038209 0.030399658
## [199] 0.021978263 0.024655184 0.029898542 0.014504849 0.022643128 0.017816768 0.025716139 0.034304175 0.029487458 0.022612825 0.026198244
## [210] 0.027339215 0.022477562 0.020237089 0.027466568 0.025711237 0.034297386 0.029484641 0.022609875 0.025063332 0.027460295 0.022805978
## [221] 0.026135079 0.036629085 0.014706104 0.025530368 0.021663161 0.026561475 0.039571468 0.035939429 0.026118476 0.026545504 0.041013052
## [232] 0.037054918 0.026930559 0.040994087 0.026555877 0.039561782 0.035934445 0.026114434 0.030401072 0.040983429 0.026250143 0.035047239
## [243] 0.029920009 0.023107165 0.034526354 0.028179799 0.034519460 0.027276062 0.041079360 0.037058752 0.027080195 0.040994099 0.042555464
## [254] 0.040983445 0.035251869 0.042573433
##
## $IJ
## I J Total
## T 0.005322251 0.003104858 0.008427108
## SOIL 0.002470836 0.002465550 0.004936385
## GDD 0.002734463 0.006302691 0.009037155
## AWSSM 0.002466903 0.002468885 0.004935788
## AET 0.011714893 0.000595449 0.012310342
## D 0.003488227 0.002547476 0.006035703
## RD 0.003231900 0.005982346 0.009214246
## VPD 0.011143960 0.002117341 0.013261302
##
## $I.perc
## ind.exp.var
## T 12.501342
## SOIL 5.803703
## GDD 6.422934
## AWSSM 5.794465
## AET 27.516910
## D 8.193437
## RD 7.591354
## VPD 26.175855
##
## $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.0091759377 0.0070578192 0.0091573380 0.0070558505 0.0002212463 0.0121725412 0.0003394020 0.0081664319 0.0130208835
## [11] 0.0092058079 0.0130208477 0.0101921723 0.0127991742 0.0091887468 0.0092502453 0.0130072508 0.0071686312 0.0123413350 0.0121768178
## [21] 0.0073307854 0.0108000447 0.0130072242 0.0101744966 0.0127950668 0.0091707703 0.0092349736 0.0123390791 0.0121768531 0.0073281008
## [31] 0.0107996317 0.0131174381 0.0008086287 0.0082007617 0.0134471780 0.0121756939 0.0091245395 0.0130443830 0.0130209319 0.0136020491
## [41] 0.0132520251 0.0138056589 0.0134913831 0.0130443040 0.0102096515 0.0128110114 0.0092159427 0.0092830238 0.0136015735 0.0132522581
## [51] 0.0138049335 0.0134913648 0.0131189940 0.0106479253 0.0102493225 0.0138449593 0.0137294233 0.0094390500 0.0130072746 0.0135912294
## [61] 0.0132437248 0.0137973435 0.0134712229 0.0124001804 0.0122040693 0.0075106829 0.0108109399 0.0140647833 0.0124035697 0.0128404556
## [71] 0.0134790620 0.0121873232 0.0128156631 0.0135907596 0.0132439632 0.0137966274 0.0134712234 0.0131186651 0.0106368395 0.0102295307
## [81] 0.0138445118 0.0137225371 0.0094307323 0.0140650282 0.0124008670 0.0128392374 0.0134790823 0.0121873855 0.0128142903 0.0139150120
## [91] 0.0135243696 0.0097400871 0.0135714925 0.0130446624 0.0136390383 0.0132692427 0.0138120291 0.0135075088 0.0136072062 0.0132594081
## [101] 0.0138249491 0.0134913879 0.0140990774 0.0140527759 0.0138757154 0.0143162146 0.0141301809 0.0138250590 0.0136384760 0.0132694507
## [111] 0.0138112834 0.0135074483 0.0131338576 0.0106516373 0.0102634815 0.0138450392 0.0137321233 0.0094612297 0.0140993013 0.0140517501
## [121] 0.0138753159 0.0143161457 0.0141306094 0.0138244208 0.0139372859 0.0139354981 0.0107017718 0.0140170044 0.0135962472 0.0132515296
## [131] 0.0138161063 0.0134712249 0.0140968034 0.0140439959 0.0138578363 0.0143140146 0.0141188073 0.0138149669 0.0140681762 0.0124871664
## [141] 0.0128660846 0.0134802875 0.0122111363 0.0128986104 0.0146697173 0.0142923190 0.0135280572 0.0137667898 0.0140970301 0.0140429782
## [151] 0.0138574520 0.0143139506 0.0141192441 0.0138143354 0.0139369933 0.0139303983 0.0106935498 0.0140154348 0.0146696501 0.0142927151
## [161] 0.0135262191 0.0137666958 0.0139503143 0.0136462812 0.0132750903 0.0138325430 0.0135076150 0.0141274235 0.0140675872 0.0139030902
## [171] 0.0143171739 0.0141357543 0.0138319166 0.0141017776 0.0140777972 0.0138790165 0.0143167991 0.0141583097 0.0138396417 0.0147418343
## [181] 0.0147532151 0.0140693230 0.0144540800 0.0141275929 0.0140665159 0.0139026120 0.0143170994 0.0141361694 0.0138312579 0.0139378887
## [191] 0.0139398286 0.0107051287 0.0140171073 0.0147417163 0.0147537258 0.0140683999 0.0144542109 0.0141266888 0.0140995771 0.0140685927
## [201] 0.0138608841 0.0143145191 0.0141484149 0.0138292510 0.0147402305 0.0147445698 0.0140586131 0.0144500047 0.0146699914 0.0143028447
## [211] 0.0135985250 0.0137698305 0.0146759240 0.0147401156 0.0147450890 0.0140576936 0.0144501408 0.0141249905 0.0146759082 0.0141289765
## [221] 0.0140948699 0.0139078711 0.0143178609 0.0141620976 0.0138475434 0.0147478773 0.0147662522 0.0140847624 0.0144550587 0.0147426511
## [231] 0.0147695798 0.0140893800 0.0144566369 0.0149005732 0.0147477385 0.0147667283 0.0140837944 0.0144551845 0.0141274320 0.0149007674
## [241] 0.0147410040 0.0147615676 0.0140785337 0.0144527856 0.0148960003 0.0146759350 0.0148961978 0.0147490109 0.0147804028 0.0141067966
## [251] 0.0144574240 0.0149068530 0.0149032280 0.0149070283 0.0148987588 0.0149090117
##
## $IJ
## I J Total
## T 0.0021381942 0.0070377435 0.0091759377
## SOIL 0.0020187559 0.0050390633 0.0070578192
## GDD 0.0021284689 0.0070288691 0.0091573380
## AWSSM 0.0020182382 0.0050376123 0.0070558505
## AET 0.0008341641 -0.0006129178 0.0002212463
## D 0.0035013999 0.0086711413 0.0121725412
## RD 0.0004521797 -0.0001127777 0.0003394020
## VPD 0.0018176108 0.0063488211 0.0081664319
##
## $I.perc
## ind.exp.var
## T 14.341622
## SOIL 13.540508
## GDD 14.276392
## AWSSM 13.537035
## AET 5.595033
## D 23.485124
## RD 3.032929
## VPD 12.191357
##
## $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: 110 unique values with summary
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.2227 0.9040 0.9532 0.9097 0.9820 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$dec))
##plot(lm)
summary(lm)##
## Call:
## lm(formula = log(proportions) ~ quants, data = estimated)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.53762 -0.21281 0.02197 0.19476 0.92328
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -7.9712 0.2174 -36.67 <2e-16 ***
## quants 7.6587 0.2358 32.48 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2563 on 97 degrees of freedom
## Multiple R-squared: 0.9158, Adjusted R-squared: 0.9149
## F-statistic: 1055 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.4527061\times 10^{-4} e^{7.6587032 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.1170 0.8892 0.9428 0.9074 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$dec))
##plot(lm)
summary(lm)##
## Call:
## lm(formula = log(proportions) ~ quants, data = estimated)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.48946 -0.14994 0.05073 0.14202 0.69295
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -8.4876 0.1968 -43.14 <2e-16 ***
## quants 8.2272 0.2137 38.49 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2189 on 97 degrees of freedom
## Multiple R-squared: 0.9386, Adjusted R-squared: 0.9379
## F-statistic: 1482 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.060083\times 10^{-4} e^{8.227219 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.9107 0.9628 0.9099 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$dec))
##plot(lm)
summary(lm)##
## Call:
## lm(formula = log(proportions) ~ quants, data = estimated)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.45075 -0.24819 0.03401 0.21845 0.78265
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -7.3061 0.2019 -36.19 <2e-16 ***
## quants 6.9122 0.2182 31.67 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2703 on 97 degrees of freedom
## Multiple R-squared: 0.9118, Adjusted R-squared: 0.9109
## F-statistic: 1003 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.7146057\times 10^{-4} e^{6.9121875 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.1941 0.8893 0.9601 0.9015 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$dec))
##plot(lm)
summary(lm)##
## Call:
## lm(formula = log(proportions) ~ quants, data = estimated)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.37246 -0.17610 -0.00724 0.19924 0.54526
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6.8126 0.1522 -44.75 <2e-16 ***
## quants 6.4467 0.1661 38.80 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2239 on 97 degrees of freedom
## Multiple R-squared: 0.9395, Adjusted R-squared: 0.9389
## F-statistic: 1506 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.0010998 e^{6.4466519 x}\end{align*}\)