如何使用全年的小时数据计算每天的总降水量？答案

【问题标题】：How to calculate total precipitation per day using hourly data for whole year?如何使用全年的小时数据计算每天的总降水量？
【发布时间】：2019-04-16 14:59:10
【问题描述】：

我有来自 ERA5 的特定年份中每一天的每小时数据。我想将该数据从每小时转换为每天。我知道要做到这一点的漫长而艰难的道路，但我需要一些可以轻松做到的东西。

哥白尼在这里有一个代码https://confluence.ecmwf.int/display/CKB/ERA5%3A+How+to+calculate+daily+total+precipitation，如果数据集只转换一天，它可以正常工作，但是当转换一整年时，我遇到了问题。

链接下载 ERA5 数据集，可在 https://cds.climate.copernicus.eu/cdsapp#!/home 获得

按照这里的步骤使用哥白尼服务器

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

此脚本仅下载 2 天（2017 年 1 月 1 日和 2 日）的 houly 数据：

#!/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 个值。

输入数据集示例： 可从此处下载数据的子集（2017 年 1 月 1 日至 2 日）https://www.dropbox.com/sh/0vdfn20p355st3i/AABKYO4do_raGHC34VnsXGPqa?dl=0。只需在此之后使用第二个脚本来检查代码。 {全年代码为>10 GB，无法上传

提前致谢

【问题讨论】：

感谢代码转储，但它不是很有用。我无法运行您的代码，这意味着您只是粘贴了一本教科书，其中包含我们需要逐行了解您在做什么的内容。请编辑您的问题并给出示例输入数据和您期望从该数据中获得的输出。您的问题应该是一个完整 MINIMUM 和可验证 示例。我可以帮助你，但不能帮助你目前拥有的内容。
@Error-SyntacticalRemorse .....感谢您查看脚本，我现在对其进行了更改。简单来说，我想将全年的每小时数据（大约 8760 个值）转换为全年的每日数据（365 个值；24 小时的总和）。 Copernicus 服务器需要大量排队时间来下载大型数据集，但您可以非常轻松地更改第一个代码中的值。我希望这会有所帮助。
请发布我们可以使用的实际数据。该问题与以前的问题相同。您发布了指向何处获取它的链接，但我如何从那里获取它？只需获取您获取的数据的子集并将其添加到问题中即可。
@Error-SyntacticalRemorse....抱歉，我误解了您对评论的解释。数据的日落现在附在链接中。这可能会有所帮助。谢谢！！
今天晚些时候我会看看这个，看看我能做些什么。

标签： python-3.x dataset analysis cdo-climate era5

【解决方案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 或每年开始）等等。时间频率代码可以在Pandas docs中找到。

【讨论】：

不错的答案...只是一件事，我认为需要使用函数 xarray.Dataset.shift 在计算平均值或总和之前将时间提前 1 小时，因为时间戳位于每小时窗口的末尾（因此当天最后一小时的 23:00-24:00 在后一天的时间戳为 00:00）。请参阅此处：confluence.ecmwf.int/display/CKB/… 我会尝试在 include this 中编辑答案，但我不熟悉 xarray，可能会出错。

【解决方案2】：

从命令行使用 CDO 的另一种快速方法是：

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

注意，根据这里的答案/讨论：Calculating ERA5 Daily Total Precipitation using CDO ERA5 每小时数据的时间步长位于每小时窗口的末尾，因此您需要在求和之前移动时间戳，我不确定 xarray 解决方案是否能处理这个问题。还要有毫米/天，我认为需要求和，而不是取平均值。

【讨论】：