随机森林建模及反演流程（遥感影像）

随机森林建模及反演流程（遥感影像）

如果需要进遥感深度学习微信交流群，或者需要代码请添加作者微信（ictnet1）无偿领取/邀您进群！

我把上一步导出来的表格放在了 E:\xianyu\landsat反演盖度和可食牧草产量\jianmo\样本\图斑样本

路径下。

下面我们搜索anaconda prompt打开，并输入E:(这里我数据放在E盘下所以输入E,如果你是放在F盘那就是输入F：)意思是进入E盘

然后采用cd命令进入我们存放样本和代码的路径下，输入jupyter notebook打开python的编辑器，它会自动弹出，以网页的形式，注意不要关闭这个黑窗。

如何运行？

下面是对代码的详细介绍

1.数据导入与准备

首先，我们需要导入必要的库，并读取数据文件。以下是详细步骤：

import pandas as pd  # 导入用于数据处理的pandas库
import numpy as np  # 导入用于数值计算的numpy库
from matplotlib import pyplot as plt  # 导入用于绘图的matplotlib库
from sklearn.ensemble import RandomForestRegressor  # 导入随机森林回归模型
from sklearn.model_selection import train_test_split  # 导入数据集拆分工具
from sklearn.metrics import mean_squared_error, r2_score  # 导入误差计算和R²评分工具
from sklearn import metrics  # 导入额外的度量工具

读取表格数据

我们有一个CSV文件new_yb0.csv，其中第一列是目标变量（即我们要预测反演的变量），其余列是特征变量（即用于预测的指标）。我们使用pandas库读取该文件：

# 读取表格数据, 将excel另存为csv文件
data = pd.read_csv(r'./new_yb0.csv', na_values=[' '])  # 读取CSV文件，空值处理为NaN
y = data.iloc[:, 0].values  # 提取第一列作为目标变量
X = data.iloc[:, 1:].values  # 提取其他列作为特征变量
#X = data[['dem', 'b4_ronghe', 'gcvi', 'gndvi', 'isr2']].values  
df = pd.DataFrame(data)  # 将数据转换为DataFrame格式

分割数据集

为了评估模型的性能，我们将数据集拆分为训练集和测试集：

random_forest_seed = np.random.randint(low=1, high=230)  # 生成随机种子
# 分割数据集为训练集和测试集，比例为7:3
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=random_forest_seed)

2.构建并训练随机森林模型

接下来，我们构建一个随机森林回归模型，并对训练集数据进行训练：

# 初始化随机森林回归器，设置树的数量为600，最大特征数为8，这个您可以自行改动
regressor = RandomForestRegressor(n_estimators=600, max_features=8)
regressor.fit(X_train, y_train)  # 使用训练数据拟合模型

模型误差分析

我们可以通过绘制误差分布图来观察模型的误差情况，如果绘制的图像呈现正态分布，说明模型很好

# 预测测试集数据
random_forest_predict = regressor.predict(X_test)
# 计算预测误差
random_forest_error = random_forest_predict - y_test
# 绘制误差分布直方图
plt.figure(1)
plt.clf()
plt.hist(random_forest_error)
plt.xlabel('Prediction Error')
plt.ylabel('Count')
plt.grid(False)
plt.show()

精度验证

通过计算R²值和均方根误差（RMSE）来评估模型精度：

random_forest_R2 = metrics.r2_score(y_test, random_forest_predict)
random_forest_RMSE = metrics.mean_squared_error(y_test, random_forest_predict)**0.5
print('R^2 = {0} and RMSE = {1}.'.format(random_forest_R2, random_forest_RMSE))

3.变量重要性分析

随机森林模型能够评估每个特征在预测中的重要性。我们可以绘制重要性图来直观地显示前20个最重要的变量：

# 计算变量重要性
random_forest_importance = list(regressor.feature_importances_)
# 选择前20个重要变量
random_forest_feature_importance = [(feature, round(importance, 9)) 
                                    for feature, importance in zip(df.columns[1:], random_forest_importance)]
# 按重要性排序
random_forest_feature_importance = sorted(random_forest_feature_importance, key=lambda x:x[1], reverse=True)

# 提取变量名和重要性
x_data = [item[0] for item in random_forest_feature_importance[:20]]
y_data = [item[1] for item in random_forest_feature_importance[:20]]

# 设置图表大小
plt.figure(figsize=(16, 10))
# 绘制柱状图
plt.bar(x_data, y_data)
# 设置横轴标签
plt.xlabel('Variable', fontsize=16)
plt.xticks(rotation=-90, fontsize=16)
# 设置纵轴标签
plt.ylabel('Importance', fontsize=20)
# 设置标题
plt.title('Variable Importances', fontsize=16)
# 显示图表
plt.show()

4.栅格数据预测

最后，我们将使用训练好的随机森林模型对栅格数据进行预测。这里我们根据变量重要性分析，重新选择了前5个重要的变量，把代码X = data.iloc[:, 1:].values 改为 X = data[[‘dem’, ‘b4_ronghe’, ‘gcvi’, ‘gndvi’, ‘isr2’]].values （然后从头开始运行代码），读取这5个变量即可实现研究区反演，我们将对这些数据进行反演预测：

import rasterio  # 导入栅格数据处理库
from tqdm import tqdm  # 导入进度条显示库

# 读取栅格数据
with rasterio.open(r'./dem.tif') as src:
    data1 = src.read(1)
    meta = src.meta
with rasterio.open(r'./b4_ronghe.tif') as src:
    data2 = src.read(1)
with rasterio.open(r'./gcvi.tif') as src:
    data3 = src.read(1)
with rasterio.open(r'./gndvi.tif') as src:
    data4 = src.read(1)
with rasterio.open(r'./isr2.tif') as src:
    data5 = src.read(1)

# 整合栅格数据为一个矩阵
X = np.stack((data1, data2, data3, data4, data5), axis=-1)

# 清洗输入数据
X_2d = X.reshape(-1, X.shape[-1])
print(np.isnan(X_2d).any())  # 检查数据中是否存在NaN值

# 将NaN值替换为0
X_2d[np.isnan(X_2d)] = 0

# 使用训练好的模型进行预测
y_pred = []
for i in tqdm(range(0, X_2d.shape[0], 10000)):
    y_pred_chunk = regressor.predict(X_2d[i:i+10000])
    y_pred.append(y_pred_chunk)
y_pred = np.concatenate(y_pred)

# 保存预测结果为新的栅格数据
with rasterio.open(r'./建模变量.tif', 'w', **meta) as dst:
    dst.write(y_pred.reshape(X.shape[:-1]), 1)
print("预测结束")

运行过程会产生一个进度条，进度条走完代表反演结束

欢迎关注我们的公众号：Python与遥感，或许里面有你想要的免费教程！