Tips：
本期向大家分享SWAT模型的气象数据库的制备方法。最终要制作5项气象数据，分别是日降水量、最高/最低气温、太阳辐射、相对湿度、平均风速，每项气象数据要有①站点文件和②站点气象数据文件，都用txt文本文件来储存。站点文件是对气象站点基本信息的汇总，站点气象数据文件是对具体的气象数据的存储。一个站点文件往往对应多个站点气象数据文件。
最终得到的结果示例：（以日降水量为例）
①站点文件：
各字段分别是：序号ID、站点名称NAME、纬度LAT、经度LONG、高程ELEVATION，各字段用逗号分隔。
② 站点气象数据文件（列举某一站点的数据）：
注意数据必须是逐日数据，除了气温数据以外，其他数据都是只有一列（气温数据需要将日最高最低气温数据以逗号分隔，按两列放在txt文件里）。文件第一行是起始日期2014年1月1日，从第二行开始是逐日的降水量，如果有异常值，用-99来表示。每个站点都要制作一个站点气象数据文件，所以一个站点文件对应多个站点气象数据。
降水站点数据示例：
气温站点数据示例：
1 原始气象数据说明
①数据表格字段
原始气象数据表共有12个字段：站号、年份、月份、日、平均风速、日照时数、平均气温、日最高气温、日最低气温、平均水气压、平均相对湿度、20-20降水量。
②量纲（度量单位）
注意这里多是以0.1为单位的，但是SWAT读取的是以1为单位的，需要除以10。比如30mm的降水，在原始气象数据中是300，所以在制作站点气象数据文件时，要把300除以10。日降水量、最高/最低气温、相对湿度、平均风速分别对应“20-20降水量”、“日最高气温”、“日最低气温”、“平均相对湿度”、“平均风速”，其中除了平均相对湿度的单位不需要改变，其他变量都要除以10。太阳辐射需要后续借助Python用日照时数来计算。
“平均风速”量纲：0.1m/s
“日照时数”量纲：0.1小时
“平均气温”量纲：0.1℃
“日最高气温”量纲：0.1℃
“日最低气温”量纲：0.1℃
“平均水气压”量纲：0.1hPa
“平均相对湿度”量纲：1%
“20-20降水量”量纲：0.1mm
③原始数据表格缩略图
可以看到有些数值高的离谱，如32766、32700、30014等，这些其实都是异常值，后续制作站点气象数据文件要用-99代替，后续用天气发生器模拟。
2 站点文件制作
在网上查找气象站点的纬经度和高程，共有5项气象数据，分别是日降水量、最高/最低气温、太阳辐射、相对湿度、平均风速，所以只需要制作5个站点文件。先在Excel中填好数据，然后另存为CSV格式，再修改为txt格式。
3 代码：批量生成站点气象数据文件
①原始气象数据的存放方式
原始气象数据按站点以Excel格式存储在同一个文件夹里
②批量生成站点气象数据文件代码
太阳辐射的Python代码参考了这位博主的R语言代码：SWAT模型学习——气象数据制备过程记录_swat年尺度模拟过程-CSDN博客
使用说明：注意上述原始气象数据的存放方式，所有Excel格式的文件放在同一个文件夹，所以要循环读取该文件夹的Excel文件，初始定义的路径是文件夹路径。此外，要将太阳辐射的站点文件复制一份放在代码所在的文件夹，用于代码读取站点纬度，进而计算太阳辐射。
import os
import pandas as pd
import math #算太阳辐射
# 定义文件夹路径
folder_path = r'C:\你的路径\存放Excel格式原始气象数据的文件夹'
# 获取文件夹中所有Excel文件的路径
excel_files = [os.path.join(folder_path, f) for f in os.listdir(folder_path) if f.endswith('.xlsx')]
# 遍历每个Excel文件
for excel_file in excel_files:# 获取当前Excel文件的名称excel_filename = os.path.basename(excel_file)
# 拿到代码之后需要按照自己气象数据修改的部分====================================================================================# 截取前6个字符，假设格式是"[54118]啦啦.xlsx"，取出 "54118"用来命名区分file_prefix = excel_filename[1:6] # 修改索引以适应实际情况# 气象数据的起始日期start_date = '20140101'# 定义p数据的输出文件名，加上后缀 "pcp、tem、w、h"pcp_output_path = f"PCP{file_prefix}.txt"tem_output_path = f"TEM{file_prefix}.txt"w_output_path = f"W{file_prefix}.txt"h_output_path = f"H{file_prefix}.txt"sr_output_path = f"SR{file_prefix}.txt"lst_output_path = [pcp_output_path, tem_output_path, w_output_path, h_output_path]# 读取Excel文件中的第一个工作表df = pd.read_excel(excel_file, sheet_name=0) # sheet_name=0 表示第一个工作表# 筛选出2014年往后的数据for index,time in enumerate(df['年']):if time==int(start_date[:4]): # start_date[:4]表示年份start=index #找到2014年的索引，后续读取从2014年开始往后的数据breakpcp=list(df['20-20降水量'][start:]/10) #把整体除以10，得到mm单位的降水temmax=list(df['日最高气温'][start:]/10)temmin=list(df['日最低气温'][start:]/10)w=list(df['平均风速'][start:]/10)h=list(df['平均相对湿度'][start:]) sd=list(df['日照时数'][start:]/10)year=list(df['年'][start:]) # 用于后续计数计算太阳辐射day=list(df['日'][start:]) # 用于后续计数计算太阳辐射# # 将处理后的数据写入到动态命名的txt文件for i in lst_output_path:with open(i, 'w', encoding='utf-8') as file: #open表示读取，'w' 表示写入模式, 'utf-8' 表示编码格式file.write(f'{start_date}\n') # 写入起始日期if i==pcp_output_path:for value in pcp:if value>1000:# value=-99value=0file.write(f"{value}\n") # 写入每天的降水量elif i==tem_output_path:for max_temp, min_temp in zip(temmax, temmin):if abs(max_temp)>1000:max_temp=-99if abs(min_temp)>1000:min_temp=-99file.write(f"{max_temp},{min_temp}\n") # 写入每天的最高温度和最低温度，用逗号分隔！！elif i==w_output_path:for value in w:if value>1000:value=-99file.write(f"{value}\n") # 写入每天的风速elif i==h_output_path:for value in h:if value>100:value=-99file.write(f"{value}\n") # 写入每天的湿度# 定义函数 Fsr，输入为 day, lat 和 SSDdef Fsr(day, lat, SSD):# gsc 是太阳常数，单位为 W/m²gsc = 118.108# j 是年内的第几天与年份周期相关的角度参数j = 2 * math.pi * (day - 1) / 365 # (day - 1) 计算当前年份中的天数，/365 将其归一化为年j2 = j * 2 # j2 是 j 的两倍，后续用于进一步的计算j3 = j * 3 # j3 是 j 的三倍，后续用于进一步的计算# e0 是地球与太阳的平均距离因子的修正项（反映了地球与太阳的距离变化对辐射的影响）e0 = 1.00011 + 0.034221 * math.cos(j) + 0.00128 * math.sin(j) + 0.000719 * math.cos(j2) + 0.000077 * math.sin(j2)# a1, a2, a3, a4, a5, a6 是太阳辐射强度模型中的系数，分别用于计算太阳辐射的角度调整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 是太阳的倾斜角（太阳的垂直角度）sita = (180 / math.pi) * (0.006918 - a1 + a2 - a3 + a4 - a5 + a6)# ws1 是天文日照角度的初步计算，ws 表示太阳光照射的角度ws1 = -math.tan((math.pi / 180) * sita) * math.tan((math.pi / 180) * lat)# ws 是天文日照角度，单位为度，用来表示太阳在天空中的角度ws = 180 / math.pi * math.acos(ws1)# h01 和 h02 是计算日照强度时的两个部分，分别表示太阳辐射的不同组件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 是辐射强度的基础值，单位为 W/m²h0 = (1 / math.pi) * gsc * e0 * (h01 + h02)# h1 是修正后的辐射强度（考虑了地表反射等因素）h1 = 0.8 * h0 # 取 80% 的太阳辐射到达地表# sl 是日照的有效时长，与太阳高度角和昼夜长度相关sl = 2 / 15 * ws # 通过太阳的角度来估算有效日照时长# SR 是计算出来的太阳辐射强度，单位为 W/m²SR = h1 * (0.248 + 0.752 * SSD / sl) # 通过日照时数和有效日照时长的关系来调整辐射强度if SR>100: # 大于1000的值认为是异常值，用-99代替,后续用于处理异常值SR=-99return SR # 返回计算得到的太阳辐射强度# 获取站点纬度latfile_path = 'SRstation.txt' # 替换为文件路径column_name_head = 'NAME' # 对应站点字段的名称是NAMEcolum_value_head = 'LAT' # 对应站点纬度的字段名称是LATwith open(file_path, 'r') as file:# 读取第一行作为列名，假设列名以逗号分隔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]==f"SR{file_prefix}":lat = float(columns[column_value_index])breakwith open(sr_output_path, 'a', 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]=1else:day[i]=day[i-1]+1SR=Fsr(day[i], lat, sd[i])file.write(f"{SR}\n")
③代码运行结果
成功生成站点气象数据文件。5项气象数据，22个站点，所以共生成110个站点气象数据文件，这比手动制备节省了大量时间。
以上是全部内容，欢迎大家评论区留言，批评指正。
