ERA5 气象数据处理:netCDF转geotiff同时归一化

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格式打包的数据有效。

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 194,088评论 5 459
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 81,715评论 2 371
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 141,361评论 0 319
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 52,099评论 1 263
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 60,987评论 4 355
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 46,063评论 1 272
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 36,486评论 3 381
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 35,175评论 0 253
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 39,440评论 1 290
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 34,518评论 2 309
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 36,305评论 1 326
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 32,190评论 3 312
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 37,550评论 3 298
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 28,880评论 0 17
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 30,152评论 1 250
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 41,451评论 2 341
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 40,637评论 2 335

推荐阅读更多精彩内容