Как рассчитать общее количество осадков за день, используя почасовые данные за весь год?

Question

У меня есть почасовые данные от ERA5 для каждого дня в определенном году. Я хочу конвертировать эти данные с почасовой на ежедневную. Я знаю длинный и сложный способ сделать это, но мне нужно что-то, что делает это легко.

У Коперника есть код для этого здесь https://confluence.ecmwf.int/display/CKB/ERA5%3A+How+to+calculate+daily+total+precipitation, который отлично работает, если набор данных конвертируется только за один день, но при конвертации на целый год у меня возникают проблемы с этим.

Ссылка для загрузки набора данных ERA5, доступного по адресу https://cds.climate.copernicus.eu/cdsapp

Следуйте инструкциям по использованию сервера copernicus здесь

https://confluence.ecmwf.int/display/CKB/How+to+download+ERA5

Этот скрипт загружает данные за час только за 2 дня (1 и 2 января 2017 года):

#!/usr/bin/env python
"""
Save as get-tp.py, then run "python get-tp.py".

Input file : None
Output file: tp_20170101-20170102.nc
"""
import cdsapi

c = cdsapi.Client()
r = c.retrieve(
    'reanalysis-era5-single-levels', {
            'variable'    : 'total_precipitation',
            'product_type': 'reanalysis',
            'year'        : '2017',
            'month'       : '01',
            'day'         : ['01', '02'],
            'time'        : [
                '00:00','01:00','02:00',
                '03:00','04:00','05:00',
                '06:00','07:00','08:00',
                '09:00','10:00','11:00',
                '12:00','13:00','14:00',
                '15:00','16:00','17:00',
                '18:00','19:00','20:00',
                '21:00','22:00','23:00'
            ],
            'format'      : 'netcdf'
    })
r.download('tp_20170101-20170102.nc')
## Add multiple days and multiple months to donload more data

Ниже скрипт создаст файл netCDF всего за один день

#!/usr/bin/env python
"""
Save as file calculate-daily-tp.py and run "python calculate-daily-tp.py".

Input file : tp_20170101-20170102.nc
Output file: daily-tp_20170101.nc
"""
import time, sys
from datetime import datetime, timedelta

from netCDF4 import Dataset, date2num, num2date
import numpy as np

day = 20170101
d = datetime.strptime(str(day), '%Y%m%d')
f_in = 'tp_%d-%s.nc' % (day, (d + timedelta(days = 1)).strftime('%Y%m%d'))
f_out = 'daily-tp_%d.nc' % day

time_needed = []
for i in range(1, 25):
    time_needed.append(d + timedelta(hours = i))

with Dataset(f_in) as ds_src:
    var_time = ds_src.variables['time']
    time_avail = num2date(var_time[:], var_time.units,
            calendar = var_time.calendar)

    indices = []
    for tm in time_needed:
        a = np.where(time_avail == tm)[0]
        if len(a) == 0:
            sys.stderr.write('Error: precipitation data is missing/incomplete - %s!\n'
                    % tm.strftime('%Y%m%d %H:%M:%S'))
            sys.exit(200)
        else:
            print('Found %s' % tm.strftime('%Y%m%d %H:%M:%S'))
            indices.append(a[0])

    var_tp = ds_src.variables['tp']
    tp_values_set = False
    for idx in indices:
        if not tp_values_set:
            data = var_tp[idx, :, :]
            tp_values_set = True
        else:
            data += var_tp[idx, :, :]

    with Dataset(f_out, mode = 'w', format = 'NETCDF3_64BIT_OFFSET') as ds_dest:
        # Dimensions
        for name in ['latitude', 'longitude']:
            dim_src = ds_src.dimensions[name]
            ds_dest.createDimension(name, dim_src.size)
            var_src = ds_src.variables[name]
            var_dest = ds_dest.createVariable(name, var_src.datatype, (name,))
            var_dest[:] = var_src[:]
            var_dest.setncattr('units', var_src.units)
            var_dest.setncattr('long_name', var_src.long_name)

        ds_dest.createDimension('time', None)
        var = ds_dest.createVariable('time', np.int32, ('time',))
        time_units = 'hours since 1900-01-01 00:00:00'
        time_cal = 'gregorian'
        var[:] = date2num([d], units = time_units, calendar = time_cal)
        var.setncattr('units', time_units)
        var.setncattr('long_name', 'time')
        var.setncattr('calendar', time_cal)

        # Variables
        var = ds_dest.createVariable(var_tp.name, np.double, var_tp.dimensions)
        var[0, :, :] = data
        var.setncattr('units', var_tp.units)
        var.setncattr('long_name', var_tp.long_name)

        # Attributes
        ds_dest.setncattr('Conventions', 'CF-1.6')
        ds_dest.setncattr('history', '%s %s'
                % (datetime.now().strftime('%Y-%m-%d %H:%M:%S'),
                ' '.join(time.tzname)))

        print('Done! Daily total precipitation saved in %s' % f_out)

Мне нужен код, который будет следовать тому же шагу, что и вышеприведенные данные, но при условии, что у меня есть входной файл с данными за один год и я преобразую его в данные за один год.

Результатом должны быть ежедневные значения для расчетной переменной (например, количества осадков и т. Д.) За весь год.

Пример: Допустим, у меня есть данные об осадках за весь год в виде 1 мм / ч каждый день, у меня будет 2928 значений за весь год.

То, что я хочу, - это 24 мм в день в течение всего года и только 365 значений для не високосного года.

Пример входного набора данных: Подмножество данных можно загрузить здесь (на 1 и 2 января 2017 г.) https://www.dropbox.com/sh/0vdfn20p355st3i/AABKYO4do_raGHC34VnsXGPqa?dl=0. Просто используйте второй сценарий после этого, чтобы проверить код. {код на весь год >10 ГБ, поэтому не может быть загружен

заранее спасибо

Answer 1

Xarray Resample это просто инструмент для вас. Он преобразует данные netCDF из одного временного разрешения (например, ежечасно) в другое (например, ежедневно) в одну строку. Используя ваш пример файла данных, мы можем создать ежедневные средства, используя следующий код:

import xarray as xr

ds = xr.open_dataset('./tp_20170101-20170102.nc')
tp = ds['tp'] # dimensions [time: 48, latitude: 721, longitude: 1440]
tp_daily = tp.resample(time='D').mean(dim='time') # dimensions (time: 2, latitude: 721, longitude: 1440)

Вы увидите, что resample команда принимает во временном коде, в этом случае 'D' что означает ежедневный, а затем мы указываем, что мы хотим вычислить среднее значение для каждого дня, используя почасовые данные этого дня с .mean(dim='time'),

Если вместо этого, например, вы хотите вычислить дневной максимум, а не дневное среднее, вы должны заменить .mean(dim='time') с .max(dim='time'), Вы также можете перейти с почасовой на ежемесячную (MS или месяц-начало), годовой (AS или годовой старт) и многое другое. Временные частотные коды можно найти в документах Панд.

Answer 2

Альтернативный быстрый метод из командной строки с использованием CDO:

cdo daysum -shifttime,-1hour era5_hourly.nc era5_daily.nc

Обратите внимание, что в соответствии с этим ответом / обсуждением здесь: Расчет суточного общего количества осадков ERA5 с использованием CDO почасовые данные ERA5 имеют временной шаг в конце часового окна, поэтому вам нужно сместить временную метку, прежде чем подсчитывать сумму, я не уверен, что Решение xarray справляется с этим. Также, чтобы иметь мм / день, я думаю, нужно суммировать, а не брать среднее.