从逐日气象数据生成生物气候变量

25 年 4 月 15 日 星期二
2789 字
14 分钟

生物气候变量可以用于最大熵模型（Maxent，Maximum Entropy Model）中的生物信息学，用于物种分布预测，通过分析物种已知分布数据和环境变量，推断物种在未知区域的分布概率。一般来说，所使用的环境变量采用的世界气候数据库中的数据(http://www.worldclim.org)，但是WorldClim里的数据并不是最新的，历史数据只包含1970-2000年，2000年以后的数据都是预测数据，用来做当年的分析多多少少不够准确

怎样才能自己用数据创建**生物气候变量？**WorldClim的原话是这样的：

To create these values yourself, you can use the ‘biovars’ function in the R package dismo

也就是说，我们需要使用R语言进行计算，函数具体信息如下：

r
# Function to create 'bioclimatic variables' from monthly climate data.

# Usage:
#prec    vector, matrix, or RasterStack/Brick of precipitation data
#tmin    vector, matrix, or RasterStack/Brick of minimum temperature data
#tmax    vector, matrix, or RasterStack/Brick of maximum temperature data
#...     Additional arguments
# Input data is normaly monthly. I.e. there should be 12 values (layers) for each variable, but the function should also work for e.g. weekly data (with some changes in the meaning of the output variables. E.g. #8 would then not be for a quater (3 months), but for a 3 week period)

biovars(prec, tmin, tmax, ...)

从文档得知，要使用这个函数，需要准备12个月每月的平均降水、平均最低温和平均最高温数据。数据哪里来？可以参看这篇文章《中国气象站点数据获取的几种方式 》。我选择了难度最低的方法，在NCEI获取数据，下面就以2024年的逐日数据为例，逐步处理最后得到生物气候变量

站点气候数据获取与粗处理

我下载的是2024年逐日数据的归档数据，范围为全球，所以首先要将中国站点的数据找到，中国目前被收录的站点编号为50136099999-59997099999，每一个文件内的格式如下：

csv
"STATION","DATE","LATITUDE","LONGITUDE","ELEVATION","NAME","TEMP","TEMP_ATTRIBUTES","DEWP","DEWP_ATTRIBUTES","SLP","SLP_ATTRIBUTES","STP","STP_ATTRIBUTES","VISIB","VISIB_ATTRIBUTES","WDSP","WDSP_ATTRIBUTES","MXSPD","GUST","MAX","MAX_ATTRIBUTES","MIN","MIN_ATTRIBUTES","PRCP","PRCP_ATTRIBUTES","SNDP","FRSHTT"
"50136099999","2024-01-01","52.9666666","122.5333333","438.0","MOHE, CH"," -17.1"," 8"," -24.6"," 8","1021.6"," 8","962.5"," 8","  8.0"," 8","  1.5"," 8","  2.7","999.9","  -0.4"," "," -32.8"," "," 0.03","G","999.9","001000"
"50136099999","2024-01-02","52.9666666","122.5333333","438.0","MOHE, CH"," -12.3"," 8"," -18.8"," 8","1018.3"," 8","958.8"," 8","  4.3"," 8","  0.8"," 8","  1.6","999.9","  -0.4"," "," -31.2"," "," 0.03","G","  7.5","001000"

将所有国内站点的数据放到新的文件夹，我这里命名为“2024中国气象站逐日数据”，然后将所有文件进行合并，以便后续分析，这里我使用Python进行处理：

python
import pandas as pd
import glob

# 读取数据
csv_files = glob.glob("2024中国气象站逐日数据/*.csv")

# 读取所有CSV文件并合并
combined_data = pd.DataFrame()
for file in csv_files:
    data = pd.read_csv(file, encoding='utf-8')
    combined_data = pd.concat([combined_data, data], ignore_index=True)
# 保存文件
combined_data.to_csv("2024中国气象站逐日数据-合并后.csv", index=False, encoding='utf-8')

随后对合并后的数据进行清洗和计算月均值，同样使用Python：

python
import pandas as pd
from enum import Enum

class MssingValue(Enum):
    """
    缺失值
    """
    MAX = 9999.9
    MIN = 9999.9
    PRCP = 99.99
# 计算每个站点的月均值,并将数据转换为中国标准
def calculate_monthly_means(df):
    """
    计算每个站点的月均值
    :param df: 分组后的数据
    :return: 月均值数据
    """
    # 根据缺失值进行数据清洗
    df = df[df['MAX'] != MssingValue.MAX.value]
    df = df[df['MIN'] != MssingValue.MIN.value]
    df = df[df['PRCP'] != MssingValue.PRCP.value]

    # 将日期列转换为日期类型
    df['DATE'] = pd.to_datetime(df['DATE'], format='%Y-%m-%d')

    # 提取月份
    df['MONTH'] = df['DATE'].dt.month

    monthly_means = df.groupby(['STATION', 'MONTH']).agg({
        'MAX': 'mean',
        'MIN': 'mean',
        'PRCP': 'sum',  # 降水量求和
    }).reset_index()
    # 将月份转换为字符串格式
    monthly_means['MONTH'] = monthly_means['MONTH'].astype(str).str.zfill(2)
    # 每个站点的经纬度
    station_info = df[['STATION', 'LATITUDE', 'LONGITUDE', 'ELEVATION', 'NAME']].drop_duplicates(subset='STATION')
    monthly_means = monthly_means.merge(station_info, on='STATION', how='left')

    monthly_means['MAX'] = (monthly_means['MAX'] - 32) / 1.8  # 华氏度转摄氏度
    monthly_means['MIN'] = (monthly_means['MIN'] - 32) / 1.8  # 华氏度转摄氏度
    monthly_means['PRCP'] = monthly_means['PRCP'] * 25.4  # 英寸转毫米

    return monthly_means

# 处理数据
def process_data(file_path):
    """
    处理数据
    :param file_path: 文件路径
    :return: 处理后的数据
    """
    # 读取数据
    df = pd.read_csv(file_path)

    # 计算每个站点的月均值
    monthly_means = calculate_monthly_means(df)

    return monthly_means
# 主函数
if __name__ == "__main__":
    # 文件路径
    file_path = "2024中国气象站逐日数据-合并后.csv"
    # 处理数据
    processed_data = process_data(file_path)
    processed_data.to_csv("2024中国气象站逐日数据-月均值.csv", index=False)

此时得到的 2024中国气象站逐日数据-合并后.csv结构应该是这样的：

csv
STATION,MONTH,MAX,MIN,PRCP,LATITUDE,LONGITUDE,ELEVATION,NAME
50136099999,01,-15.600427350427351,-36.113247863247864,5.08,52.9666666,122.5333333,438.0,"MOHE, CH"
50136099999,02,-10.016460905349795,-34.619341563786,6.096,52.9666666,122.5333333,438.0,"MOHE, CH"

再将这个文件按月拆分开，同样使用Python：

python
import pandas as pd

# 读取数据
df = pd.read_csv('2024中国气象站逐日数据-月均值.csv')

# 按月分割数据
def split_monthly_data(df):
    # 将 'MONTH' 列转换为整数类型
    df['MONTH'] = df['MONTH'].astype(int)
    # 遍历每个月份
    for month in range(1, 13):
        # 过滤出当前月份的数据
        monthly_data = df[df['MONTH'] == month]
        # 将数据保存到 CSV 文件中，文件名为 'month_XX.csv'
        monthly_data.to_csv(f'2024中国气象站逐月数据/month_{month:02}.csv', index=False)
        print(f"Month {month:02} data saved to month_{month:02}.csv")

if __name__ == "__main__":
    split_monthly_data(df)

这个时候我们的工作目录结构应该是这样的：

sh
.
├── 按月分割月均值数据.py
├── 合并气候数据.py
├── 处理气候数据.py
├── 2024中国气象站逐日数据
│   ├── 50136099999.csv
│   ├── 50353099999.csv
│   ├── 50434099999.csv
│   ├── 50468099999.csv
│   ├── 50527099999.csv
│   └── ............csv
├── 2024中国气象站逐月数据
│   ├── month_01.csv
│   ├── month_02.csv
│   ├── month_03.csv
│   ├── month_04.csv
│   ├── month_05.csv
│   ├── month_06.csv
│   ├── month_07.csv
│   ├── month_08.csv
│   ├── month_09.csv
│   ├── month_10.csv
│   ├── month_11.csv
│   └── month_12.csv
├── 2024中国气象站逐日数据-合并后.csv
└── 2024中国气象站逐日数据-月均值.csv

月度气候数据插值

现在得到的数据都是点类型的数据，并没有对国内进行覆盖，没有数据的坐标点就需要进行插值计算获取。这里我们使用R语言进行普通克里金法插值：

为什么要插值？插值是由有限数量的采样点数据估计栅格中的单元的值。它可以用来估计任何地理点数据的未知值：高程、降雨、化学污染程度、噪声等级等等。在研究区域内，测量某种现象每个点的高度、等级或集聚程度一般是非常困难，同时也是很昂贵的。相反，用户可以选择一些离散的样本点进行测量，通过插值得出采样点的值。采样点可以是随机的、分层的或者规则的格网点，包含高度、污染程度或者等级等信息。具体概念和实现可以看这个论坛:Python气象数据空间插值与绘图

r
# 加载必要的包
library(raster)
library(sp)
library(rgdal)
library(gstat)
library(maptools)

output_dir <- "逐月数据插值结果"
input_dir <- "2024中国气象站逐月数据"

# 创建输出目录
dir.create(output_dir)

# 读取中国边界数据
china_map <- readOGR("China.shp") #需要自备中国范围的矢量文件

# 创建覆盖中国的规则网格
# 转换坐标系统为更适合中国区域的投影
china_proj <- spTransform(
    china_map,
    CRS("+proj=laea +lat_0=36 +lon_0=104 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs")
)

# 创建规则网格(使用相同的分辨率)
grid_res <- 0.04165 # 分辨率，单位为度
china_extent <- extent(china_map)
x_seq <- seq(china_extent@xmin, china_extent@xmax, by = grid_res)
y_seq <- seq(china_extent@ymin, china_extent@ymax, by = grid_res)
grid <- expand.grid(x = x_seq, y = y_seq)
coordinates(grid) <- ~ x + y
gridded(grid) <- TRUE
proj4string(grid) <- proj4string(china_map)

# 只保留中国边界内的网格点
grid_in_china <- grid[china_map, ]

# 定义插值处理函数
process_monthly_data <- function(month_num) {
    # 月份编号格式化(补零)
    month_str <- sprintf("%02d", month_num)

    cat(paste0("处理第 ", month_str, " 月数据...\n"))

    # 读取气象数据
    file_path <- paste0(input_dir, "/month_", month_str, ".csv")
    climate_data <- read.csv(file_path)

    # 转换为空间点数据框
    coordinates(climate_data) <- ~ LONGITUDE + LATITUDE
    proj4string(climate_data) <- CRS("+proj=longlat +datum=WGS84")

    # 创建输出文件名前缀
    output_prefix <- paste0(output_dir, "/month_", month_str, "_")

    # 用普通克里金法插值平均温度(TEMP)
    cat("  插值平均温度(TEMP)...\n")
    v_temp <- variogram(TEMP ~ 1, climate_data)
    v_temp_fit <- fit.variogram(v_temp, model = vgm("Sph"))
    temp_kriging <- krige(TEMP ~ 1, climate_data, grid_in_china, v_temp_fit)
    temp_raster <- raster(temp_kriging)
    temp_raster_masked <- mask(temp_raster, china_map)
    writeRaster(temp_raster_masked, paste0(output_prefix, "temp.tif"),
        format = "GTiff", overwrite = TRUE
    )

    # 对最高温度(MAX)插值
    cat("  插值最高温度(MAX)...\n")
    v_max <- variogram(MAX ~ 1, climate_data)
    v_max_fit <- fit.variogram(v_max, model = vgm("Sph"))
    max_kriging <- krige(MAX ~ 1, climate_data, grid_in_china, v_max_fit)
    max_raster <- raster(max_kriging)
    max_raster_masked <- mask(max_raster, china_map)
    writeRaster(max_raster_masked, paste0(output_prefix, "max.tif"),
        format = "GTiff", overwrite = TRUE
    )

    # 对最低温度(MIN)插值
    cat("  插值最低温度(MIN)...\n")
    v_min <- variogram(MIN ~ 1, climate_data)
    v_min_fit <- fit.variogram(v_min, model = vgm("Sph"))
    min_kriging <- krige(MIN ~ 1, climate_data, grid_in_china, v_min_fit)
    min_raster <- raster(min_kriging)
    min_raster_masked <- mask(min_raster, china_map)
    writeRaster(min_raster_masked, paste0(output_prefix, "min.tif"),
        format = "GTiff", overwrite = TRUE
    )

    # 对降水量(PRCP)插值
    cat("  插值降水量(PRCP)...\n")
    v_prcp <- variogram(PRCP ~ 1, climate_data)
    v_prcp_fit <- fit.variogram(v_prcp, model = vgm("Sph"))
    prcp_kriging <- krige(PRCP ~ 1, climate_data, grid_in_china, v_prcp_fit)
    prcp_raster <- raster(prcp_kriging)
    # 将降水量负值替换为0
    prcp_raster[prcp_raster < 0] <- 0
    prcp_raster_masked <- mask(prcp_raster, china_map)
    writeRaster(prcp_raster_masked, paste0(output_prefix, "prcp.tif"),
        format = "GTiff", overwrite = TRUE
    )

    cat(paste0("第 ", month_str, " 月数据处理完成\n"))
}

# 批量处理所有月份
start_time <- Sys.time()

for (month in 1:12) {
    process_monthly_data(month)
}

end_time <- Sys.time()
processing_time <- difftime(end_time, start_time, units = "mins")

cat(paste0("\n所有月份数据处理完成！\n"))
cat(paste0("总耗时: ", round(processing_time, 2), " 分钟\n"))
cat(paste0("结果文件已保存到", output_dir, "目录\n"))

运行后得到全国的逐月三种气候数据：

sh
.
├── ...
└── 逐月数据插值结果
    ├── month_01_max.tif
    ├── month_01_min.tif
    ├── month_01_prcp.tif
    └── ..............tif

计算生物气候变量

现在使用R语言dismo包中的biovars方法计算生物气候变量：

r
# 加载必要的包
library(raster)
library(dismo)
library(rgdal)

# 设置输入和输出目录
input_dir <- "逐月数据插值结果"
output_dir <- "逐月数据插值结果计算bioclims"
dir.create(output_dir, showWarnings = FALSE, recursive = TRUE)

cat("开始读取每月气候栅格数据...\n")

# 读取每月温度和降水栅格数据
tmax_stack <- stack()
tmin_stack <- stack()
prcp_stack <- stack()

# 按顺序读取12个月的数据
for (month in 1:12) {
  month_str <- sprintf("%02d", month)
  # 构建文件路径
  max_file <- file.path(input_dir, paste0("month_", month_str, "_max.tif"))
  min_file <- file.path(input_dir, paste0("month_", month_str, "_min.tif"))
  prcp_file <- file.path(input_dir, paste0("month_", month_str, "_prcp.tif"))
  # 检查文件是否存在并读取栅格
  if (file.exists(max_file)) {
    tmax_stack <- addLayer(tmax_stack, raster(max_file))
  } else {
    stop(paste("文件不存在:", max_file))
  }
  if (file.exists(min_file)) {
    tmin_stack <- addLayer(tmin_stack, raster(min_file))
  } else {
    stop(paste("文件不存在:", min_file))
  }
  if (file.exists(prcp_file)) {
    prcp_stack <- addLayer(prcp_stack, raster(prcp_file))
  } else {
    stop(paste("文件不存在:", prcp_file))
  }
}

# 设置栅格名称以表示月份
names(tmax_stack) <- paste0("tmax_", 1:12)
names(tmin_stack) <- paste0("tmin_", 1:12)
names(prcp_stack) <- paste0("prcp_", 1:12)

cat("正在计算19个标准生物气候变量...\n")

# 使用dismo包中的biovars函数计算19个生物气候变量
bioclim_vars <- biovars(prcp_stack, tmin_stack, tmax_stack)

# 设置生物气候变量的名称
bioclim_names <- c(
  "年平均温度",
  "平均昼夜温差",
  "等温性",
  "温度季节性变化",
  "最热月最高温度",
  "最冷月最低温度",
  "年温度变化范围",
  "最湿季平均温度",
  "最干季平均温度",
  "最热季平均温度",
  "最冷季平均温度",
  "年降水量",
  "最湿月降水量",
  "最干月降水量",
  "降水季节性变化",
  "最湿季降水量",
  "最干季降水量",
  "最热季降水量",
  "最冷季降水量"
)

# 设置栅格名称
names(bioclim_vars) <- paste0("BIO", 1:19)

# 保存各个生物气候变量为单独的GeoTIFF文件
cat("正在保存19个生物气候变量栅格...\n")

for (i in 1:19) {
  var_name <- paste0("BIO", i)
  file_name <- file.path(output_dir, paste0(var_name, ".tif"))
  writeRaster(bioclim_vars[[i]], file_name, format = "GTiff", overwrite = TRUE)
  cat(paste(" -", var_name, "(", bioclim_names[i], ")", "已保存\n"))
}

cat("\n19个标准生物气候变量计算完成！\n")
cat(paste0("所有结果已保存到: ", output_dir, "\n"))

到此，我们就得到了2024年的19个生物气候变量的数据：

sh
.
├── ...
└── 逐月数据插值结果计算bioclims
    ├── BIO1.tif
    ├── BIO10.tif
    ├── BIO11.tif
    ├── BIO12.tif
    └── ......tif

文章标题：从逐日气象数据生成生物气候变量

文章作者：Riceneeder

文章链接：https://gankun.cn.ma/posts/2025-04-15[复制]

最后修改时间：

商业转载请联系站长获得授权，非商业转载请注明本文出处及文章链接，您可以自由地在任何媒体以任何形式复制和分发作品，也可以修改和创作，但是分发衍生作品时必须采用相同的许可协议。
本文采用CC BY-NC-SA 4.0进行许可。