[SWAT·站点气象数据制备] 批量生成站点气象数据文件代码

特别感谢CSDN博主“摆烂老大”的SWAT教程帖子对我的指导和帮助！

首先先简短讲下SWAT站点气象数据的制备（不包括天气发生器的构建）：

注：SWAT气象部分工作=站点气象数据制备+天气发生器构建（构建所需参数可用SwatWeather程序计算）

SWAT模型的站点气象数据库的制备，最终要制作5项气象数据，分别是日降水量（PCP）、最高/最低气温（TMP）、太阳辐射（SLR）、相对湿度（HMD）、平均风速（WND），按要素类型建好五个文件夹，每项气象要素（每个文件夹）中包含的内容要有①站点文件和②站点气象数据文件，都用txt文本文件来储存。站点文件是对气象站点基本信息的汇总，站点气象数据文件是对具体的气象数据的存储。一个站点文件往往对应多个站点气象数据文件。

以降水某一站点为例，形成的该站点的气象数据文件内容如下：

详细的SWAT教程在这位博主（摆烂老大）的帖子中，包括本部分的气象数据制备内容：

https://blog.csdn.net/weixin_67437706/category_12886961.html?fromshare=blogcolumn&sharetype=blogcolumn&sharerId=12886961&sharerefer=PC&sharesource=MENGCCSF&sharefrom=from_link

这位博主介绍的很详细，非常推荐，可以跟着这个博主的教程做下去

制备站点气象数据先从原始气象数据excel中整理出来按站点分类的气象数据的几个excel，每个excel是一个站点里对应原始气象数据的各列指标，然后用代码读取这些excel批量一次性生成这些气象数据文件（.txt），最后复制粘贴到五个文件夹（PCP、TMP、SLR、HWD、WND）中

PS：由于要生成SLR太阳辐射，但是原始数据中只有日照时数，我们是用代码通过日照时数对SLR进行计算的，然后再生成SLR的txt，详细情况见代码。另外空值和缺失值自动换为-99，该代码还可以生成质量报告，质量报告能告诉你各个站点的那天存在异常值，告诉你替换的是哪些天的哪个值。

以下为代码：

import os
import pandas as pd
import math  # 算太阳辐射
import numpy as np  # === 新增：用于判定 NaN

# 定义文件夹路径
folder_path = r'D:\你存放各站点excel的文件夹位置'
# ===【输出文件夹设置】===
outdir = r"D:\你输出txt结果的文件夹位置"
os.makedirs(outdir, exist_ok=True)

# === 过滤掉 ~$.xlsx 锁文件，避免 PermissionError
excel_files = [
    os.path.join(folder_path, f)
    for f in os.listdir(folder_path)
    if f.endswith('.xlsx') and not f.startswith('~$')
]

# === 质控阈值（可按需调整）
QC = {
    'pcp_min': 0.0,     # mm/d，降水不能为负
    'pcp_max': 500.0,   # 极端上限；若你想更宽松可调大
    't_min':  -60.0,    # ℃
    't_max':   60.0,    # ℃
    'wnd_min': 0.0,     # m/s
    'wnd_max': 60.0,    # m/s（极端强风上限）
    'hmd_min': 0.0,     # %
    'hmd_max': 100.0,   # %
    'ssd_min': 0.0,     # h
    'ssd_max': 24.0,    # h
    # 兼容你原先“>1000判异常”的规则（继续保留）
    'HUGE': 1000.0
}

# 遍历每个Excel文件
for excel_file in excel_files:
    excel_filename = os.path.basename(excel_file)
    file_prefix = os.path.splitext(excel_filename)[0]  # 去掉扩展名，直接用整名

    # === 新增：每站点一份质控日志，收集到列表，最后写 CSV
    qc_log = []  # 每条记录：{'station','date','var','orig','reason'}

    # 气象数据的起始日期
    start_date = '19940101'

    # 定义输出文件名
    pcp_output_path = os.path.join(outdir, f"PCP{file_prefix}.txt")
    tem_output_path = os.path.join(outdir, f"TEM{file_prefix}.txt")
    wnd_output_path = os.path.join(outdir, f"WND{file_prefix}.txt")
    hmd_output_path = os.path.join(outdir, f"HMD{file_prefix}.txt")
    slr_output_path = os.path.join(outdir, f"SLR{file_prefix}.txt")
    lst_output_path = [pcp_output_path, tem_output_path, wnd_output_path, hmd_output_path]

    # 读取Excel文件中的第一个工作表
    df = pd.read_excel(excel_file, sheet_name=0, engine='openpyxl')  # === 修改：显式指定引擎更稳

    # ===尝试获取“月”列用于更准确日期；没有则允许缺省
    month_col_exists = '月' in df.columns

    # 筛选起始年份
    for index, time in enumerate(df['年']):
        if time == int(start_date[:4]):
            start = index
            break
    else:
        raise ValueError(f"{excel_filename} 中找不到年份 {start_date[:4]} 的数据。")

    # === 新增：一个小函数，生成日志里用的日期字符串
    def date_str_at(offset_idx: int) -> str:
        """offset_idx 是从 start 开始的相对索引"""
        y = int(df['年'].iloc[start + offset_idx])
        d = int(df['日'].iloc[start + offset_idx])
        if month_col_exists:
            m = int(df['月'].iloc[start + offset_idx])
            return f"{y:04d}-{m:02d}-{d:02d}"
        else:
            # 没有“月”列时，尽量提示
            return f"{y:04d}-??-{d:02d}"

    # 取数据（单位换算后）
    pcp    = list(df['20-20时累计降水量(0.1mm)'][start:] / 10)   # mm
    temmax = list(df['日最高气温(0.1℃)'][start:] / 10)         # ℃
    temmin = list(df['日最低气温(0.1℃)'][start:] / 10)         # ℃
    w      = list(df['平均风速(0.1m/s)'][start:] / 10)          # m/s
    h      = list(df['平均相对湿度(1%)'][start:])               # %
    sd     = list(df['日照时数(0.1小时)'][start:] / 10)        # h
    year   = list(df['年'][start:])
    day    = list(df['日'][start:])

    # === 新增：一个统一判定函数（None/NaN/空 → 缺失；越界/不合理 → 异常），返回新值与日志
    def qc_value(value, varname, i, extra_check=None):
        """
        value: 换算后的值
        varname: 'pcp'|'tmax'|'tmin'|'wnd'|'hmd'|'ssd'
        i: 从 start 起的相对索引，用于定位日期
        extra_check: 额外规则函数，返回 (ok:bool, reason:str)
        """
        date_s = date_str_at(i)
        # 缺失（空单元格/NaN）
        if value is None or (isinstance(value, float) and (np.isnan(value))):
            qc_log.append({'station': file_prefix, 'date': date_s, 'var': varname,
                           'orig': 'NaN', 'reason': 'missing'})
            return -99

        # 原大阈值（兼容旧规则）
        if isinstance(value, (int, float)) and abs(value) > QC['HUGE']:
            qc_log.append({'station': file_prefix, 'date': date_s, 'var': varname,
                           'orig': value, 'reason': f'abs(value)>{QC["HUGE"]}'})
            return -99

        # 一般物理阈值
        if varname == 'pcp' and (value < QC['pcp_min'] or value > QC['pcp_max']):
            qc_log.append({'station': file_prefix, 'date': date_s, 'var': varname,
                           'orig': value, 'reason': f'out_of_range[{QC["pcp_min"]},{QC["pcp_max"]}]'})
            return -99
        if varname in ('tmax', 'tmin') and (value < QC['t_min'] or value > QC['t_max']):
            qc_log.append({'station': file_prefix, 'date': date_s, 'var': varname,
                           'orig': value, 'reason': f'out_of_range[{QC["t_min"]},{QC["t_max"]}]'})
            return -99
        if varname == 'wnd' and (value < QC['wnd_min'] or value > QC['wnd_max']):
            qc_log.append({'station': file_prefix, 'date': date_s, 'var': varname,
                           'orig': value, 'reason': f'out_of_range[{QC["wnd_min"]},{QC["wnd_max"]}]'})
            return -99
        if varname == 'hmd' and (value < QC['hmd_min'] or value > QC['hmd_max']):
            qc_log.append({'station': file_prefix, 'date': date_s, 'var': varname,
                           'orig': value, 'reason': f'out_of_range[{QC["hmd_min"]},{QC["hmd_max"]}]'})
            return -99
        if varname == 'ssd' and (value < QC['ssd_min'] or value > QC['ssd_max']):
            qc_log.append({'station': file_prefix, 'date': date_s, 'var': varname,
                           'orig': value, 'reason': f'out_of_range[{QC["ssd_min"]},{QC["ssd_max"]}]'})
            return -99

        # 额外规则（如 Tmax<Tmin）
        if extra_check is not None:
            ok, why = extra_check(value)
            if not ok:
                qc_log.append({'station': file_prefix, 'date': date_s, 'var': varname,
                               'orig': value, 'reason': why})
                return -99

        return value

    # === 修改：写 4 个基础要素时，加入缺失/异常判定与记录
    for i_path in lst_output_path:
        with open(i_path, 'w', encoding='utf-8') as file:
            file.write(f'{start_date}\n')

            if i_path == pcp_output_path:
                for i, v in enumerate(pcp):
                    newv = qc_value(v, 'pcp', i)
                    file.write(f"{newv}\n")

            elif i_path == tem_output_path:
                for i, (tx, tn) in enumerate(zip(temmax, temmin)):
                    # 先分别做基本 QC
                    tx_new = qc_value(tx, 'tmax', i)
                    tn_new = qc_value(tn, 'tmin', i)
                    # 再做 Tmax>=Tmin 的一致性检查
                    if (tx_new != -99) and (tn_new != -99) and (tx_new < tn_new):
                        # 视为异常 → 两个都置 -99（也可只置 Tmax）
                        qc_log.append({'station': file_prefix, 'date': date_str_at(i),
                                       'var': 'tmax,tmin', 'orig': f'{tx_new},{tn_new}', 'reason': 'tmax<tmin'})
                        tx_new, tn_new = -99, -99
                    file.write(f"{tx_new},{tn_new}\n")

            elif i_path == wnd_output_path:
                for i, v in enumerate(w):
                    newv = qc_value(v, 'wnd', i)
                    file.write(f"{newv}\n")

            elif i_path == hmd_output_path:
                for i, v in enumerate(h):
                    newv = qc_value(v, 'hmd', i)
                    file.write(f"{newv}\n")

    # === 原函数不变（仅在之后调用时做 SSD 的 QC）
    def Fsr(day, lat, SSD):
        gsc = 1361 * 0.0864 #太阳常数1361 W·m⁻²，乘0.0864将太阳常数换算成MJ/m²
        j = 2 * math.pi * (day - 1) / 365
        j2 = j * 2
        j3 = j * 3
        e0 = 1.00011 + 0.034221 * math.cos(j) + 0.00128 * math.sin(j) + 0.000719 * math.cos(j2) + 0.000077 * math.sin(j2)
        a1 = 0.399912 * math.cos(j)
        a2 = 0.070257 * math.sin(j)
        a3 = 0.006758 * math.cos(j2)
        a4 = 0.000907 * math.sin(j2)
        a5 = 0.002697 * math.cos(j3)
        a6 = 0.00148 * math.sin(j3)
        sita = (180 / math.pi) * (0.006918 - a1 + a2 - a3 + a4 - a5 + a6)
        ws1 = -math.tan((math.pi / 180) * sita) * math.tan((math.pi / 180) * lat)
        # === 新增：夹断，避免 acos 域外
        ws1 = max(-1.0, min(1.0, ws1))
        ws = 180 / math.pi * math.acos(ws1)
        h01 = math.cos((math.pi / 180) * lat) * math.cos((math.pi / 180) * sita) * math.sin((math.pi / 180) * ws)
        h02 = (math.pi / 180) * math.sin((math.pi / 180) * lat) * math.sin((math.pi / 180) * sita) * ws
        h0 = (1 / math.pi) * gsc * e0 * (h01 + h02)
        h1 = 0.8 * h0
        sl = 2 / 15 * ws
        if sl <= 0:
            return -99
        SLR = h1 * (0.248 + 0.752 * SSD / sl)
        if SLR > QC['HUGE']:
            return -99
        return SLR

    # 获取站点纬度lat
    file_path = r"F:\MMG hydrological situation\Data All\气象数据\老师给1994-2020\批量生成站点气象数据过程材料\SLRstation.txt"
    column_name_head = 'ID'
    colum_value_head = 'LAT'
    with open(file_path, 'r', encoding='utf-8') as file:  # === 修改：加 encoding 更安全
        header = file.readline().strip().split(',')
        column_name_index = header.index(column_name_head)
        column_value_index = header.index(colum_value_head)
        for line in file:
            columns = line.strip().split(',')
            if columns[column_name_index] == file_prefix:
                lat = float(columns[column_value_index])
                break

    # 写 SLR（太阳辐射），对 SSD 先做 QC
    with open(slr_output_path, 'w', encoding='utf-8') as file:
        file.write(f'{start_date}\n')
        for i in range(len(day)):
            if i > 0:
                if year[i] != year[i - 1]:
                    day[i] = 1
                else:
                    day[i] = day[i - 1] + 1

            # === 新增：对日照时数做 QC；若为 -99 则直接写 -99 并记录
            ssd_new = qc_value(sd[i], 'ssd', i)
            if ssd_new == -99:
                # 这里已经在 qc_value 里记过日志了
                file.write("-99\n")
                continue

            SLR = Fsr(day[i], lat, ssd_new)
            if SLR == -99:
                # 记录：通常是几何或>HUGE造成
                qc_log.append({'station': file_prefix, 'date': date_str_at(i), 'var': 'slr',
                               'orig': 'calc', 'reason': 'Fsr_invalid'})
            file.write(f"{SLR}\n")

    # === 新增：把本站点的质控日志写出
    qc_path = os.path.join(outdir, f"QC_{file_prefix}.csv")
    if len(qc_log) > 0:
        pd.DataFrame(qc_log).to_csv(qc_path, index=False, encoding='utf-8-sig')
    else:
        # 没有任何替换，仍然给一行提示
        pd.DataFrame([{'station': file_prefix, 'date': '', 'var': '', 'orig': '', 'reason': 'no_replacement'}]) \
          .to_csv(qc_path, index=False, encoding='utf-8-sig')

PS：气象数据制备的核心思路（怎么又要搜集日尺度的气象数据，又要弄气象发生器？算月尺度指标？）：

在 SWAT 制备气象数据的过程中，我需要五项气象数据，分别是日降水量、最高/最低气温、太阳辐射、相对湿度、平均风速，这五项都是按日分辨率获取的数据，将数据制作成 txt 文本，缺失值和异常值用 -99 代替，将这些数据输入到 SWAT 中。因为有缺失值和异常值，所以要用SWAT的气候发生器进行补充，气候发生器只在逐日驱动文件中“被标为 -99 的变量—日期对”上出手补值，其他仍使用观测值。将日分辨率的数据（平均气压、平均风速、平均气温、日最高气温、日最低气温、平均相对湿度、20-20 时降水量、小型蒸发量、大型蒸发量、日照时数）输入气候发生器得到月尺度的指标，基于多年逐日观测，计算各月统计量（如降水发生与转移概率、各月的 Tmax/Tmin、太阳辐射、风速、相对湿度的月均与标准差；如有需要，日照可用于反推辐射）。这些“逐月统计量”供气候发生器使用，并非直接作为逐日驱动。平均气压及小/大型蒸发并非 SWAT 必需项，也不是气候发生器的必要输入。然后补充参数库 SWAT2012.mdb 中的 WGEN_user，再在 SWAT 中调用参数库 SWAT2012.mdb 中的 WGEN_user。综上，输入气象数据实际上是上面两部分：① 主体是“逐日五要素 txt 驱动”；② 作为缺测补值依据的“气候发生器逐月统计参数库（WGEN_user）”。两者并行配合，互不替代。