Расчет суточных осадков ERA5 с использованием CDO

Question

По сути, это репост этого вопроса: https://confluence.ecmwf.int/pages/viewpage.action?pageId=149341027

Я загрузил ERA5 из CDS. Входной файл имеет 24 часовых шага (0, 1, 2, 3, 4, .., 23) для каждого календарного дня, начиная с 1 января по De c 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

Но я не уверен, что это правильно - есть предложения?

Answer 1

Для подобных проблем полезной командой в CDO является shifttime , которая, по сути, выполняет то, что написано на банке, и смещает отметку времени.

Проблема такого рода часто возникает с любым потоком или накопленным полем, где отметка времени, назначенная значению данных, указывает на END периода накопления времени или «окна» для Например, с 3-часовыми данными TRMM последние три часа дня имеют отметку 00 на дату после, и такие функции, как daymean или daysum, примененные напрямую, рассчитают среднее значение за 21 час за один день и 3 часа с предыдущего дня , неправильно. Сдвиг временной метки на три часа, чтобы время указывало на начало окна (или на 1,5, указывая на середину) перед выполнением вычисления, разрешит эту проблему.

Так что для вашего конкретного c вопроса где у вас есть длинный ряд почасовых данных с 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 )

Answer 2

# 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')