利用GSOD数据集计算参考作物蒸散量ET0

本文参考了根据gsod数据集计算潜在蒸发量（Matlab）这篇博客。

GSOD数据集地址

直接选择红色框中的下载csv

现在我想计算出库尔勒地区2008年和2009年的ET0值，找到2008和2009年的数据

以2008年为例，点进去全是一串串的数字，我无法识别哪个是库尔勒气象站的文件。但我了解到数字的前六位是气象站的USAF编号，后面的99999是wban编号，所以现在的问题就是要找到库尔勒气象站的USAF编号。

我发现一个stationaRy的R语言包可以解决这个问题。github地址：stationaRy

下面用R代码找到库尔勒气象站的USAF编号

library(stationaRy)

stations_KORLA <- 
  get_station_metadata() %>%
  dplyr::filter(name == "KORLA") # 用名字过滤，注意名字需要全部大写

print(stations_KORLA)

可以看到库尔勒气象站的USAF编号和一些基本信息

接下来找到对应的CSV文件就可以下载了

接下来，将下载好的CSV文件导入mysql数据库，方便后面查询和相关操作。我创建了GSOD数据库，然后导入了KORLA_2008和KORLA_2009数据表。

数据表每列所代表的含义可以去GSOD数据集地址中查看

readme截图，注意每列参数的单位，与我们常用的单位有些不同，如温度都是华氏温度，风速是节konts，1konts = 0.5144m/s，后面计算时需要单位转换。

下面我定义一个connect_db.py文件来连接数据库以及获取数据表中的某一列，这里为了简便我没有使用连接池，没有考虑性能要求。

import pymysql


def fetch_column_from_table(table_name, column_name):
    """
    查询数据表中的某列
    :param table_name: 表的名称
    :param column_name: 列名
    :return: 返回一个数组
    """
    # 连接到 MySQL 数据库
    connection = pymysql.connect(
        host="localhost",
        user="root",
        password="123456",
        database="GSOD"
    )
    # 创建游标对象
    cursor = connection.cursor()

    try:
        # 执行查询
        query = f"SELECT {column_name} FROM {table_name}"
        cursor.execute(query)

        # 获取结果
        column_data = cursor.fetchall()

        # 转换结果为数组
        result = [row[0] for row in column_data]
        return result

    finally:
        # 关闭游标和连接
        cursor.close()
        connection.close()

然后创建一个get_data.py文件用来接收所需原始气象数据方便后续计算ET0。

from connect_db import fetch_column_from_table

# 获取平均温度、最高温度、最低温度、风速、露点温度数据
temp_mean_2008 = fetch_column_from_table('KORLA_2008', 'TEMP')
temp_max_2008 = fetch_column_from_table('KORLA_2008', 'MAX')
temp_min_2008 = fetch_column_from_table('KORLA_2008', 'MIN')
wind_speed_2008 = fetch_column_from_table('KORLA_2008', 'WDSP')
temp_dew_2008 = fetch_column_from_table('KORLA_2008', 'DEWP')

temp_mean_2009 = fetch_column_from_table('KORLA_2009', 'TEMP')
temp_max_2009 = fetch_column_from_table('KORLA_2009', 'MAX')
temp_min_2009 = fetch_column_from_table('KORLA_2009', 'MIN')
wind_speed_2009 = fetch_column_from_table('KORLA_2009', 'WDSP')
temp_dew_2009 = fetch_column_from_table('KORLA_2009', 'DEWP')

接下来创建一个GSOD_toCalculate_ET0.py文件计算ET0。计算ET0采用的是彭曼公式，具体可以参考其它一些资料，这里不赘述。

import math

import numpy as np

from get_data import (temp_max_2008, temp_min_2008, temp_dew_2008,
                      wind_speed_2008, temp_mean_2008, temp_mean_2009,
                      temp_dew_2009, temp_max_2009, temp_min_2009, wind_speed_2009)


def calculate_et0(temp_mean, temp_max, temp_min, wind_speed,
                  temp_dew, day_of_year, lat=41.73, elevation=903):
    """

    :param temp_mean:  平均温度，单位是华氏温度
    :param temp_max:   最高温度，单位是华氏温度
    :param temp_min:   最低温度，单位是华氏温度
    :param wind_speed: 风速，单位是节kont
    :param lat:  站点纬度，单位为度
    :param temp_dew:  露点温度，单位是华氏温度
    :param elevation:  高程（m）
    :param day_of_year: 一年中的第几天，日序
    :return: 返回潜在蒸散量值（mm/day）
    """

    # 华氏温度转为摄氏度
    temp_mean = (temp_mean - 32) / 1.8
    temp_max = (temp_max - 32) / 1.8
    temp_min = (temp_min - 32) / 1.8
    temp_dew = (temp_dew - 32) / 1.8

    # 将摄氏度转换为开尔文温度
    T_mean = temp_mean + 273.16
    T_max = temp_max + 273.16
    T_min = temp_min + 273.16
    T_dew = temp_dew + 273.16

    # 风速单位转换为m/s
    wind_speed = 0.5144 * wind_speed

    # 纬度单位转换为弧度
    lat = lat * math.pi / 180

    # 计算平均饱和蒸气压
    es_max = 0.6108 * math.exp((17.27 * temp_max) / (temp_max + 237.3))
    es_min = 0.6108 * math.exp((17.27 * temp_min) / (temp_min + 237.3))
    es = (es_max + es_min) / 2

    # 计算实际蒸气压
    ea = 0.6108 * math.exp(17.27 * temp_dew / (temp_dew + 273.3))

    # 计算饱和蒸气压斜率
    delta = (4098 * es_max) / ((temp_max + 273.3) ** 2)

    # 计算gamma
    P = 101.3 * ((293 - 0.0065 * elevation) / 293) ** 5.26
    gamma = 0.665 * 10 ** (-3) * P

    # 计算晴空太阳辐射Rso
    Gsc = 0.082  # 太阳常数
    dr = 1 + 0.033 * math.cos(2 * math.pi * day_of_year / 365)  # 日地间相对距离的倒数
    Sdelta = 0.408 * math.sin(2 * math.pi * day_of_year / 365 - 1.39)  # 太阳磁偏角，单位为弧度
    ws = math.acos(-math.tan(lat) * math.tan(Sdelta))  # 日落时角，单位为弧度
    Ra = 24 * 60 / math.pi * Gsc * dr * (ws * math.sin(lat) * math.sin(Sdelta) +
                                         math.cos(lat) * math.cos(Sdelta) * math.sin(ws))
    Krs = 0.162
    Rs = Krs * Ra * (temp_max - temp_min) ** 0.5  # 太阳辐射Rs
    Rso = (0.75 + 2 * 10 ** (-5) * elevation) * Ra

    # 净太阳辐射Rns
    alpha = 0.23
    Rns = (1 - alpha) * Rs

    # 净长波辐射Rnl
    sigma = 4.903 * 10 ** (-9)
    air_correction = 0.34 - 0.14 * ea ** 0.5  # 空气湿度的订正项
    cloud_correction = 1.35 * Rs / Rso - 0.35
    Rnl = sigma * (T_max ** 4 + T_min ** 4) / 2 * air_correction * cloud_correction

    # 净辐射Rn
    Rn = Rns - Rnl

    # 土壤热通量
    G = 0

    # 2m高处风速，由Zm高处风速转来
    u2 = wind_speed * (4.87 / math.log(67.8 * 10 - 5.42))

    # 计算ET0
    ET0 = (0.408 * delta * (Rn - G) + gamma * 900 / (temp_mean + 273) * u2 * (es - ea)) / (
            delta + gamma * (1 + 0.34 * u2))
    return ET0


ET0_2008 = []
ET0_2009 = []
for i in range(len(temp_mean_2008)):
    ET0_2008.append(calculate_et0(temp_mean_2008[i], temp_max_2008[i], temp_min_2008[i],
                                  wind_speed_2008[i], temp_dew_2008[i], day_of_year=i + 1))
for j in range(len(temp_mean_2009)):
    ET0_2009.append(calculate_et0(temp_mean_2009[j], temp_max_2009[j], temp_min_2009[j],
                                  wind_speed_2009[j], temp_dew_2009[j], day_of_year=j + 1))

# with open('ET0_2008.txt', 'a') as file:
#     for m in range(len(ET0_2008)):
#         file.write(f'{ET0_2008[m]}\n')
#
# with open('ET0_2009.txt', 'a') as file:
#     for n in range(len(ET0_2009)):
#         file.write(f'{ET0_2009[n]}\n')

现在我们就获取到了库尔勒地区2008年和2009年的ET0值。

创建一个plot_ET0.py文件绘制2008年库尔勒地区全年ET0

import matplotlib.pyplot as plt
import pandas as pd
from GSOD_toCalculate_ET0 import ET0_2008

# 利用pandas库生成日期
dates = pd.date_range(start='2008-01-01', end='2008-12-31', freq='D')

# 创建一个Pandas DataFrame，将日期和数据组合起来
df = pd.DataFrame({'Date': dates, 'Data': ET0_2008})

# 将日期列设置为DataFrame的索引
df.set_index('Date', inplace=True)

# 绘制折线图
plt.figure(figsize=(10, 6))
plt.plot(df.index, df['Data'], marker='o', linestyle='-')
plt.title('Daily ET0 Data in 2008')
plt.xlabel('Date')
plt.ylabel('ET0(mm/day)')
plt.show()

结果展示：