ERA5气象数据可以用来输入模型与研究目标建立关系,因此,在对于netcdf格式的ERA5数据进行处理的同时,进行归一化处理方便输入模型,代码如下:
import netCDF4 as nc
import glob
import os
import numpy as np
from osgeo import gdal, osr
from sklearn import preprocessing
# normalization for each array
def min_max_scaler(nc_array):
a = []
for j in range(len(nc_array[:])):
x = nc_array[j, :, :]
min_max_scaler = preprocessing.MinMaxScaler()
x1 = min_max_scaler.fit_transform(x)
a.append(x1)
a = np.array(a)
return a
# convert netcdf to geotiff
def nc2geotiff(nc_array, out_dir, out_basename, Num_Lon, Num_lat, LonMin, LatMax):
for i in range(len(nc_array[:])):
# create .tif
driver = gdal.GetDriverByName('GTiff')
out_tif_name = os.path.join(out_dir, out_basename + '_' + str(i+1) + '.tif')
print(out_tif_name)
out_tif = driver.Create(out_tif_name, Num_Lon, Num_lat, 1, gdal.GDT_Float32)
# set image boundary
geotransform = (LonMin, 0.25, 0, LatMax, 0, -0.25)
out_tif.SetGeoTransform(geotransform)
# set srs
srs = osr.SpatialReference()
srs.ImportFromEPSG(4326)
out_tif.SetProjection(srs.ExportToWkt())
# write to disk
out_tif.GetRasterBand(1).WriteArray(nc_array[i]) # write raster to memory
out_tif.FlushCache() # write to disk
out_tif = None # close .tif
out_dir = r'D:\DATA\ERA5\out'
# read in all netcdf files
nc_list = glob.glob(r'D:\DATA\ERA5\tmp\*.nc')
print(nc_list)
for nc_file in nc_list:
print('Reading ... ', nc_file)
# read the variable
nc_data_obj = nc.Dataset(nc_file)
nc_var = nc_data_obj.variables['t2m'] # change the variable name
nc_var_array = np.asarray(nc_var)
data = nc_var_array
# correction for netcdf packed data
scale_factor = nc_var.scale_factor
add_offset = nc_var.add_offset
FillValue = nc_var._FillValue
missing_value = nc_var.missing_value
data = nc_var_array[:]
data[data == FillValue] = np.nan
data[data == missing_value] = np.nan
data = data*scale_factor+add_offset
# normalization
print('Normalizing ...')
data = min_max_scaler(data)
# lat and lon
Lon = nc_data_obj.variables['longitude'][:]
Lat = nc_data_obj.variables['latitude'][:]
LonMin, LatMax, LonMax, LatMin = [Lon.min(), Lat.max(), Lon.max(), Lat.min()]
Num_Lat = len(Lat)
Num_Lon = len(Lon)
# define output file basename
basename = os.path.splitext(os.path.basename(nc_file))
out_basename = basename[0]
# convert nc to geotiff
nc2geotiff(data, out_dir, out_basename, Num_Lon, Num_Lat, LonMin, LatMax)
print('... ... ... COMPLETED ... ... ...')
由于ERA5数据有hourly的,需要计算daily的值,比如每日平均温度,每日平均风速等,以计算2020年每日2米温度(t2m)为例代码如下,同时需要注意的是如果是需要求出日均值再归一化代入模型的,将上一段代码中的归一化放入此段代码即可:
import glob
import os
import numpy as np
import pandas as pd
import rasterio
out_dir = r'D:\DATA\ERA5\t2m' # change the output directory
in_list = glob.glob(r'D:\DATA\ERA5\out\*.tif')
# normalization for each array
# def min_max_scaler(nc_array):
# a = []
# for j in range(len(nc_array[:])):
# x = nc_array[j, :, :]
# min_max_scaler = preprocessing.MinMaxScaler()
# x1 = min_max_scaler.fit_transform(x)
# a.append(x1)
# a = np.array(a)
# return a
# read in .tif files
def read_raster(file):
with rasterio.open(file) as src:
return (src.read(1))
# obtain metadata from one of the files
with rasterio.open(in_list[0]) as src:
meta = src.meta
meta.update(dtype=rasterio.float32)
# make a list for output file names
dates = pd.date_range('20200101', '20201231').strftime('%Y-%m-%d').to_list()
# caculate daily mean (every 24)
averaging_list = []
for i, j in zip(range(0, len(in_list), 24), dates):
print(i, j)
# make a list of every 24 files(of a single day)
averaging_list = in_list[i:i+24]
# read in all the files for a day and calculate the mean with nan ignored
array_list = [read_raster(x) for x in averaging_list]
array_mean = np.nanmean(array_list, axis=0)
# normalization
# print('Normalizing ...')
# array_mean = min_max_scaler(array_mean)
# make output file name
out_fn = os.path.join(out_dir, j + '.tif')
# write to disk
with rasterio.open(out_fn, 'w', **meta) as dst:
dst.write(array_mean.astype(rasterio.float32), 1)
print('... ... ... COMPLTED ... ... ...')
ERA5可以选择输出grib或者netcdf格式的文件,grib文件最简便的处理方法是在linux下通过CDO(Climate Data Operator)进行转换输出netcdf格式的文件:
安装CDO
sudo apt install cdo
安装完毕后使用命令进行转换
cdo -f nc copy file.grib file.nc
需要注意的是由grib格式转化成的netcdf格式是不含有scale factor, offset, fillvalue和missingvalue这些信息的,故在第一段代码中的correction部分仅对原生netcdf格式打包的数据有效。