python——GDAS1数据处理并转为netCDF格式

GDAS1数据来源介绍

美国国家环境预报中心（NCEP）使用全球资料同化系统（GDAS）将来自各种观测系统和仪器的数据，其中一个操作系统是GDAS（全球数据同化系统）。有关该模型的信息可在NCEP网站上找到。而在NOAA的空气资源实验室（ARL），NCEP模型输出用于空气质量传输和扩散建模。

GDAS数据每三小时存档一次，其中UTC时间00、06、12和18为模型分析数据，03、09、15和21为模型预测数据，该模型将处理数据从光谱系数形式转换为1度经纬度（360乘181）网格，模型输出采用GRIB格式。其中，ARL每天四次保存连续的分析和3小时的预测，以生成连续的数据存档。需要注意的是，一些字段（如降水量和地表通量）在分析时不可用，因此这些字段取自6小时预报文件。

GDAS1数据下载

一个比较简便的方法是通过ftp进行下载
其中ftp网址为：https://www.ready.noaa.gov/data/archives/gdas1/
该网址内存储了05年开始至今的数据

GDAS1数据命名

以 gdas1.nov22.w3 文件为例：
gdas1：表示数据类型
nov：为11月份英文缩写，表示11月
22：为2022年，年份缩写
w3：为一个月内15-21号
其中：w1至w5代表意义如下：
#=1 - days 1-7
#=2 - days 8-14
#=3 - days 15-21
#=4 - days 22-28
#=5 - days 29 - rest of the month

故一个数据文件中会存放几天的数据。

GDAS1数据说明

GDAS1中所包含的字段，缩写及单位如下：

其中：
Field：为气象字段名称
Units：为该变量单位
Label：为字段简写
Data Order：代表字段代号，其中S开头的为地面监测字段；U开头的为高空监测字段

GDAS1中包含的不同观测高度信息如下：

数据处理

在用python处理GDAS1时，发现该数据并不是标准的grib数据，故在github上找到了一个处理库包进行后续数据读取及写出等工作，参考网址为：https://github.com/martin-rdz/ARLreader
该库包目前是只能使用python36环境进行安装，并且可能因为自己的网不行，在用官方提供的安装方法时一直报错，后来的处理方式是下载zip文件后，利用python setup.py install进行离线安装。安装过程中的报错一般是因为库包版本不匹配导致，在新建的环境中进行重新安装即可。

本代码目的是根据选择文件，进行某一字段的均值合成，最后进行nc文件写出。
故代码分为两部分：

1. 原始单个文件内某一天数据读取，并计算该天内所需字段的均值
2. 计算完毕后进行结果输出，其中头文件信息不改变，输出结果为均值结果，输出格式为nc

python处理代码如下：

1. 原始数据内某一天均值合成

    def GDAS_alg_calc(d,level,field):
        data_day = []
        for h in np.arange(0, 24, 3):
            # 根据每三小时数据进行统计
            recinfo, grid, data = Ar.reader(file).load_heightlevel(
                d, h, level, field
            )
			
			# fc为-1时，代表数据不能使用
            if recinfo.fc != -1:
                data_day.append(data)

        # 所有数据都不能用的情况
        if data_day == []:
            continue

        # 根据算法求出当天结果
        data_out = np.mean(data_day, axis=0)

    return data_out

2. nc文件输出

'''
其中ref_file为原始文件
bname为输出文件名
field和alg_type根据上一步结果进行自定义，如RH2M和mean
'''
# 原始文件元信息读取
ref_gdas = Ar.reader(ref_file)
lats = ref_gdas.grid.lats
lons = ref_gdas.grid.lons

# nc文件输出
file_name = bname + ".nc"
ncfile = Dataset(file_name, mode="w", format="NETCDF4_CLASSIC")
lat_dim = ncfile.createDimension("lat", len(lats))
lon_dim = ncfile.createDimension("lon", len(lons))

ncfile.source = "GDAS"
ncfile.subtitle = bname
ncfile.griddef = ref_gdas.headerinfo["griddef"]
ncfile.Nx = ref_gdas.headerinfo["Nx"]
ncfile.Ny = ref_gdas.headerinfo["Ny"]
ncfile.resolution = ref_gdas.resolution

# 创建variables
lat = ncfile.createVariable("lat", np.float64, ("lat",))
lat.units = "degrees_north"
lat.long_name = "latitude"
lon = ncfile.createVariable("lon", np.float64, ("lon",))
lon.units = "degrees_east"
lon.long_name = "longitude"

locals()[f"{field}_{alg_type}"] = ncfile.createVariable(
    field, np.float64, ("lat", "lon")
)

# write data
nlats = len(lats)
nlons = len(lons)

lat[:] = lats
lon[:] = lons

data_arr = data.T
locals()[f"{field}_{alg_type}"][:, :] = data_arr

ncfile.close()

print(f"{file_name}输出完成")

注：GDAS1数据介绍详见：https://www.ready.noaa.gov/gdas1.php