Расчет суточного количества осадков ERA5 с использованием CDO
По сути, это репост этого вопроса: https://confluence.ecmwf.int/pages/viewpage.action?pageId=149341027
Я загрузил ERA5 с компакт-диска. Входной файл имеет 24-часовой шаг (0, 1, 2, 3, 4,...,23) для каждого календарного дня, начиная с 1 января по 31 декабря каждого рассматриваемого года.
ECMWF заявляет здесь https://confluence.ecmwf.int/display/CKB/ERA5%3A+How+to+calculate+daily+total+precipitation, что дневное общее количество осадков должно быть рассчитано путем накопления осадков, например, за 1 января 1979 г. путем суммирования шаги 1, 2,...,23 от 1 января И шаг 0 от 2 января. Это означает, что шаг 0 от 1 января 1979 года не включен в расчет общего количества осадков за этот день. Для расчета общего количества осадков за 2 января 1979 г. мы также используем шаги 1, 2, 3,...,23 этого дня плюс шаг 0 3 января и так далее.
Кажется, есть вариант сделать это в python следующим образом:
import xarray as xr # import xarray library
ds_nc = xr.open_dataset('name_of_your_file.nc') # read the file
daily_precipitation = ds_nc.tp.resample(time='24H').sum('time')*1000 # calculate sum with frequency of 24h and multiply by 1000
daily_precipitation.to_netcdf('daily_prec.nc') # save as netCDF
Теперь мне интересно, возможно ли это также с помощью операторов климатических данных (CDO) простым способом. Обычно я бы сделал любой такой расчет, используяdaysum в CDO, но я не уверен, что это правильно.
Кто-то предложил использовать:
cdo -f nc copy out.nc aux.nc
cdo -delete,timestep=1, aux.nc aux1.nc
cdo -b 32 timselsum,24 aux1.nc aux2.nc
cdo -expr,'ppt=tp*1000' -setmissval,-9999.9 -remapbil,r240x120 aux2.nc era5_ppt_prev-0_1979-2018.nc
Но я не уверен, что это правильно - есть предложения?
3 ответа
Для такого рода проблем полезной командой в CDO является shifttime, которая, по сути, выполняет то, что написано на банке, и сдвигает отметку времени.
Этот вид проблемы часто возникает с любым типом потока или накопленного поля, где отметка времени, назначенная значению данных, указывает на КОНЕЦ периода накопления времени или "окно", например, с 3-часовыми данными TRMM за последние три часа день имеет отметку 00 на дате после, и такие функции, как daymean или daysum, применяемые напрямую, неправильно вычисляют среднее значение 21 часа за один день и 3 часа с предыдущего дня. Сдвиг временной метки на три часа так, чтобы время указывало на начало окна (или, действительно, на 1,5, указывая на середину), перед выполнением вычисления решит эту проблему.
Поэтому для вашего конкретного вопроса, когда у вас есть длинная серия почасовых данных из ERA5 и вы хотите получить ежедневную сумму, вы можете сделать:
cdo shifttime,-1hour in.nc shift.nc # now step 0 on Jan 2 has Jan 1, 23:00 stamp
cdo daysum shift.nc daysum.nc
или соединены вместе:
cdo daysum -shifttime,-1hour in.nc daysum.nc
(ПРИМЕЧАНИЕ. Эта процедура отличается от тех, кто использует потоки из более ранней версии ERA-Interim, где потоки накапливаются в течение короткого периода прогноза. Для ERA5 "деаккумуляция" уже выполнена за вас. С ERA-Interim вам необходимо разница между последовательными временными шагами для преобразования из накопленного поля, и здесь есть сообщение, которое показывает, как это сделать с CDO или python: Лучшая дисперсия накопленных временных шагов netcdf с CDO)
# Correction to above python example to account for the time shift, as in the CDO example. Input file always needs to have the following day to the last day for which you want to compute daily sums/averages
import xarray as xr
ds_nc = xr.open_dataset('name_of_your_file.nc') # read the file
sds= ds_nc.shift(time=-1).dropna(dim='time',how='all') # shift to account for time shift for accumulated variables
daily_precipitation = sds.tp.resample(time='24H').sum('time')*1000 # calculate sum with frequency of 24h and multiply by 1000
# need to figure start_time and end_time for separately or slice differently.
sdaily=daily_precipitation.sel(time=slice("<start_time>", "<end_time>)") # drop the last value because values aren't complete.
sdaily.to_netcdf('daily_prec.nc')
Если вы отображаете данные ERA 5 за любые два дня, вы можете заметить, что tp в 0000 2 января (скажем) - это уже накопленные осадки за последние 24 часа (с 01:00 1 января до 24:00 (0000 2 января) 1 января. ). Таким образом, вам нужно выбрать только значения осадков на временном шаге 0000.