Harjoitus 4 For-255 Metsien kehityksen simulointi

Mustikkasatomalli Kurttila et al. 2018

Saatavilla: https://www.mdpi.com/1999-4907/9/7/417

rm(list=ls())
library(lmfor)
data("spati")

Kopioidaan Kurttilan artikkelin mustikkasatomalli, joka ennustaa vuosittaista mustikkasadon määrää kuviolla (kg/ha) D = rinnankorkeuden keskiläpimitta (cm) P_koivu = koivun osuus kuviolla P_kuusi = kuusen osuus kuviolla G = kuvion pohjapinta-ala OMT, VT, CT kasvupaikat, dummy-muuttujia

Malli luultavasti ennustaa mäntyvaltaiselle metsälle mustikkasadon MT kasvupaikalle perusoletuksena (kun kasvupaikat ja puulajien osuudet saavat arvon nolla)

Verrataan kasvupaikan vaikutusta mustikkasatoon vertaamalla kolmea eri mallinnusta eri kasvupaikoille vaihtelemalla kasvupaikkaa puuston pysyessä samana (MT, VT ja CT (mänty ei sovellu OMT, pelkälle mäntyvaltaiselle puustolle, koska data sisältää vain mäntyä)

malli <- function(D,P_koivu,P_kuusi,G,OMT,VT,CT){exp(1.723+0.376*sqrt(D)-0.543*P_koivu-0.288*P_kuusi-0.083*G/log(D+1)+0.000519*D*G-0.000263*G^2-1.542*OMT-0.626*VT-1.511*CT)}

Kasvupaikat MT, VT ja CT

ennusteet_MT <- c()
ennusteet_VT <- c()
ennusteet_CT <- c()
for (i in 1:66) {
  kuvio <- subset(spati, plot==i)
  tulos_MT <- malli(mean(kuvio$d),0,0,mean(kuvio$G),0,0,0)
  tulos_VT <- malli(mean(kuvio$d),0,0,mean(kuvio$G),0,1,0)
  tulos_CT <- malli(mean(kuvio$d),0,0,mean(kuvio$G),0,0,1)
  ennusteet_MT <- c(ennusteet_MT, tulos_MT)
  ennusteet_VT <- c(ennusteet_VT, tulos_VT)
  ennusteet_CT <- c(ennusteet_CT, tulos_CT)
}

Luodaan kuvio

kuviot <- c(1:66)
plot(c(0,66), c(0,34), type="n", xlab = "Kuvion numero", ylab = "Mustikkasato kg/ha")
points(kuviot, ennusteet_MT, col=2)
points(kuviot, ennusteet_VT, col=3)
points(kuviot, ennusteet_CT, col=4)
legend("topright", c("MT", "VT", "CT"), lty=c(1), col=c(2,3,4))

Simulointi

Simuloidaan kolme eri käsittelyvaihtoehtoa:

Kasvatus 5 vuotta
Alaharvennus 5 vuoden päästä
Yläharvennus 5 vuoden päästä

Valmiiksi annettu koodi

# Kaytto-ohje:
# Lataa Moodlesta tama ja erillinen, simulaattorin sisalta liitetiedosto.
# - Ilmoita alla liitetiedoston hakemistopolku.
# - Vaihda halutessasi simuloitavaa kuviota esim.kuvio-kohdassa
# - Aja koodi ja tutustu loppuosassa simulaattorin tuloksiin (halutessasi itse simulaattoriin)


#simulaattorin.sijainti <- "sijainti, jossa simulaattori sijaitsee, piilotettu"

#### KUVIOTIEDOT

# Poimitaan esimerkiksi spati-datan rivi numero 1; sen kuviotunnukset ppa, runkoluku ja keskipuun tunnukset
esim.kuvio <- spati[1, c(1,4,5,7,8)]

# Kuviotietorakenne - voidaan olettaa, etta kaikki kuviot ovat puhtaita mannikoita, kivennais-VT:lla ja niiden lamposumma on 1100
standdata <- list('peat'=0, 'sitetype'="VT", 
                  'TS'=1100,
                  'Dgmean'=esim.kuvio$Dg,
                  'Hgmean'=esim.kuvio$Hg,
                  'Gsp'=c(esim.kuvio$G, 0, 0),
                  'Nsp'=c(esim.kuvio$N, 0, 0),
                  'Dsp'=c(esim.kuvio$Dg, 0, 0),
                  'plot'=esim.kuvio$plot)
#### SIMULOINTI
# Ladataan itse tehty simulaattori eri tiedostosta
source(simulaattorin.sijainti)
#### SIMULOINNIN TULOKSET

# Simuloitiin kolme kasittelyvaihtoehtoa:
# 1. Lepo
# 2. Alaharvennus
# 3. Ylaharvennus

# Piirra eri vaihtoehtojen runkolukusarjat

plot(dclass.mean, Reduce('+',d.freqs), type="l", xlab="Lpm, cm", ylab="Runkoja, 1/ha")
lines(dclass.mean, Reduce('+',d.freqs.5), col="grey")
lines(dclass.mean, Reduce('+',d.freqs.below.5), col="grey", lty=2, lwd=1)
lines(dclass.mean, Reduce('+',d.freqs.above.5), col="grey", lty=3, lwd=2)
legend("topright", c("Alkup. mallinnettu", "Kasvatettu 5 v", "Alaharvennettu, kasv. 5v", "Ylaharvennettu, kasv. 5v"), lty=c(1,1,2,3), lwd=c(1,1,1,2), col=c(1, rep("grey",3)))

# Tulosta tietoa vaihtoehdoista (ks. vaihtoehdon numerot ylta)
tulosta(vaihtoehto=1) # lepo
tulosta(vaihtoehto=2) # alaharvennus
tulosta(vaihtoehto=3) # ylaharvennus

# Lisaa mahdollisesti hyodyllista tietoa (keskitunnuksia)
stand.attr # kasittelematon puusto
stand.attr.5 # 5 vuotta kasvatettu
stand.attr.5.thin.below # alaharvennettu 5 v paasta
stand.attr.5.thin.above # ylaharvennettu 5 v paasta

Tehdään simulointi omalle aineistolle simuloimalla kolme eri käsittelyvaihtoehtoa jokaiselle kuviolle

#luodaan tyhjät vektorit joihin tallennetaan simuloinnin tulokset
kasvatettu_Dmean <- c()
kasvatettu_G <- c()
kasvatettu_plot <- c()
alaharvennettu_Dmean <- c()
alaharvennettu_G <- c()
alaharvennettu_plot <- c()
yläharvennettu_Dmean <- c()
yläharvennettu_G <- c()
yläharvennettu_plot <- c()

for (i in 1:66) {
  otetaan_kuvio <- subset(spati, plot==i) #valitaan simulointiin jokainen kuvio yksitellen
  esim.kuvio <- otetaan_kuvio[1, c(1,4,5,7,8)] #valitaan kuvion ensimmäinen rivi ja tältä halutut tunnukset (edustaa koko kuviota)
  standdata <- list('peat'=0, 'sitetype'="VT", 
                    'TS'=1100, 
                    'Dgmean'=esim.kuvio$Dg,
                    'Hgmean'=esim.kuvio$Hg,
                    'Gsp'=c(esim.kuvio$G, 0, 0),
                    'Nsp'=c(esim.kuvio$N, 0, 0),
                    'Dsp'=c(esim.kuvio$Dg, 0, 0),
                    'plot'=esim.kuvio$plot)
  source(simulaattorin.sijainti)
  #tallennetaan simuloinnin tulokset
  kasvatettu_Dmean <- c(kasvatettu_Dmean, stand.attr.5$Dmean)
  kasvatettu_G <- c(kasvatettu_G, stand.attr.5$G)
  kasvatettu_plot <- c(kasvatettu_plot, stand.attr.5$plot)
  alaharvennettu_Dmean <- c(alaharvennettu_Dmean, stand.attr.5.thin.below$Dmean)
  alaharvennettu_G <- c(alaharvennettu_G, stand.attr.5.thin.below$G)
  alaharvennettu_plot <- c(alaharvennettu_plot, stand.attr.5.thin.below$plot)
  yläharvennettu_Dmean <- c(yläharvennettu_Dmean, stand.attr.5.thin.above$Dmean)
  yläharvennettu_G <- c(yläharvennettu_G, stand.attr.5.thin.above$G)
  yläharvennettu_plot <- c(yläharvennettu_plot, stand.attr.5.thin.above$plot)
}

## Error in if (d.freqs[[sp.cur]][dcl] > 0) {: missing value where TRUE/FALSE needed

Aiheuttaa “errorin” jostain syystä, mutta ei luultavasti vaikuta tuloksiin (syntyy tuloksia).

Luodaan taulukko simuloinnin tuloksista, josta voidaan poimia myöhemmin tunnuksia mustikkasadon laskentaa varten

taulukko_kasvatus <- data.frame(kasvatettu_plot,kasvatettu_Dmean,kasvatettu_G)
taulukko_alaharvennus <- data.frame(alaharvennettu_plot,alaharvennettu_Dmean,alaharvennettu_G)
taulukko_yläharvennus <- data.frame(yläharvennettu_plot,yläharvennettu_Dmean,yläharvennettu_G)

Tulostetaan taulukko_kasvatus esimerkiksi:

taulukko_kasvatus

##    kasvatettu_plot kasvatettu_Dmean kasvatettu_G
## 1                1         16.42079     22.79192
## 2                2         18.50816     21.02316
## 3                3         13.28003     21.52153
## 4                4         12.53742     24.57062
## 5                5         16.95586     22.63581
## 6                6         20.67838     23.34044
## 7                7         23.20410     20.96215
## 8                8         27.81783     27.03608
## 9                9         25.79143     29.31916
## 10              10         21.06195     28.64678
## 11              11         21.49750     22.23984
## 12              12         24.64057     29.19826
## 13              13         25.67161     30.67827
## 14              14         14.52161     24.81266
## 15              15         17.65818     24.80606
## 16              16         12.20025     23.53305
## 17              17         17.07585     18.89825
## 18              18         26.68441     23.58867
## 19              19         12.98066     20.73746
## 20              20         11.96416     20.88429
## 21              21         22.99741     24.66303

Jostain syystä simulointi kohdistuu vain 21 ensimmäiseen kuvioon (en saanut syytä selville).

Mallinnetaan mustikkasadot kasvatetulle ja käsitellyille metsiköille ja verrataan alkutilanteeseen. Tarkastellaan vain VT kasvupaikan osalta (simulointi oletti kasvupaikan olevan VT)

#luodaan tyhjät vektorit, joihin tallennetaan saadut ennusteet mustikkasadolle
ennusteet_kasvatus <- c()
ennusteet_alaharvennus <- c()
ennusteet_yläharvennus <- c()
for (i in 1:21) {
  kuvio1 <- subset(taulukko_kasvatus, kasvatettu_plot==i) #valitaan kuviot yksitellen
  kuvio2 <- subset(taulukko_alaharvennus, alaharvennettu_plot==i)
  kuvio3 <- subset(taulukko_yläharvennus, yläharvennettu_plot==i)
  tulos_kasvatus <- malli(mean(kuvio1$kasvatettu_Dmean),0,0,mean(kuvio1$kasvatettu_G),0,1,0) #ennustetaan mustikkasato kuvioittain kullekin käsittelyvaihtoehdolle erikseen
  tulos_alaharvennus <- malli(mean(kuvio2$alaharvennettu_Dmean),0,0,mean(kuvio2$alaharvennettu_G),0,1,0)
  tulos_yläharvennus <- malli(mean(kuvio3$yläharvennettu_Dmean),0,0,mean(kuvio3$yläharvennettu_G),0,1,0)
  ennusteet_kasvatus <- c(ennusteet_kasvatus, tulos_kasvatus) #tallennetaan tulokset omiin vektoreihin
  ennusteet_alaharvennus <- c(ennusteet_alaharvennus, tulos_alaharvennus)
  ennusteet_yläharvennus <- c(ennusteet_yläharvennus, tulos_yläharvennus)
}

Luodaan kuvio, joka havainnollistaa käsittelyvaihtoehtojen vaikutusta mustikkasatoon

plots <- c(1:21)
plot(c(0,21), c(0,30), type="n", xlab = "Kuvion numero", ylab = "Mustikkasato kg/ha")
points(kuviot, ennusteet_VT, col=3)
points(plots, ennusteet_kasvatus, col=4)
points(plots, ennusteet_alaharvennus, col=5)
points(plots, ennusteet_yläharvennus, col=6)
legend("topleft", c("Alkup. mallinnettu", "Kasvatettu 5v", "Alaharvennettu 5v päästä", "Yläharvennettu 5v päästä"), lty=c(1), col=c(3,4,5,6))

Harjoitus 4 For-255 Metsien kehityksen simulointi

Johannes

12/13/2021

Mustikkasatomalli Kurttila et al. 2018

Simulointi