【maxENT】最大熵模型（Maximum Entropy Model）R语言实现

🟢🟠先看：【maxENT】最大熵模型（Maximum Entropy Model）介绍与使用（maxENT软件）

ASCII文件太大，使用R语言直接输入环境的栅格文件。

一、相关package介绍

1.1 dismo 包

专门用于物种分布建模（Species Distribution Modeling, SDM），提供与MaxEnt模型的接口，支持数据预处理、模型训练、预测及评估。

核心函数详解
1. maxent()
功能: 训练MaxEnt模型（基于Java实现）。
参数:

maxent(x, p, removeDuplicates=TRUE, path, args, ...)

• x

  • 类型: RasterStack 或 RasterBrick
  • 作用: 包含环境变量（如温度、海拔）的多层栅格数据，作为模型输入。

• p

  • 类型: 两列数据框（data.frame）
  • 作用: 物种存在点的坐标，列名需为 x（经度）和 y（纬度）。
  • 示例:

    presence_points <- data.frame(x = c(100, 101), y = c(30, 31))
    

• removeDuplicates

  • 类型: 逻辑值（TRUE/FALSE）
  • 默认: TRUE
  • 作用: 是否移除重复坐标点，避免过拟合。

• path

  • 类型: 字符串
  • 作用: 保存MaxEnt输出文件的目录路径（包括HTML报告、模型参数、预测结果图）。
  • 示例: path = "G:/sdm_output"

• args

  • 类型: 字符向量
  • 作用: 控制MaxEnt模型行为的命令行参数，常用选项：
    • betamultiplier: 正则化系数（默认1），增大可简化模型（防过拟合）。
    • linear/quadratic/product/hinge/threshold: 启用（TRUE）或禁用（FALSE）特定特征类型。
      • linear: 线性特征（适合单调响应）
      • hinge: 铰链特征（适合非线性响应）
    • threads: 使用的CPU线程数（加速计算，如 threads=12）。
    • responsecurves: 是否生成响应曲线（TRUE/FALSE）。
  • 示例:

    args = c("betamultiplier=2", "hinge=FALSE", "threads=8")
    

2. predict()
功能: 应用训练好的模型预测物种分布概率。
参数:

predict(object, x, ext=NULL, filename="", progress='text', ...)

• object

  • 类型: MaxEnt 模型对象
  • 作用: 训练好的MaxEnt模型。

• x

  • 类型: RasterStack
  • 作用: 输入的环境变量栅格数据（需与训练数据层数、名称一致）。

• filename

  • 类型: 字符串
  • 作用: 直接保存预测结果到指定路径（支持.tif、.grd格式）。
  • 示例: filename = "G:/prediction.tif"

• progress

  • 类型: 字符串
  • 作用: 显示进度条（'text' 或 'window'）。

1.2 raster包

处理栅格数据（如地理空间栅格图像），支持数据读写、裁剪、计算、统计及可视化。

核心函数详解
1. rasterOptions()
功能: 全局配置栅格数据处理参数。
常用参数:

rasterOptions(maxmemory=1e+09, tmpdir="path/to/temp", chunksize=1e+08)

• maxmemory

  • 类型: 数值
  • 作用: 设置R可使用的最大内存量（单位：字节）。
  • 示例: maxmemory=10e9 表示分配10GB内存。

• tmpdir

  • 类型: 字符串
  • 作用: 指定临时文件存储目录（需位于高速存储设备以提升性能）。
  • 示例: tmpdir = "G:/temp"

• chunksize

  • 类型: 数值
  • 作用: 定义分块处理数据的大小，优化大文件处理速度。

2. stack()
功能: 将多个栅格文件合并为RasterStack对象（多层栅格）。
参数:

stack(..., bands=NULL, quick=FALSE)

• ...

  • 类型: 字符向量或栅格文件路径
  • 作用: 输入待合并的栅格文件列表。
  • 示例:

    stack("G:/climate/temp.tif", "G:/climate/precip.tif")
    

• bands

  • 类型: 整数
  • 作用: 选择多波段文件中的特定波段（默认为所有波段）。

3. writeRaster()
功能: 将栅格数据写入文件（如GeoTIFF）。
参数:

writeRaster(x, filename, format="GTiff", overwrite=FALSE, ...)

• x

  • 类型: RasterLayer 或 RasterStack
  • 作用: 待保存的栅格数据。

• filename

  • 类型: 字符串
  • 作用: 输出文件路径及名称。
  • 示例: filename = "G:/output/prediction.tif"

• format

  • 类型: 字符串
  • 作用: 指定文件格式（如 "GTiff"、"HFA"（Erdas Imagine））。

• overwrite

  • 类型: 逻辑值
  • 作用: 是否覆盖同名文件（默认FALSE）。

4. predict()
功能: 与dismo包中的predict()类似，但更通用，支持多种模型（如GLM、随机森林）。
参数:

predict(object, newdata, filename="", ...)

• newdata
  • 类型: RasterStack
  • 作用: 输入的环境变量栅格数据。

1.3 常见问题与解决

1. Java内存不足

   • 错误提示: Java.lang.OutOfMemoryError
   • 解决:

     options(java.parameters = "-Xmx20g")  # 分配更大内存
     

2. 临时文件占用过大

   • 现象: 临时目录爆满导致程序崩溃。
   • 解决:

     rasterOptions(tmpdir = "D:/fast_disk/temp")  # 指定大容量高速存储
     

3. 预测结果全为NA

   • 原因: 环境变量栅格与存在点坐标的空间参考不一致。
   • 解决:

     crs(env_stack) <- crs(presence_points)  # 检查并统一坐标系
     

二、代码示例

只是示例，你自己问AI即可。

options(java.parameters = "-Xmx10g")  # 分配10GB内存给Java

library(dismo)
library(raster)

rasterOptions(maxmemory = 10e9)  # 10GB内存
rasterOptions(tmpdir = "xxxxx")  # 设置高速存 储上的临时目录

# 读取环境变量数据
env_files <- list.files(
  path = "xxx",
  pattern = ".tif$",
  full.names = TRUE
)
env_stack <- stack(env_files)

# 读取物种存在点
species_data <- read.csv("xxx点位_albers.csv")
presence_points <- species_data[, c("x", "y")]

# 创建列表存储每次的预测结果
prediction_list <- list()

# 一次为例
for (i in 1:1) {
  # 设置不同的随机种子确保数据分割不同
  set.seed(i)

  # 分割训练集和测试集 (80/20)
  train_indices <- sample(1:nrow(presence_points), 0.8 * nrow(presence_points))
  train_points <- presence_points[train_indices, ]
  test_points <- presence_points[-train_indices, ]

  # 训练MaxEnt模型
  maxent_model <- maxent(
    x = env_stack,
    p = train_points,
    removeDuplicates = TRUE,
    path = paste0("xxx/run_", i),  # 为每次运行创建独立文件夹
    args = c(
      "betamultiplier=1",
      "linear=TRUE",
      "quadratic=TRUE",
      "product=TRUE",
      "hinge=TRUE",
      "threshold=FALSE",
      "threads=12"
    )
  )

  # 预测分布概率
  prediction_map <- predict(
    maxent_model,
    env_stack,
    progress = 'text',
    filename = paste0("xxx/predict_", i, ".tif")  # 保存预测结果
  )

  # 将预测结果存入列表
  prediction_list[[i]] <- prediction_map
}

# # 合并所有预测结果并计算平均值
# prediction_stack <- stack(prediction_list)
# average_prediction <- mean(prediction_stack)
# 
# # 保存平均结果
# writeRaster(
#   average_prediction,
#   filename = "G:/人兽冲突/R_code/out/average_prediction.tif",
#   overwrite = TRUE
# )
# 
# # 可视化平均预测
# plot(average_prediction, main = "Average Species Distribution Probability")
# points(presence_points, pch = 19, cex = 0.5, col = "red")