使用R语言和APSIM模型模拟不同经纬度土壤条件下的作物生长

library(data.table)

# 链接apsim.exe
apsim_exe <- 'D:/APSIM/APSIM710-r4158/Model/APSIM.exe'

# 读取土壤信息文件
df_soil <- fread('D:/APSIM_work/scenario_585/soil_info.csv')

# 循环遍历每个土壤信息
for (i in 1:nrow(df_soil)) {
  # 读取apsim文件
apsim <- readLines('D:/APSIM_work/scenario_585/SW/chuan_HT33.apsim')
  
  # 替换经纬度和土壤信息
apsim[4] <- gsub('57237', df_soil$station_id[i], apsim[4])
apsim[12] <- gsub('0.150', df_soil$clay_pct[i], apsim[12])
apsim[13] <- gsub('0.080', df_soil$sand_pct[i], apsim[13])
apsim[14] <- gsub('0.770', df_soil$om_pct[i], apsim[14])
apsim[15] <- gsub('120', df_soil$bd[i], apsim[15])
apsim[16] <- gsub('6.5', df_soil$ph[i], apsim[16])
  
  # 设置met文件路径
met_site <- paste0('D:/APSIM_work/scenario_585/met_585/', df_soil$station_id[i], '_SSP585_BCCC_r1i1p1f1.MET')
  
  ## 替换气象数据####
met <- readLines(met_site, n = 33)
b <- which(grepl('! AMP', met))
  
  # 替换CO2浓度
apsim[583] <- gsub('!', '', met[b + 6]) # C02
  
  # 读取met文件数据
met_data <- fread(met_site, skip = 35)
  # 固定1961-2019
met_data1 <- subset(met_data, met_data$V1 %in% c(1961:2019))
  # 修改2020-2100
met_data2 <- subset(met_data, met_data$V1 %in% c(2020:2100))
met_data2 <- subset(met_data2, met_data2$V2 != 366)
  
  # 替换2020-2099年的辐射数据
a <- met_data$V3[which(met_data$V1 == 2020 & met_data$V2 == 1):which(met_data$V1 == 2029 & met_data$V2 == 365)]
met_data2$V3[which(met_data2$V1 == 2020 & met_data$V2 == 1):which(met_data2$V1 == 2099 & met_data$V2 == 365)] <- rep(a, 8)
  
  # 补全闰年
for (year in seq(2020, 2096, 4)) {
    n <- which(met_data2$V1 == year & met_data2$V2 == 365)
    df_up <- met_data2[1:n, ]
    df_bottom <- met_data2[(n + 1):nrow(met_data2), ]
    insert_row <- df_up[nrow(df_up), ]
    insert_row$V2 <- 366
    met_data2 <- rbind(df_up, insert_row, df_bottom)
  }
  
  # 将修改后的met数据写入文件
writeLines(met, met_site)
write.table(met_data2, met_site, append = T, row.names = F, col.names = F)
  
  # 将修改后的apsim文件写入新文件
file_apsim <- paste0('D:/APSIM_work/scenario_585/SW/chuan_HT33_run_', df_soil$station_id[i], '.apsim')
writeLines(apsim, file_apsim)

  # 运行模型
system(paste(apsim_exe, file_apsim), wait = TRUE, ignore.stdout = TRUE, show.output.on.console = FALSE)
  
  # 将输出文件复制到新目录
file.copy('D:/APSIM_work/scenario_585/SW/chuan_HT33.out', paste0('D:/APSIM_work/result585_radn/SW/HT33/', df_soil$station_id[i], '.out'))
  
  print(paste0('Finished running APSIM for station ', df_soil$station_id[i]))
}

# 将所有输出文件合并为一个CSV文件
out_files <- list.files('D:/APSIM_work/result585_radn/SW/HT33', pattern = '\.out$', full.names = TRUE)
out_data <- do.call(rbind, lapply(out_files, fread, skip = 9))
write.csv(out_data, 'D:/APSIM_work/result585_radn/SW/HT33_all_out.csv', row.names = FALSE)