Forestat
forestat: Forest Carbon Sequestration and Potential Productivity Calculation 森林碳汇计量和潜力计算
Install / Use
/learn @caf-ifrit/ForestatREADME
<div align="center"><strong>森林碳汇计量和潜力计算</strong></div>
<p align="right"><strong>Forestat version:</strong> 1.1.0</p> <p align="right"><strong>Date:</strong> 10/10/2023 </p> <br>forestat 是基于中国林业科学研究院资源信息研究所（Institute of Forest Resource Information Techniques, Chinese Academy of Forestry）提出的基于林分潜在生长量的立地质量评价方法与应用<sup>[1]</sup>和A basal area increment-based approach of site productivity evaluation for multi-aged and mixed forests<sup>[2]</sup>开发的R包。可依据林分高生长，划分立地等级；并以全林整体模型为基础，建立不同立地等级下的非线性混合效应生物量模型，实现更精准的碳汇计量；尤其提出了一种基于林分潜在生长量的碳汇潜力计算方法。该套算法适用于天然林和人工林，能够定量回答一定立地条件下的潜在生产力、现实生产力、提升空间有多大，可用于立地质量评价、树种适宜性评价、退化林评价等多个方面。
<div align="center">1 概述</div>
forestat 包实现了碳汇潜力计算和退化林评价，其中碳汇潜力计算包括了基于林分高生长的立地等级划分，树高模型、断面积生长模型、生物量生长模型的建立，林分现实生产力与潜在生产力的计算。树高模型可用Richard模型、Logistic模型、korf模型、Gompertz模型、Weibull模型和Schumacher模型构建，断面积生长模型和生物量生长模型仅可用Richard模型构建。碳汇潜力计算依赖于给定林分类型（树种）的多个样地数据，退化林的评价依赖于多个样木和样地数据，forestat 包中带有一些样例数据。
1.1 forestat 流程图
<div align="center"> <img width="70%" src="forestat/vignettes/img/flowchart-1.png"> <p>图 1.1 碳汇计算工作流程图</p> </div> <div align="center"> <img width="50%" src="forestat/vignettes/img/flowchart-2.png"> <p>图 1.2 退化林评价工作流程图</p> </div>1.2 forestat 依赖的R包
| Package | Download Link | | ----------- | ------------------------------------------ | | dplyr | https://CRAN.R-project.org/package=dplyr | | ggplot2 | https://CRAN.R-project.org/package=ggplot2 | | nlme | https://CRAN.R-project.org/package=nlme |
<div align="center">2 安装</div>
2.1 从CRAN或GitHub安装
在 R 中使用以下命令从 CRAN 安装 forestat ：
# 安装forestat
install.packages("forestat")
当然，你也可以在 R 中使用以下命令从 GitHub 安装 forestat ：
# 安装devtools
install.packages("devtools")
# 安装forestat
devtools::install_github("caf-ifrit/forestat/forestat")
2.2 加载 forestat
library(forestat)
<div align="center">3 快速开始</div>
本部分展示的是快速完成立地分级、潜在生产力和现实生产力的完整步骤，使用的数据是包中自带的forestData样例数据。
# 加载包中 forestData 样例数据
data("forestData")
# 基于 forestData 数据建立模型，返回一个 forestData 类对象
forestData <- class.plot(forestData, model = "Richards",
interval = 5, number = 5,
H_start=c(a=20,b=0.05,c=1.0))
# 绘制树高生长模型散点图
plot(forestData,model.type = "H",plot.type = "Scatter",
title = "桦木阔叶树高生长模型散点图")
# 计算 forestData 对象的潜在生产力
forestData <- potential.productivity(forestData)
# 计算 forestData 对象的现实生产力
forestData <- realized.productivity(forestData)
# 获取 forestData 对象的汇总数据
summary(forestData)
本部分展示的是快速完成退化林评价的完整步骤，使用的数据是包中自带的tree_1、tree_2、 tree_3、plot_1、plot_2、plot_3样例数据。
# 加载包中tree_1 tree_2 tree_3 plot_1 plot_2 plot_3样例数据
data(tree_1)
data(tree_2)
data(tree_3)
data(plot_1)
data(plot_2)
data(plot_3)
# 对退化林数据进行预处理
plot_data <- degraded_forest_preprocess(tree_1, tree_2, tree_3,
plot_1, plot_2, plot_3)
# 计算退化林评价指标
res_data <- calc_degraded_forest_grade(plot_data)
# 查看计算结果
res_data
<div align="center">4 碳汇潜力计算</div>
<details> <summary style="font-size:21px;"><strong>4.1 建立模型</strong></summary> <br> <details> <summary style="font-size:18px;"><strong>4.1.1 自定义数据</strong></summary>为了建立一个准确的模型，好的数据是不可或缺的，在 forestat 包中内置了一个经过清洗的样例数据，可以通过如下命令加载查看样例数据：
# 加载包中 forestData 样例数据
data("forestData")
# 筛选 forestData 样例数据中ID、code、AGE、H、S、BA 和 Bio字段
# 并查看前6行数据
head(dplyr::select(forestData,ID,code,AGE,H,S,BA,Bio))
# 输出
ID code AGE H S BA Bio
1 1 1 13 2.0 152.67461 4.899382 32.671551
2 2 1 15 3.5 68.23825 1.387268 5.698105
3 3 1 20 4.2 128.32683 3.388492 22.631467
4 4 1 19 4.2 204.93928 4.375324 18.913886
5 5 1 13 4.2 95.69713 1.904063 6.511951
6 6 1 25 4.7 153.69393 4.129810 28.024739
当然，你也可以选择加载自定义数据：
# 加载自定义数据
forestData <- read.csv("/path/to/your/folder/your_file.csv")
自定义数据中ID（样地ID）、code（样地林分类型代码）、AGE（林分平均年龄）、H（林分平均高）是必须字段，用以建立树高模型（H-model），并绘制相关示例图。
S（林分密度指数）、BA（林分断面积）、Bio（林分生物量）是可选的字段，用以建立断面积生长模型（BA-model）与生物量生长模型（Bio-model）。
在后续的潜在生产力和现实生产力计算中，断面积生长模型与生物量生长模型是必须的。也就是自定义数据如果缺少S、BA和Bio字段将无法计算潜在生产力和现实生产力。
数据加载后，forestat 将使用class.plot()函数构建林分生长模型，如果自定义数据中同时包含ID、code、AGE、H、S、BA、Bio字段，则会同时构建树高模型、断面积生长模型、生物量生长模型，如果只包含ID、code、AGE、H字段，则只会构建树高模型。
# 选用 Richards 模型构建林分生长模型
# interval=5表示初始树高分类的林分年龄区间设置为5，number=5表示初始树高分类数最多为5，maxiter=1000表示拟合模型的最大次数为1000
# 树高模型拟合的初始参数H_start默认为c(a=20,b=0.05,c=1.0)
# 断面积生长模型拟合的初始参数BA_start默认为c(a=80, b=0.0001, c=8, d=0.1)
# 生物量生长模型拟合的初始参数Bio_start默认为c(a=450, b=0.0001, c=12, d=0.1)
forestData <- class.plot(forestData, model = "Richards",
interval = 5, number = 5, maxiter=1000,
H_start=c(a=20,b=0.05,c=1.0),
BA_start = c(a=80, b=0.0001, c=8, d=0.1),
Bio_start=c(a=450, b=0.0001, c=12, d=0.1))
其中，model为构建树高模型时选用的模型，可在"Logistic"、"Richards"、"Korf"、"Gompertz"、"Weibull"、"Schumacher"模型中任选一个，断面积生长模型和生物量生长模型默认使用Richard模型构建。interval为初始树高分类的林分年龄区间，number为初始树高分类数的最大值，maxiter为最大拟合次数。H_start 是拟合树高模型的初始参数，BA_start 是拟合断面积生长模型的初始参数，H_start 是拟合生物量生长模型的初始参数，当拟合出现错误时，可以多尝试一些初始参数作为尝试。
由class.plot()函数返回的结果为forestData 对象，包括Input（输入数据和树高分级结果）、Hmodel（树高模型结果）、BAmodel（断面积模型结果）、Biomodel（生物量模型结果）以及output（所有模型的表达式、参数及精度）。
为了解模型的建立情况，可以使用summary(forestData)函数获取forestData对象汇总数据。该函数返回summary.forestData对象并将相关数据输出至屏幕。
输出的第一段为输入数据的汇总，第二、三、四段分别为H-model、BA-model、Bio-model的参数及其精简报告。
summary(forestData)
# 输出
# 第一段
H S BA Bio
Min. : 2.00 Min. : 68.24 Min. : 1.387 Min. : 5.698
1st Qu.: 8.10 1st Qu.: 366.37 1st Qu.: 9.641 1st Qu.: 52.326
Median :10.30 Median : 494.76 Median :13.667 Median : 78.502
Mean :10.62 Mean : 522.53 Mean :14.516 Mean : 90.229
3rd Qu.:12.90 3rd Qu.: 661.84 3rd Qu.:18.750 3rd Qu.:115.636
Max. :19.10 Max. :1540.13 Max. :45.749 Max. :344.412
# 第二段
H-model Parameters:
Nonlinear mixed-effects model fit by maximum likelihood
Model: H ~ 1.3 + a * (1 - exp(-b * AGE))^c
Data: data
AIC BIC logLik
728.4366 747.2782 -359.2183
Random effects:
Formula: a ~ 1 | LASTGROUP
a Residual
StdDev: 3.767163 0.7035752
Fixed effects: a + b + c ~ 1
Value Std.Error DF t-value p-value
a 12.185054 1.7050081 313 7.146625 0
b 0.037840 0.0043682 313 8.662536 0
c 0.761367 0.0769441 313 9.895060 0
Correlation:
a b
b -0.110
c -0.093 0.946
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-3.858592084 -0.719253472 0.007120413 0.761123585 3.375793806
Number of Observations: 320
Number of Groups: 5
Concise Parameter Report:
Model Coefficients:
a1 a2 a3 a4 a5 b c
7.013778 9.575677 11.90324 14.67456 17.75801 0.03783956 0.7613666
Model Evaluations:
pe RMSE R2 Var TRE AIC BIC logLik
-0.006484677 0.6980625 0.9543312 0.4887767 0.3960163 728.4366 747.2782 -359.2183
Model Formulas:
Func Spe
model1:H ~ 1.3 + a * (1 - exp(-b * AGE))^c model1:pdDiag(a ~ 1)
# 第三段（与第二段数据格式相似）
BA-model Parameters:
# 此处省略
......
# 第四段（与第二段数据格式相似）
Bio-model Parameters:
# 此处省略
......
经过4.1.2 class.plot()函数构建林分生长模型后，就可以使用plot()函数绘制图像。
其中，model.type为绘图使用的模型，可以选择H（树高模型）、BA（断面积生长模型）或者Bio（生物量生长模型）。plot.type为绘图的类型，可以选择Curve（曲线图）、Residual（残差图）、Scatter_Curve（散点曲线图）、Scatter（散点图）。xlab、ylab、legend.lab、title分别为x轴标题、y轴标题、图例、图像标题。
# 绘制树高模型的曲线图
plot(forestData,model.type="H",
plot.type="Curve",
xlab="年龄(year)",ylab="树高(m)",legend.lab="立地等级",
title="桦木阔叶混树高模型曲线图")
# 绘制断面积生长模型散点图
plot(forestData,model.type="BA",
plot.type="Scatter",
xlab="年龄(year)",ylab=expression(paste("断面积( ",m^2,"/",hm^2,")")),legend.lab="立地等级",
title="桦木阔叶混断面积生长模型散点图")
不同的plot.type绘制的样图如图4所示：
经过4.1.2 class.plot()函数构建林分生长模型后，就可以使用potential.productivity()函数计算林分潜在生产力。在计算之前，要求forestData 对象中BA-model和Bio-model已经建立。
forestData <- potential.productivity(forestData, code=1,
age.min=5,age.max=150,
left=0.05, right=100,
e=1e-05, maxiter = 50)
其中，参数code为计算潜在生产力使用的林分类型代码。age.min和age.max分别为林分年龄的最小值和最大值，潜在生产力的计算会在最小值和最大值的区间中进行。left和right为拟合模型的初始参数，当拟合出现错误时，可以多尝试一些初始参数作为尝试。e为拟合模型的精度，当残差低于e时，认为模型收敛并停止拟合。maxiter为拟合模型的最大次数，当拟合次数等于maxiter时，认为模型收敛并停止拟合。
计算结束后，可以使用如下命令查看并输出结果：
library(dplyr)
forestData$potential.productivity
