逐日気象データから生物気候変数を生成する

生物気候変数は最大エントロピーモデル（Maxent、Maximum Entropy Model）における生物情報学に使用され、種分布予測に役立ちます。これは、種の既知の分布データと環境変数を分析することによって、未知の地域における種の分布確率を推測します。一般的に、使用される環境変数は世界気候データベースのデータを使用します (http://www.worldclim.org) が、WorldClim のデータは最新のものではなく、歴史データは 1970-2000 年のみを含み、2000 年以降のデータは予測データであり、当年の分析には多かれ少なかれ正確さに欠けます。

自分でデータを使用して生物気候変数を作成するにはどうすればよいでしょうか？WorldClim の原文は次のようになります：

これらの値を自分で作成するには、_R_パッケージ dismo の‘biovars’関数を使用できます。

つまり、R 言語を使用して計算する必要があります。関数の具体的な情報は以下の通りです：

# 月次気候データから「生物気候変数」を作成する関数。

# 使用法：
#prec    降水データのベクトル、行列、またはRasterStack/Brick
#tmin    最低気温データのベクトル、行列、またはRasterStack/Brick
#tmax    最高気温データのベクトル、行列、またはRasterStack/Brick
#...     追加の引数
# 入力データは通常月次です。つまり、各変数に対して12の値（レイヤー）が必要ですが、関数は例えば週次データにも対応するはずです（出力変数の意味にいくつかの変更が必要です。例えば、#8は四半期（3ヶ月）ではなく、3週間の期間になります）。

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"))

実行後、全国の逐月 3 種類の気候データが得られます：

.
├── ...
└── 逐月データ補間結果
    ├── 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