[ERA5 재분석 자료] 2. 시간(hourly) 자료를 일(daily) 자료로 변환하고 저장하기

- 개요

저번 포스트에서는 API로 ERA5 재분석 자료를 다운로드 받는 법을 알아보았습니다.

[ERA5 재분석 자료] 1. API를 활용한 다운로드: https://climate-data-science.tistory.com/17

 

저는 일(daily) 자료로 분석할 일이 있는데 ERA5 재분석 자료는 시간(hourly) 자료입니다.

즉 시간 자료를 일 자료로 변환해야합니다.

 

이번 포스트에서는 시간 자료를 일 자료로 변환하여 nc파일로 저장하는 법을 알아보겠습니다.

 

- ERA5 자료 읽기

import os
import netCDF4 as nc
import numpy as np

"""
읽어야 할 파일 이름 얻기
제 경우 INPATH 폴더 내에 모든 파일이 ERA5 자료입니다.
본인 상황에 따라 이 부분을 수정해주세요.
"""
INPATH = './data/reanalysis/ERA5/' # 파일 경로 입력
fnlist = sorted(os.listdir(INPATH)) # 이름 정렬을 해주면 시간 순으로 정렬이 됨
files = [os.path.join(INPATH, fnlist[i]) for i in range(len(fnlist))]


"""
파일을 읽는 부분
geopotential의 연직 1번째 index를 저장하고자 합니다.
geopotential의 longname은 'z'이고
제가 받은 ERA5 파일의 연진 첫번째 층은 500 hPa입니다.
ERA5 자료의 시간은 UTC로 저는 한국 지역을 분석할 것이기 때문에
KST 시간이 필요합니다.
"""
combined_times = [] # 시간 부분을 저장
combined_data = [] # 변수 부분을 저장
var_name = 'z' # geopotential의 longname
for file in files:
    ds = nc.Dataset(file)
    time_var = ds.variables['time']
    times = nc.num2date(time_var[:], units=time_var.units)
    
    times_utc = pd.to_datetime([pd.Timestamp(t.isoformat()) for t in times])
    times_kst = times_utc + pd.Timedelta(hours=9) # UTC to KST 변환
    
    var_data = ds.variables[var_name][:]
    var_data = var_data[:, 1, :, :] # 1은 연직층의 index
    
    combined_times.extend(times_kst)
    combined_data.append(var_data)

    ds.close()


combined_times = np.array(combined_times)

# axis=0은 시간으로 시간을 기준으로 자료를 연결
combined_data = np.concatenate(combined_data, axis=0)

 

- 시간 자료를 일 자료로 변환

전 pandas 라이브러리를 이용해서 시간 자료를 일 자료로 변환했습니다.

import pandas as pd

"""
pd.DataFrame 선언시 index를 시간(datetime 타입)으로 설정해줘야 
resample 메서드를 사용할 수 있습니다.
geopotential의 일 자료는 24시간 평균이기 때문에 mean 메서드를 사용했습니다만
만약 강수 자료라든가 다른 변수라면 sum을 사용해야 할 수도 있습니다.
"""
df = pd.DataFrame(data=combined_data.reshape(len(combined_times), -1), index=combined_times)
df_daily = df.resample('D').mean() 

ds = nc.Dataset(files[0])
# lat, lon 읽기
lat = ds.variables['latitude'][:]
lon = ds.variables['longitude'][:]
lat_dim = ds.variables['latitude'].size
lon_dim = ds.variables['longitude'].size
ds.close()

# DataFrame을 (시간, 위도, 경도) array로 변환
daily_data = df_daily.values.reshape(-1, lat_dim, lon_dim)

 

- 일 자료를 nc파일로 저장

daily_data 변수를 nc 파일로 저장하면 끝납니다.

"""
저장할 파일 이름을 정의합니다.
varname은 앞서 'z'로 정의했습니다.
앞서 언급했듯 전 500 hPa 자료를 추출한 것이라 이름 끝에 500을 붙였습니다.
본인의 상황에 맞게 수정해주세요.
"""
output_fn = f'./data/reanalysis/ERA5_daily/{var_name}500.nc'

with nc.Dataset(output_fn, 'w', format='NETCDF4') as new_ds:
    # 변수의 차원(시간, 위도, 경도) 정의
    new_ds.createDimension('time', daily_data.shape[0])
    new_ds.createDimension('latitude', lat_dim)
    new_ds.createDimension('longitude', lon_dim)
    
    # 저장할 변수 정의
    times_var = new_ds.createVariable('time', 'f8', ('time',))
    lat_var = new_ds.createVariable('latitude', 'f4', ('latitude',))
    lon_var = new_ds.createVariable('longitude', 'f4', ('longitude',))
    data_var = new_ds.createVariable(var_name, 'f4', ('time', 'latitude', 'longitude',), zlib=True)
    
    # 시간 속성 정의 (처음 다운로드 받은 ERA5 자료 속성에 적혀있음)
    times_var.units = 'hours since 1900-01-01 00:00:00'
    times_var.calendar = 'gregorian'
    
    # 변수 저장
    times_var[:] = nc.date2num(df_daily.index.to_pydatetime(), units=times_var.units, calendar=times_var.calendar)
    lat_var[:] = lat
    lon_var[:] = lon
    data_var[:] = daily_data