從逐日氣象數據生成生物氣候變量

生物氣候變量可以用於最大熵模型（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 語言進行計算，函數具體信息如下：

# 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，每一個文件內的格式如下：

"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 進行處理：

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：

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結構應該是這樣的：

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：

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)

這個時候我們的工作目錄結構應該是這樣的：

.
├── 按月分割月均值數據.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 氣象數據空間插值與繪圖

# 加載必要的包
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"))

運行後得到全國的逐月三種氣候數據：

.
├── ...
└── 逐月數據插值結果
    ├── month_01_max.tif
    ├── month_01_min.tif
    ├── month_01_prcp.tif
    └── ..............tif

計算生物氣候變量#

現在使用 R 語言dismo包中的biovars方法計算生物氣候變量：

# 加載必要的包
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 個生物氣候變量的數據：

.
├── ...
└── 逐月數據插值結果計算bioclims
    ├── BIO1.tif
    ├── BIO10.tif
    ├── BIO11.tif
    ├── BIO12.tif
    └── ......tif