利用Python对栅格数据进行EOF并输出nc文件（含EOF分解）

最近比较苦恼如何将Python处理好的多维数组重新保存为nc文件

现在我就将利用xarray库完成nc文件读写流程分享给大家

我们以计算气温的EOF为例，来教大家如何对NC文件进行读写操作

数据来源：NCEP/NCAR 的温度数据
下载链接：https://downloads.psl.noaa.gov/Datasets/ncep.reanalysis.derived/pressure/air.mon.mean.nc

首先第一步载入相关库：

import xarray as xr
import numpy as np
import eofs

EOF模块：
这一部分我就不详细介绍了，网上有很多文章供大家参考

def eofs_(values):
    solver = eofs.xarray.Eof(values) 
    EOFs = solver.eofs(neofs=10, eofscaling=2)
    pcs = solver.pcs(npcs=10, pcscaling=1)
    print(np.array(EOFs).shape)
    return EOFs, pcs

nc文件读取模块：

def read_file(path):
    ds = xr.open_dataset(path) #打开一个nc文件，注意：尽量不要是中文的容易报错
    values = ds.air  # 读取air这个变量
    lat = ds.lat # 读取lat 纬度变量
    lon = ds.lon # 读取lon 经度变量
    time = ds.time # 读取时间序列
    levels = ds.level # 读取层次信息
    return values, lat, lon, time,levels

那么我们应该如何查看一个nc文件所含有的变量名称呢？

很简单使用print(ds.keys)就能查看到一个nc文件下的所有变量名称

nc文件写入模块：

def create_file(outpath, values, lat, lon,levels):
    ds = xr.Dataset() #穿件一个Dataset
    ds['EOFS'] = (('models','levels','lat', 'lon'), values.data) #将生成的EOFs写入到Dataset中
    #注意，models、levels、lat、lon这几个名称务必与values中的每个纬度对应上否则可能会生成一些奇怪的东西或者直接报错
    ds.coords['lat'] = lat.data #写入纬度信息
    ds.coords['lon'] = lon.data #写入经度信息
    ds.coords['models'] = np.arange(1,11,1) # 我在EOF模块中设置的是取前10个模态，所以这里我是用1~10来代替
    ds.coords['levels'] = np.array(levels) #向Dataset中写入各层信息
    ds.to_netcdf(outpath) #保存为nc文件

完整代码如下：

import xarray as xr
import numpy as np
import eofs

'''
示例文件下载地址：https://downloads.psl.noaa.gov/Datasets/ncep.reanalysis.derived/pressure/air.mon.mean.nc
'''

def eofs_(values):
    solver = eofs.xarray.Eof(values)
    EOFs = solver.eofs(neofs=10, eofscaling=2)
    pcs = solver.pcs(npcs=10, pcscaling=1)
    print(np.array(EOFs).shape)
    return EOFs, pcs

def read_file(path):
    ds = xr.open_dataset(path)
    values = ds.air
    lat = ds.lat
    lon = ds.lon
    time = ds.time
    levels = ds.level
    return values, lat, lon, time,levels

def create_file(outpath, values, lat, lon,levels):
    ds = xr.Dataset()
    ds['EOFS'] = (('models','levels','lat', 'lon'), values.data)
    ds.coords['lat'] = lat.data 
    ds.coords['lon'] = lon.data
    ds.coords['models'] = np.arange(1,11,1)
    ds.coords['levels'] = np.array(levels)
    ds.to_netcdf(outpath)


if __name__ == '__main__':
    air, lat, lon, time,levels = read_file(r'D:/air.mon.mean.nc')
    EOFs, pcs = eofs_(air)
    outpath = r'D:/air.mon.mean.nc.eofs.nc'
    create_file(outpath, EOFs, lat, lon,levels)