03_NKNoManagement.Rmd

---
title: "Nunery Keeton Revisited"
author: "Nikolaus Bates-Haus"
output:
  pdf_document: default
  html_document: default
editor_options:
  markdown:
    wrap: 72
---

```{r setup, include=FALSE}
library(tidyverse)
library(reshape2) # for melt()
library(RSQLite)
library(measurements)
```

```{r source-functions}
source('R/functions.R')
```

# Reproducing the NoManagement Scenario

## Configuring FVS

FVS Software with `Release date: 20240401` was used.

FIA data version:

```{r}
fia <- DBI::dbConnect(RSQLite::SQLite(), 'data/raw/SQLite_FIADB_ENTIRE.db')
fia_ref_fiadb_version <- tbl(fia, 'REF_FIADB_VERSION')
version <- fia_ref_fiadb_version |>
  filter(INSTALL_TYPE == 'RELEASE') |>
  arrange(desc(CREATED_DATE)) |>
  head(1) |>
  collect()
DBI::dbDisconnect(fia)
remove(fia, fia_ref_fiadb_version)
knitr::kable(version)
```

### Stands

```{r nk_to_fia}
nk_to_fia <- read_csv(
  "data/intermediate/nk_to_fia.csv", 
  col_types = cols(
    `FIA plot code` = col_character(), 
    NK_INVYR = col_character(),
    STATECD = col_integer(), 
    INVYR = col_integer(),
    CYCLE = col_character(), 
    SUBCYCLE = col_character(),
    CONDID = col_character(), 
    UNITCD = col_character(),
    COUNTYCD = col_integer(), 
    PLOT = col_integer()
  )
)
```

#### Via the UI

Stands can be found in the FVS UI using FVS_STAND_ID from the table
above.

-   Inventory Data Tables = FVS_StandInit_Plot
-   Variants = ne: Northeast
-   Groups: All_FIA_Plots
-   Find stand(s): = 239901700110

Clicking `Find` will then filter to the single matching stand, which can
be added to the run.

#### Automatic

What FVS calls "StandCN" is what FIA calls "PLOT.CN". Augment nk_to_fia
to include this information.

```{r stand_cn}
fia <- DBI::dbConnect(RSQLite::SQLite(), 'data/raw/SQLite_FIADB_ENTIRE.db')

stand_cns <- tbl(fia, 'PLOT') |>
  left_join(nk_to_fia, by = join_by(STATECD, COUNTYCD, PLOT, INVYR), copy = TRUE) |>
  select(CN, STATECD, COUNTYCD, PLOT, INVYR, MEASYEAR) |>
  rename(STAND_CN = CN) |>
  collect()

nk_to_fia <- nk_to_fia |>
  left_join(stand_cns, by = join_by(STATECD, COUNTYCD, PLOT, INVYR))

DBI::dbDisconnect(fia)
rm(fia, stand_cns)
```

### Time

FVS will not run prior to the inventory year; it's unclear what NK did
to get their chart to go back prior to 2005. We will instead use 2005 as
the common starting year.

NK ends in 2164; we adjust this to 2165 so that all cycles are 10 years
long. Note that FVS does not include the ending year unless it is at the
start of a new cycle, so we set the common ending year to 2166.

All events in NK take place on 10 year intervals, therefore we select 10
years as the reporting interval.

### Regeneration

NK regeneration rates are in NK Table 4, in seedlings per hectare.

In addition to scenario-specific regeneration, NK 2.4 states,
"Background regeneration rates (intermediate to shade tolerant species
only), input at 10 year intervals, emulated natural regeneration within
stands, independent of forest management activities."

This background regeneration rate is in NK Table 4.

Note that if *any* establishment is specified, whether Natural or Plant,
FVS defaults to doing site prep. It generates a mix of mechanical and burn site
prep based on site characteristics. For grow-only plots, it's important to turn
this off by explicitly

```{r}
table4 <- read_csv(
  "data/raw/NK_Table_4.csv",
  col_types = cols(
    `Management scenario` = col_character(),
    .default = col_number(),
  )
)

# Table4 uses scientific names; translate to get FVS Species Code
fvs_species_codes <- read_csv(
    "data/raw/FVSne_Overview_Table_3.2.1.csv",
    col_types = "iiciccc"
  ) |>
  select("Scientific Name", "Species Code", "Common Name")

# table4 has one column per species. We wish to add observations per species
# with different types, so pivot the table to have one row per species.
table4_rot <- table4 |>
  pivot_longer(cols = !`Management scenario`) |>
  pivot_wider(names_from = `Management scenario`) |>
  rename(`Scientific Name` = name) |>
  left_join(fvs_species_codes, by = join_by(`Scientific Name`)) |>
  # Values in NK are seedlings per hectare; FVS needs seedlings per acre.
  # To convert rates, do the reverse conversion, i.e. to convert from
  # per-hectare to per-acre, we convert using the inverse of the ratio
  # of hectare-to-acre, which is the same as the ratio from acre-to-hectare.
  mutate(Clearcut = round(conv_unit(Clearcut, 'acre', 'hectare'))) |>
  mutate(Shelterwood = round(conv_unit(Shelterwood, 'acre', 'hectare'))) |>
  mutate(`ITS_Low Retention` = round(conv_unit(`ITS_Low Retention`, 'acre', 'hectare'))) |>
  mutate(`ITS_High Retention` = round(conv_unit(`ITS_High Retention`, 'acre', 'hectare'))) |>
  mutate(Background = round(conv_unit(Background, 'acre', 'hectare')))


write_csv(table4_rot, 'data/intermediate/nk_regen.csv')

knitr::kable(
  table4_rot |>
    filter(!is.na(Background)) |>
    select(`Scientific Name`, Background, `Species Code`, `Common Name`) |>
    arrange(desc(Background))
)
```

#### Via the UI

In FVS, go to Simulate -\> Components -\> Management From Categories,
select "Planting & Natural Regeneration" From Components, select
"Plant/Natural with Partial Estab Model"

- Component title: "Baseline Regen"
- Schedule the date of disturbance: Schedule by condition
- Create a condition: Every cycle, every even cycle, ...
- Condition title: Every cycle
- Years before condition can become true again: 0
- The modulus of cycle number: 1 = every cycle
- Number of years after condition is found true: 0
- Sprouting: On

Note that this leaves Percent survival at 100, and average age and
height empty.

In freeform, this looks like:

```         
Estab              0
Sprout
Natural            0        SM       200      100.                             0
Natural            0        AB       100      100.                             0
Natural            0        EH        25      100.                             0
Natural            0        RS        25      100.                             0
Natural            0        YB        25      100.                             0
Natural            0        RM        25      100.                             0
```

Switch to freeform and add all the rows, fix the titles.

Select `Save in run` to apply to `Grp: All_FIA_Plots`.

#### Automatic

```{r background_regen}
background_regen <- table4_rot |>
  rename(species = `Species Code`, density = Background) |>
  filter(!is.na(density)) |>
  fvs_Estab()
```

### Carbon

To set the carbon calculations to be metric (Metric tons of carbon per
hectare), go to Simulate -\> Components -\> Keywords, and select:

-   Extensions: Fire and Fuels Extension
-   Keywords: CarbCalc: Set the carbon accounting parameters.
-   Component title: CarbCalc: Metric
-   Biomass predictions: 1 = Use Jenkins biomass predictions
-   Units: 1 = Metric (metric tons carbon/hectare)
-   Note: Annual root decay rate (proportion per year) remains at its
    default value of 0.0425.

Select `Save in run` to apply to `Grp: All_FIA_Plots`.

## Keywordfile Generation

```{r keywordfile_section}
keywordfile_section <- function(Title, MgmtId, StandID, StandCN, FirstYear, LastYear) {
  TimeConfig <- fvs_TimeConfig(FirstYear, LastYear, 10)
  c(
    fvs_kwd0("StdIdent"),
    paste0(StandID, " ", Title),
    fvs_kwd0("StandCN"),
    StandCN,
    fvs_kwd0("MgmtId"),
    MgmtId,
    TimeConfig,
    fvs_kwd0("FMIn"), # Fire and Fuels Extension
    fvs_kwd0("CarbRept"),
    fvs_kwd0("CarbCut"),
    fvs_kwd5("CarbCalc", 1, 1, 0.0425, 9, 11),
    fvs_kwd0("FuelOut"),
    fvs_kwd0("FuelRept"),
    fvs_kwd0("End"), # FMIn
    fvs_kwd0("Database"), # Database extension
    fvs_kwd0("DSNIn"),
    "SQLite_FIADB_ENTIRE.db",
    fvs_kwd0("StandSQL"),
    "SELECT * FROM FVS_StandInit_Plot WHERE Stand_CN = '%Stand_CN%'",
    fvs_kwd0("EndSQL"), # StandSQL
    fvs_kwd0("TreeSQL"),
    "SELECT * FROM FVS_TreeInit_Plot WHERE Stand_CN = '%Stand_CN%'",
    fvs_kwd0("EndSQL"), # TreeSQL
    fvs_kwd0("DSNOut"),
    paste0("FVS_", Title, "_", MgmtId, ".db"),
    fvs_kwd1("Summary",  2),
    fvs_kwd2("Computdb", 0, 1),
    fvs_kwd1("MisRpts",  2),
    fvs_kwd1("CarbReDB", 2),
    fvs_kwd1("FuelReDB", 2),
    fvs_kwd1("FuelsOut", 2),
    fvs_kwd0("End"), # Database
    background_regen,
    fvs_kwd0("Process"),
    recursive = TRUE
  )
}

```

```{r FVS_NKByPlot_NONE.key, include=FALSE}
Title <- "NKByPlot"
MgmtId <- "NONE"
filename <- paste0("data/fvs/FVS_", Title, "_", MgmtId, ".key")
unlink(filename)
apply(
  nk_to_fia,
  1,
  function(row) {
    write_lines(
      keywordfile_section(
        Title,
        MgmtId,
        row['FIA plot code'],
        row['STAND_CN'],
        row['MEASYEAR'],
        2165
      ),
      filename,
      append = TRUE
    )
  }
)
write_lines("Stop", filename, append = TRUE)
```


## Run FVS

The main output will be Carbon and fuels. Stand visualization, Tree
lists, and Inventory statistics may also be interesting.

```
REM Delete destination databases. FVS will otherwise append to them, leaving previous rows.
for %%p in (FVS_%1_%2*.key) do del %%~np.db

REM Start FVSne in the background on each partition
for %%p in (FVS_%1_%2*.key) do start \FVS\FVSbin\FVSne.exe --keywordfile=%%p
```

## Carbon Projection

### By Plot

Read FVS output for carbon storage by plot. FVS column names are a bit
idiosyncratic, so clean those up.

NK states that Fig 2 sums carbon from aboveground live, standing dead,
and down dead. We plot this, as well as just aboveground live and stand total
carbon for comparison.

```{r FVS_Carbon}
fvs_out_db <- DBI::dbConnect(
  RSQLite::SQLite(),
  paste0('data/fvs/FVS_', Title, '_', MgmtId, '.db')
)

FVS_Carbon <- tbl(fvs_out_db, 'FVS_Carbon') |> collect()
FVS_Summary2_East <- tbl(fvs_out_db, 'FVS_Summary2_East') |> collect()

dbDisconnect(fvs_out_db)
remove(fvs_out_db)

NoManagement_Carbon_ByPlot <- FVS_Carbon |>
  rename(Aboveground_Live = Aboveground_Total_Live) |>
  rename(Down_Dead = Forest_Down_Dead_Wood) |>
  mutate(Aboveground_Dead = Standing_Dead + Down_Dead) |>
  mutate(Aboveground_Carbon = Aboveground_Live + Aboveground_Dead) |>
  select(StandID, Year, Total_Stand_Carbon, Aboveground_Carbon,
         Aboveground_Live, Aboveground_Dead, Standing_Dead, Down_Dead)
knitr::kable(NoManagement_Carbon_ByPlot)
```

Plot these against NK Fig 2 "NoManagement" scenario. NK Fig 2 has three
different scales on the horizontal axis:

1.  From 1995 to 2005, major ticks are every 5 years
2.  From 2005 to 2155, major ticks are every 10 years
3.  From 2155 to 2164, major ticks are every 9 years

NK states that data prior to 2005 is projected from mean growth rate. To
simplify the chart, we omit data prior to 2005, place major tick marks
every 10 years, and nudge the end date from 2164 to 2165.

Note that NK Fig2 year 2164 shows an anomalous reduction in carbon
storage in the NoManagement scenario; it might be the right thing to do
to omit 2164 / 2165 as well.

```{r}
fig2 <- read_csv(
  "data/intermediate/nk_fig2_reconstructed.csv",
  col_types = cols(
    Year = col_integer(),
    NoManagement = col_double(),
    ClearcutHigh = col_double(),
    ClearcutLow = col_double(),
    ShelterwoodHigh = col_double(),
    ShelterwoodLow = col_double(),
    ITS_LowHigh = col_double(),
    ITS_LowLow = col_double(),
    ITS_HighHigh = col_double(),
    ITS_HighLow = col_double()
  )
)

NoManagement_Carbon_ByPlot_ByYear <- NoManagement_Carbon_ByPlot |>
  group_by(Year) |>
  summarize(
    Total_Carbon=mean(Total_Stand_Carbon),
    Aboveground_Carbon=mean(Aboveground_Carbon),
    Aboveground_Live=mean(Aboveground_Live),
    .groups = "keep"
  ) |>
  filter(Year >= 2005) |>
  # Join in the NoManagement scenario from fig2 for comparison
  full_join(fig2 |> select(Year,NoManagement), by=join_by(Year)) |>
  rename(NK_Fig2_NoManagement = NoManagement)

ggplot(
    data = melt(NoManagement_Carbon_ByPlot_ByYear, id.vars = "Year"),
    mapping = aes(x = Year, y = value, color = variable, linetype = variable)
  ) +
  ggtitle("FVS Carbon Projection for NK Plots with no management") +
  ylab(bquote("Carbon " ~(Mg %*% ha^{-1}))) +
#  theme(legend.title = element_blank()) +
  theme(legend.title = element_blank(), legend.position = "bottom") +
  scale_linetype_manual(values = c("solid", "solid", "solid", "dashed")) +
  geom_line() +
  coord_cartesian(xlim = c(2005, 2165), ylim = c(0, 320)) +
  scale_x_continuous(breaks=seq(2005,2165,20)) +
  scale_y_continuous(breaks=seq(0,320,20))
```

From the chart it seems evident that aboveground live carbon is the best match
for the NoManagement scenario in NK Fig 2. We confirm this by computing
RMS error for the residuals between our FVS runs and the NK Fig 2 data:

```{r}
Plots_RMSE <- data.frame(
  Total_Carbon = sqrt(mean(
    (NoManagement_Carbon_ByPlot_ByYear$Total_Carbon - NoManagement_Carbon_ByPlot_ByYear$NK_Fig2_NoManagement)^2
  )),
  Aboveground_Carbon = sqrt(mean(
    (NoManagement_Carbon_ByPlot_ByYear$Aboveground_Carbon - NoManagement_Carbon_ByPlot_ByYear$NK_Fig2_NoManagement)^2
  )),
  Aboveground_Live = sqrt(mean(
    (NoManagement_Carbon_ByPlot_ByYear$Aboveground_Live -
       NoManagement_Carbon_ByPlot_ByYear$NK_Fig2_NoManagement)^2
  ))
)
ggplot(
    data = melt(Plots_RMSE, id = NULL, variable.name = "Series", value.name = "RMSE"),
    mapping = aes(x = Series, y = RMSE)
  ) +
  ggtitle("FVS Carbon Projection by FIA Plot") +
  geom_col()
```

Note that NK describe using Aboveground Carbon, including both live and dead,
for their chart, but the best fit is from Aboveground Live, with RMSE \~= 12.

```{r}
Plots_RMSE |> select(Aboveground_Live)
```

### By Subplot

To confirm that NK used the FVS_StandInit_Plot table, we examine the
RMSE between projections from other tables and the NK Fig 2 values.

We repeat the above FVS run, but select table FVS_PlotInit_Plot and
select all subplots using the translated stand IDs. We apply the same
configuration to all stands as we did for FVS_StandInit_Plot. We load
and clean the results in the same manner.

```{r}
NoManagement_Carbon_BySubplot_ByYear <- read_csv(
    "data/fvs/FVS_NKBySubplot_NONE_Carbon.csv",
    col_types = cols(StandID = col_character())
  ) |>
  rename(Total_Carbon = Total_Stand_Carbon) |>
  rename(Aboveground_Live = Aboveground_Total_Live) |>
  rename(Down_Dead = Forest_Down_Dead_Wood) |>
  mutate(Aboveground_Dead = Standing_Dead + Down_Dead) |>
  mutate(Aboveground_Carbon = Aboveground_Live + Aboveground_Dead) |>
  select(StandID, Year, Total_Carbon, Aboveground_Carbon,
         Aboveground_Live, Aboveground_Dead, Standing_Dead, Down_Dead) |>
  group_by(Year) |>
  summarize(
    Total_Carbon=mean(Total_Carbon),
    Aboveground_Carbon=mean(Aboveground_Carbon),
    Aboveground_Live=mean(Aboveground_Live),
    .groups = "keep"
  ) |>
  filter(Year >= 2005) |>
  full_join(fig2 |> select(Year,NoManagement), by=join_by(Year)) |>
  rename(NK_Fig2_NoManagement = NoManagement)

ggplot(
    data = melt(NoManagement_Carbon_BySubplot_ByYear, id.vars = "Year"),
    mapping = aes(x = Year, y = value, color = variable)
  ) +
  ggtitle("FVS Carbon Projection by FIA Subplot") +
  ylab(bquote("Carbon " ~(Mg %*% ha^{-1}))) +
  theme(legend.title = element_blank()) +
  geom_line() +
  coord_cartesian(xlim = c(2005, 2165), ylim = c(0, 320)) +
  scale_x_continuous(breaks=seq(2005,2165,20)) +
  scale_y_continuous(breaks=seq(0,320,20))
```

Examine RMSE for the residuals:

```{r}
Subplots_RMSE <- data.frame(
  Total_Carbon = sqrt(mean(
    (NoManagement_Carbon_BySubplot_ByYear$Total_Carbon - NoManagement_Carbon_BySubplot_ByYear$NK_Fig2_NoManagement)^2
  )),
  Aboveground_Carbon = sqrt(mean(
    (NoManagement_Carbon_BySubplot_ByYear$Aboveground_Carbon - NoManagement_Carbon_BySubplot_ByYear$NK_Fig2_NoManagement)^2
  )),
  Aboveground_Live = sqrt(mean(
    (NoManagement_Carbon_BySubplot_ByYear$Aboveground_Live - NoManagement_Carbon_BySubplot_ByYear$NK_Fig2_NoManagement)^2
  ))
)
ggplot(
    data = melt(Subplots_RMSE, id = NULL, variable.name = "Series", value.name = "RMSE"),
    mapping = aes(x = Series, y = RMSE)
  ) +
  ggtitle("FVS Carbon Projection by FIA Subplot") +
  geom_col()
```

Note that all RMSE for the subplot projections are significantly worse
than for the plot projections.

### By Condition

We repeat the exercise using the FVS_StandInit_Cond table, which creates
FVS stands for corresponding FIA conditions.

We select all conditions using the translated stand IDs. Note that plot
360304303966 has 3 conditions; we include all three. We apply the same
configuration to all stands as we did for FVS_StandInit_Plot. We load
and clean the results in the same manner.

```{r}
NoManagement_Carbon_ByCondition_ByYear <- read_csv(
    "data/fvs/FVS_NKByCondition_NONE_Carbon.csv",
    col_types = cols(StandID = col_character())
  ) |>
  rename(Total_Carbon = Total_Stand_Carbon) |>
  rename(Aboveground_Live = Aboveground_Total_Live) |>
  rename(Down_Dead = Forest_Down_Dead_Wood) |>
  mutate(Aboveground_Dead = Standing_Dead + Down_Dead) |>
  mutate(Aboveground_Carbon = Aboveground_Live + Aboveground_Dead) |>
  select(StandID, Year, Total_Carbon, Aboveground_Carbon,
         Aboveground_Live, Aboveground_Dead, Standing_Dead, Down_Dead) |>
  group_by(Year) |>
  summarize(
    Total_Carbon=mean(Total_Carbon),
    Aboveground_Carbon=mean(Aboveground_Carbon),
    Aboveground_Live=mean(Aboveground_Live),
    .groups = "keep"
  ) |>
  filter(Year >= 2005) |>
  full_join(fig2 |> select(Year,NoManagement), by=join_by(Year)) |>
  rename(NK_Fig2_NoManagement = NoManagement)

ggplot(
    data = melt(NoManagement_Carbon_ByCondition_ByYear, id.vars = "Year"),
    mapping = aes(x = Year, y = value, color = variable)
  ) +
  ggtitle("FVS Carbon Projection by FIA Condition") +
  ylab(bquote("Carbon " ~(Mg %*% ha^{-1}))) +
  theme(legend.title = element_blank()) +
  geom_line() +
  coord_cartesian(xlim = c(2005, 2165), ylim = c(0, 320)) +
  scale_x_continuous(breaks=seq(2005,2165,20)) +
  scale_y_continuous(breaks=seq(0,320,20))
```

Examine RMSE for the residuals:

```{r}
Conditions_RMSE <- data.frame(
  Total_Carbon = sqrt(mean(
    (NoManagement_Carbon_ByCondition_ByYear$Total_Carbon - NoManagement_Carbon_ByCondition_ByYear$NK_Fig2_NoManagement)^2
  )),
  Aboveground_Carbon = sqrt(mean(
    (NoManagement_Carbon_ByCondition_ByYear$Aboveground_Carbon - NoManagement_Carbon_ByCondition_ByYear$NK_Fig2_NoManagement)^2
  )),
  Aboveground_Live = sqrt(mean(
    (NoManagement_Carbon_ByCondition_ByYear$Aboveground_Live - NoManagement_Carbon_ByCondition_ByYear$NK_Fig2_NoManagement)^2
  ))
)
ggplot(
    data = melt(Conditions_RMSE, id = NULL, variable.name = "Series", value.name = "RMSE"),
    mapping = aes(x = Series, y = RMSE)
  ) +
  ggtitle("FVS Carbon Projection by FIA Condition") +
  geom_col()
```

RMSE for Aboveground Live by Condition is \~13, which is close to that
for Aboveground Live by Plot.

```{r}
Conditions_RMSE |> select(Aboveground_Live)
```