# Daily total rainfall from AUS2200 output

**URL:** <https://forum.access-hive.org.au/t/daily-total-rainfall-from-aus2200-output/2088>\
**Category:** AUS2200\
**Tags:** python, rainfall\
**Created:** [9 May 2024 06:18 UTC](https://forum.access-hive.org.au/t/daily-total-rainfall-from-aus2200-output/2088 "2024-05-09T06:18:01Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![cchambers](https://avatars.discourse-cdn.com/v4/letter/c/3ec8ea/32.png) [@cchambers](https://forum.access-hive.org.au/u/cchambers)\
**Post date:** [9 May 2024 06:18 UTC](https://forum.access-hive.org.au/t/daily-total-rainfall-from-aus2200-output/2088/1 "2024-05-09T06:18:01Z")

</div>

Hi AUS2200/ACCESS/Unified Model experts,

I’m trying to plot daily total accumulated rainfall from the AUS2200 model output. I think I need to sum the variable fld\_s04i201 “LARGE SCALE RAIN AMOUNT KG/M2/TS” over a day but I’m having a bit of trouble doing this and I’m not sure if this is the correct variable to be summing. Unlike most other variables in the output, fld\_s04i201 has a “time\_0 = 6” dimension which I think causes issues for me when working with the variable.

I’m wondering if anyone has made total rainfall calculations and if you would be willing to share your method here or any tips at all would be much appreciated. I am quite new at working with Python.

Thanks.  
Chris

---

<div class="post-metadata">

**Author:** ![Scott](https://sea2.discourse-cdn.com/flex020/user_avatar/forum.access-hive.org.au/scott/32/37_2.png) [@Scott](https://forum.access-hive.org.au/u/Scott)\
**Post date:** [9 May 2024 06:27 UTC](https://forum.access-hive.org.au/t/daily-total-rainfall-from-aus2200-output/2088/2 "2024-05-09T06:27:25Z")

</div>

Documentation for s4i201 (/g/data/access/projects/access/umdir/vn13.5/ctldata/STASHmaster/STASHmaster-meta.conf):

```none
[stashmaster:code(4201)]
description=LARGE SCALE RAIN AMOUNT KG/M2/TS
help=Large scale rain amount kg/m2/ts
    =This is the mass per metre squared of rain (liquid precipitation) which
    =falls on the surface during a single timestep as diagnosed by the large
    =scale precipitation scheme. It does not include a convective
    =contribution (see diagnostic 1 5 201).

```

You’ll need to look at what time processing of this field has already been done by the model - it’s probably already been accumulated over some interval.

Also consider field s5i214

```none
[stashmaster:code(5214)]
description=TOTAL RAINFALL RATE: LS+CONV KG/M2/S
help=Total rainfall rate at the surface in kg/m2/s.
    =This is the sum of the large scale and convective rainfall at the surface.

```

---

<div class="post-metadata">

**Author:** ![cchambers](https://avatars.discourse-cdn.com/v4/letter/c/3ec8ea/32.png) [@cchambers](https://forum.access-hive.org.au/u/cchambers)\
**Post date:** [10 May 2024 00:53 UTC](https://forum.access-hive.org.au/t/daily-total-rainfall-from-aus2200-output/2088/3 "2024-05-10T00:53:52Z")

</div>

Thanks for the suggestions Scott!

I’ve taken a look for the second rain variable (s5i214) but I can’t seem to find it in any of our output files.

Perhaps this is because it has not been output, possibly because our simulations are resolving the precipitation explicitly  
and do not use a cumulus scheme I think.

For finding the period of rainfall accumulation for variable s4i201,  
do you know where I would look to determine what time processing of this field has been done by the model?  
Is there some documentation that details this?

An example of our output files are at if anyone wants to take a look:  
/g/data/v45/cc6171/experiments/ECL-2016\_smallshift/netcdf/u-cs142/20160603T0000Z/aus2200/d0198/RA3/um

If anyone has ever plotted any kind of accumulated precipitation using ACCESS/UM/AUS2200 please any further tips would be much appreciated.

The three-hourly ACCESS accumulated precipitation is plotted by the BOM here [http://www.bom.gov.au/australia/charts/viewer/index.shtml](http://www.bom.gov.au/australia/charts/viewer/index.shtml)  
so I’m wondering if there might be a contact at the BOM who has done the processing for producing this.

---

<div class="post-metadata">

**Author:** ![Scott](https://sea2.discourse-cdn.com/flex020/user_avatar/forum.access-hive.org.au/scott/32/37_2.png) [@Scott](https://forum.access-hive.org.au/u/Scott)\
**Post date:** [10 May 2024 01:12 UTC](https://forum.access-hive.org.au/t/daily-total-rainfall-from-aus2200-output/2088/4 "2024-05-10T01:12:00Z")

</div>

Guessing from your output path your suite ID is u-cs142. You can get the details in the ‘STASH requests’ pages in rose edit.

LS rain is being output using time profile `10mA0-E`:

 ![image](https://us1.discourse-cdn.com/flex020/uploads/access1/original/2X/8/817a087c3818f49f5a23a57642310bccc6bb7462.png)

This is a time accumulation sampled every timestep accumulated over 10 minutes and output every 10 minutes:

 ![image](https://us1.discourse-cdn.com/flex020/uploads/access1/original/2X/9/9ef99de9d338a8f31c2585814e41416dd0cfadb0.png)

In that case, just sum the values over a model day to get the accumulated daily value.

---

<div class="post-metadata">

**Author:** ![cchambers](https://avatars.discourse-cdn.com/v4/letter/c/3ec8ea/32.png) [@cchambers](https://forum.access-hive.org.au/u/cchambers)\
**Post date:** [10 May 2024 03:43 UTC](https://forum.access-hive.org.au/t/daily-total-rainfall-from-aus2200-output/2088/5 "2024-05-10T03:43:11Z")

</div>

That is good to confirm that it is the 10-minute accumulated rainfall and I do just have to sum it up. So my problem is really just a technical and probably very simple one, to do with the “time\_0” dimension that this variable has. Most other variables simply use a “time” dimension but this one uses a “time\_0” dimension as below

```
float fld_s04i201(time_0, lat, lon) ;
	fld_s04i201:_FillValue = 1.e+20f ;
	fld_s04i201:standard_name = "stratiform_rainfall_amount" ;
	fld_s04i201:long_name = "LARGE SCALE RAIN AMOUNT KG/M2/TS" ;
	fld_s04i201:units = "kg m-2" ;
	fld_s04i201:um_stash_source = "m01s04i201" ;
	fld_s04i201:missing_value = 1.e+20f ;
	fld_s04i201:cell_methods = "time_0: mean" ;
	fld_s04i201:grid_mapping = "latitude_longitude" ;

```

I’m just wondering why does this variable have a “time\_0” dimension rather than just a “time” dimension? Because of this dimension, if I merge the files together, for example below using cdo, it only outputs the first 6 time steps. This command works fine for other variables such as radar, sea level pressure etc that have the “time” dimension.

cdo -L -O -mergetime -apply,-selname,fld\_s04i201 [umnsa/\*.nc] thisone.nc

I’m sure the solution is very simple but I have brain trouble trying to think about these things.

---

<div class="post-metadata">

**Author:** ![Scott](https://sea2.discourse-cdn.com/flex020/user_avatar/forum.access-hive.org.au/scott/32/37_2.png) [@Scott](https://forum.access-hive.org.au/u/Scott)\
**Post date:** [10 May 2024 03:51 UTC](https://forum.access-hive.org.au/t/daily-total-rainfall-from-aus2200-output/2088/6 "2024-05-10T03:51:04Z")

</div>

Probably there’s some difference between `time` and `time_0`, perhaps one has a timestep 0 output and not the other, or one is indicating the middle of the time interval instead of the end (e.g. [1200Z, 1210Z, 1220Z] vs. [1205Z, 1215Z, 1225Z] for your 10-minute field). The metadata of the time dimensions may help clarify.

I won’t be able to help with the CDO command sorry, hopefully someone else can chip in.

---

<div class="post-metadata">

**Author:** ![Aidan](https://sea2.discourse-cdn.com/flex020/user_avatar/forum.access-hive.org.au/aidan/32/42_2.png) [@Aidan](https://forum.access-hive.org.au/u/Aidan)\
**Post date:** [10 May 2024 05:02 UTC](https://forum.access-hive.org.au/t/daily-total-rainfall-from-aus2200-output/2088/7 "2024-05-10T05:02:41Z")

</div>

Some tips for getting assistance:

- Provide more information. You’ve said you’re interested in the `fld_s04i201` variable, but not said in which file to find this variable. The more work you ask someone to do the less likely you will get help.

- You can’t assume everyone can access your data files. Not everyone is a member of the same projects, in this case `v45`. If you want assistance from a wider pool of people then you can copy data files to a directory in `/scratch/public` with permissions for anyone to access it.

That said, I am a member of `v45`.

It appears the variable you’re interested in is in the `umnsa_spec_*.nc` files.

The `cdo` command you tried is accessing all the files in the directory. Is there a reason for that?

I’m not a `cdo` expert, and I am sure it can be done with that or `nco` and `ncrcat`. What you want can also be accomplished with python and `xarray`

```ipython
In [4]: import xarray as xr                                             
                                                                                                                                             
In [5]: ds = xr.open_mfdataset('umnsa_spec_*.nc')                                                                                            
                                                                        
In [6]: ds                                                              
Out[6]:                                                                 
<xarray.Dataset> Size: 1TB Dimensions: (time: 36, lat: 2120, lon: 2600, pseudo_level: 5,
                               lon_0: 2600, lat_0: 2121, time_0: 36, bnds: 2)
Coordinates:                                                            
  * time (time) datetime64[ns] 288B 2016-06-03T00:10:00 ...
  * lat (lat) float64 17kB -48.79 -48.77 ... -6.852 -6.832
  * lon (lon) float64 21kB 114.3 114.3 ... 165.7 165.7
  * pseudo_level (pseudo_level) int32 20B 1 2 3 4 5
  * lon_0 (lon_0) float64 21kB 114.3 114.3 ... 165.7 165.7  
  * lat_0 (lat_0) float64 17kB -48.8 -48.78 ... -6.822
  * time_0 (time_0) datetime64[ns] 288B 2016-06-03T00:05:0...
    theta_level_height float64 8B 5.0   
    model_theta_level_number int32 4B 1     
    sigma_theta float64 8B 0.9994
    height float64 8B 10.0
    height_0 float64 8B 1.5
Dimensions without coordinates: bnds
Data variables: (12/39)
    fld_s00i010 (time_0, time, lat, lon) float32 29GB dask.array<chunksize=(6, 1, 1060, 1300), meta=np.ndarray>
    latitude_longitude (time_0, time) float64 10kB -2.147e+09 ... -2.1...
    lat_bnds (time_0, time, lat, bnds) float64 44MB dask.array<chunksize=(6, 36, 2120, 2), meta=np.ndarray>
    lon_bnds (time_0, time, lon, bnds) float64 54MB dask.array<chunksize=(6, 36, 2600, 2), meta=np.ndarray>
    theta_level_height_bnds (time_0, time, bnds) float64 21kB dask.array<chunksize=(6, 36, 2), meta=np.ndarray>
    sigma_theta_bnds (time_0, time, bnds) float64 21kB dask.array<chunksize=(6, 36, 2), meta=np.ndarray>
    ... ...
    fld_s21i101 (time_0, time, lat, lon) float64 57GB dask.array<chunksize=(6, 1, 1060, 1300), meta=np.ndarray>
    fld_s21i104 (time, time_0, lat, lon) float32 29GB dask.array<chunksize=(36, 2, 1060, 1300), meta=np.ndarray>
    fld_s30i403 (time_0, time, lat, lon) float32 29GB dask.array<chunksize=(6, 1, 1060, 1300), meta=np.ndarray>
    fld_s30i404 (time_0, time, lat, lon) float32 29GB dask.array<chunksize=(6, 1, 1060, 1300), meta=np.ndarray>
    fld_s30i405 (time_0, time, lat, lon) float32 29GB dask.array<chunksize=(6, 1, 1060, 1300), meta=np.ndarray>
    fld_s30i406 (time_0, time, lat, lon) float32 29GB dask.array<chunksize=(6, 1, 1060, 1300), meta=np.ndarray>
Attributes:
    history: File /home/563/cc6171/cylc-run/u-cs142/share/cycle/20160603...
    Conventions: CF-1.6
    source: Data from Met Office Unified Model
    um_version: 12.2
In [7]: ds['fld_s04i201']                                                                                                                     
Out[7]:                                                                                                                                       
<xarray.DataArray 'fld_s04i201' (time: 36, time_0: 36, lat: 2120, lon: 2600)> Size: 29GB                                                      
dask.array<concatenate, shape=(36, 36, 2120, 2600), dtype=float32, chunksize=(36, 2, 1060, 1300), chunktype=numpy.ndarray>                    
Coordinates:                                                                                                                                  
  * time (time) datetime64[ns] 288B 2016-06-03T00:10:00 ...                                                              
  * lat (lat) float64 17kB -48.79 -48.77 ... -6.852 -6.832                                                              
  * lon (lon) float64 21kB 114.3 114.3 ... 165.7 165.7                                                                  
  * time_0 (time_0) datetime64[ns] 288B 2016-06-03T00:05:0...                                                              
    theta_level_height float64 8B 5.0                                                                                                  
    model_theta_level_number int32 4B 1                                                                                                      
    sigma_theta float64 8B 0.9994                                                                                               
    height float64 8B 10.0                                                                                                 
    height_0 float64 8B 1.5                                                                                                  
Attributes:                                                                                                                                   
    standard_name: stratiform_rainfall_amount                                                                                              
    long_name: LARGE SCALE RAIN AMOUNT KG/M2/TS                                                                                    
    units: kg m-2                                                                                                                  
    um_stash_source: m01s04i201                                                                                                              
    cell_methods: time_0: mean                                                                                                            
    grid_mapping: latitude_longitude                                                                                                      
                                                                                                                                              
In [8]: ds['time_0']                                                                                                                          
Out[8]:                                                                                                                                       
<xarray.DataArray 'time_0' (time_0: 36)> Size: 288B                                                                                           
array(['2016-06-03T00:05:00.000000000', '2016-06-03T00:15:00.000000000',
       '2016-06-03T00:24:59.999999744', '2016-06-03T00:35:00.000000256',
       '2016-06-03T00:45:00.000000000', '2016-06-03T00:54:59.999999744',
       '2016-06-03T01:05:00.000000256', '2016-06-03T01:15:00.000000000',
       '2016-06-03T01:25:00.000000000', '2016-06-03T01:35:00.000000000',
       '2016-06-03T01:45:00.000000000', '2016-06-03T01:54:59.999999744',
       '2016-06-03T02:05:00.000000256', '2016-06-03T02:15:00.000000000',
       '2016-06-03T02:24:59.999999744', '2016-06-03T02:35:00.000000256',
       '2016-06-03T02:45:00.000000000', '2016-06-03T02:55:00.000000000',
       '2016-06-03T03:05:00.000000000', '2016-06-03T03:15:00.000000000',
       '2016-06-03T03:24:59.999999744', '2016-06-03T03:35:00.000000256',
       '2016-06-03T03:45:00.000000000', '2016-06-03T03:54:59.999999744',
       '2016-06-03T04:05:00.000000256', '2016-06-03T04:15:00.000000000',
       '2016-06-03T04:25:00.000000000', '2016-06-03T04:35:00.000000000',
       '2016-06-03T04:45:00.000000000', '2016-06-03T04:54:59.999999744',
       '2016-06-03T05:05:00.000000256', '2016-06-03T05:15:00.000000000',
       '2016-06-03T05:24:59.999999744', '2016-06-03T05:35:00.000000256',
       '2016-06-03T05:45:00.000000000', '2016-06-03T05:55:00.000000000'],
      dtype='datetime64[ns]')
Coordinates:
  * time_0 (time_0) datetime64[ns] 288B 2016-06-03T00:05:0...
    theta_level_height float64 8B 5.0
    model_theta_level_number int32 4B 1
    sigma_theta float64 8B 0.9994
    height float64 8B 10.0
    height_0 float64 8B 1.5
Attributes:
    axis: T
    bounds: time_0_bnds
    standard_name: time

In [11]: ds['fld_s04i201'].resample(time_0='1D').sum()
Out[11]: 
<xarray.DataArray 'fld_s04i201' (time: 36, time_0: 1, lat: 2120, lon: 2600)> Size: 794MB
dask.array<transpose, shape=(36, 1, 2120, 2600), dtype=float32, chunksize=(36, 1, 1060, 1300), chunktype=numpy.ndarray>
Coordinates:
  * time (time) datetime64[ns] 288B 2016-06-03T00:10:00 ...
  * lat (lat) float64 17kB -48.79 -48.77 ... -6.852 -6.832
  * lon (lon) float64 21kB 114.3 114.3 ... 165.7 165.7
    theta_level_height float64 8B 5.0
    model_theta_level_number int32 4B 1
    sigma_theta float64 8B 0.9994
    height float64 8B 10.0
    height_0 float64 8B 1.5
  * time_0 (time_0) datetime64[ns] 8B 2016-06-03
Attributes:
    standard_name: stratiform_rainfall_amount
    long_name: LARGE SCALE RAIN AMOUNT KG/M2/TS
    units: kg m-2
    um_stash_source: m01s04i201
    cell_methods: time_0: mean
    grid_mapping: latitude_longitude

```

Note that the final object is a lazy `dask` array, and you need to add `.compute()` on the end, or evaluate it in some other way to force the computation. When I tried this on a login node it was killed, probably exceeded memory, so you’ll need to do this on an ARE session or similar.

---

<div class="post-metadata">

**Author:** ![cchambers](https://avatars.discourse-cdn.com/v4/letter/c/3ec8ea/32.png) [@cchambers](https://forum.access-hive.org.au/u/cchambers)\
**Post date:** [10 May 2024 05:45 UTC](https://forum.access-hive.org.au/t/daily-total-rainfall-from-aus2200-output/2088/8 "2024-05-10T05:45:01Z")

</div>

Thanks a lot for this python code Aidan! Could you tell me what this is calculating as I’m not too familiar with python, currently learning. Is this grabbing the 10 minute rainfall, summing it over the 6 times in each hourly file and then adding up to create a 24 hour total accumulated rain?

On the cdo command the reason for me accessing all the files in the directory is that I have temporarily copied the spec files to this other directory (for other reasons), so this command is only accessing those files.

As an update to the earlier cdo command to merge variables with time\_0 dimensions into 1 file then “merge” rather than “mergetime” seems to work e.g.:

cdo -L -O -merge -apply,-selname,fld\_s04i201 [umnsa/umnsa\_spec\*.nc] merge.nc

Using this merged file you can calculate the daily total rainfall using:

cdo daysum merge.nc daily.nc

---

<div class="post-metadata">

**Author:** ![Aidan](https://sea2.discourse-cdn.com/flex020/user_avatar/forum.access-hive.org.au/aidan/32/42_2.png) [@Aidan](https://forum.access-hive.org.au/u/Aidan)\
**Post date:** [10 May 2024 09:14 UTC](https://forum.access-hive.org.au/t/daily-total-rainfall-from-aus2200-output/2088/9 "2024-05-10T09:14:35Z")

</div>

> [@cchambers](#):
>
> Could you tell me what this is calculating as I’m not too familiar with python, currently learning. Is this grabbing the 10 minute rainfall, summing it over the 6 times in each hourly file and then adding up to create a 24 hour total accumulated rain?

In this case I’m using `xarray` to read in all the data into a single virtual dataset (called `ds` for convenience), so the `time-0` dimension has length 36, which was all the data that was available. If they’re 10 minutes apart then that is 6 hours, not 24.

So the python command below is accessing just the `fld_s04i201` variable by accessing it directly as a key value for the `ds` object, then calling the `.resample` method to convert to daily (`1D`) time intervals and then calling the `.sum()` as the aggregating function for the resample operation

```python
ds['fld_s04i201'].resample(time_0='1D').sum()

```

Here are some resources about the resample operation

> **[Xarray Interpolation, Groupby, Resample, Rolling, and Coarsen](https://earth-env-data-science.github.io/lectures/xarray/xarray-part2.html#resample)**
>
> In this lesson, we cover some more advanced aspects of xarray. Interpolation: In the previous lesson on Xarray, we learned how to select data based on its dimension coordinates and align data with ...

> **[Time series data](https://docs.xarray.dev/en/stable/user-guide/time-series.html#resampling-and-grouped-operations)**
>
> A major use case for xarray is multi-dimensional time-series data. Accordingly, we’ve copied many of features that make working with time-series data in pandas such a joy to xarray. In most cases, ...

> [@Efficient daily resampling using xarray](https://forum.access-hive.org.au/t/efficient-daily-resampling-using-xarray/233):
>
> Good day. I often resample hourly netcdf data (eg. ERA5) into daily means with xarray. However, when doing this for large 4D dataset (eg. ERA5 file) it is very slow. Does anyone have a more efficient method that is effective and returns an xarray object? An example of code below: import xarray as xr tfile = ‘/g/data/rt52/era5/pressure-levels/reanalysis/t/2022/t\_era5\_oper\_pl\_20220301-20220331.nc’ ft=xr.open\_dataset(tfile).resample(time=‘1D’).mean(keep\_attrs=True)

Glad you got the `cdo` version working. The advantage of `python` is that it doesn’t have to create lots of intermediate data, which wastes space and can rapidly fill directories with lots of derivative data.

Also using python in jupyter notebooks provides a nicer data exploration interface, and provides some record of you work, which can be useful to refer to later on.

---

<div class="post-metadata">

**Author:** ![cchambers](https://avatars.discourse-cdn.com/v4/letter/c/3ec8ea/32.png) [@cchambers](https://forum.access-hive.org.au/u/cchambers)\
**Post date:** [12 May 2024 02:40 UTC](https://forum.access-hive.org.au/t/daily-total-rainfall-from-aus2200-output/2088/10 "2024-05-12T02:40:52Z")

</div>

Thanks so much! This is hugely helpful to for both this total rainfall calculation and for my future data processing with python. I’ll be taking a thorough look through this next week.

A follow up note:

For calculating the total precipitation it is our understanding that the sum of the LARGE SCALE RAINFALL and the LARGE SCALE SNOW variables is appropriate, particularly for our wintertime case.

---

<div class="post-metadata">

**Author:** ![cchambers](https://avatars.discourse-cdn.com/v4/letter/c/3ec8ea/32.png) [@cchambers](https://forum.access-hive.org.au/u/cchambers)\
**Post date:** [14 May 2024 01:35 UTC](https://forum.access-hive.org.au/t/daily-total-rainfall-from-aus2200-output/2088/11 "2024-05-14T01:35:00Z")

</div>

As far as I can tell the python method using resample does not correctly calculate the accumulated rainfall. In the script below I am bringing in 4 datasets for four 6-hour periods and summing them to create a 24-hour total. The reason I am doing it this way because it prevents a kernel crash that occurs (even on XXLarge) if I bring them all in as one virtual dataset. Perhaps this can be solved by using dask in some way but I’m not familiar with that yet.

```py
import xarray as xr                                                                                                                                                                                         
#for a 24 hour segment
ds1 = xr.open_mfdataset('/g/data/v45/cc6171/experiments/ECL-2016_smallshift/netcdf/u-cs142/20160605T0000Z/aus2200/d0198/RA3/um/umnsa_spec_*.nc')
ds2 = xr.open_mfdataset('/g/data/v45/cc6171/experiments/ECL-2016_smallshift/netcdf/u-cs142/20160605T0600Z/aus2200/d0198/RA3/um/umnsa_spec_*.nc')
ds3 = xr.open_mfdataset('/g/data/v45/cc6171/experiments/ECL-2016_smallshift/netcdf/u-cs142/20160605T1200Z/aus2200/d0198/RA3/um/umnsa_spec_*.nc')
ds4 = xr.open_mfdataset('/g/data/v45/cc6171/experiments/ECL-2016_smallshift/netcdf/u-cs142/20160605T1800Z/aus2200/d0198/RA3/um/umnsa_spec_*.nc')
rain=ds1['fld_s04i201'].resample(time_0='1D').sum().values+ds2['fld_s04i201'].resample(time_0='1D').sum().values+ds3['fld_s04i201'].resample(time_0='1D').sum().values+ds4['fld_s04i201'].resample(time_0='1D').sum().values

```

The plot below is of this “rain” variable calculated using resample.  
 ![python_daily](https://us1.discourse-cdn.com/flex020/uploads/access1/original/2X/5/5adb6b315ac85ac3f2dd7fed6750a4dbaa7eaaa1.png)

and this plot is the same day accumulated rainfall calculated using CDO  
 ![cdo_daily](https://us1.discourse-cdn.com/flex020/uploads/access1/original/2X/a/a910d6a8bf977b753eba76b6d73fb6d306b3240f.png)

So it looks to me like the resample method may be summing a 1-hourly rainfall total from each 6 hourly period and so missing 5 hours of rainfall in each period. This is similar to what happened when I tried using nco. Probably this comes back to the time\_0 dimension used by this variable. I still have not found out why the rainfall uses this dimension rather than the “time” dimension used by other variables in the same file so if anyone knows why this is please share here. For the moment I will carry on using the CDO method. As a warning to others, if you are doing a similar calculation for total rainfall, based on my experience it is possible to mistakenly produce a total rainfall plot using python and nco that could be misinterpreted as the total rainfall but is not. My full script is below, the final section produces a plot from the CDO calculated total rainfall in a file “daily.nc”. This file was calculated using the CDO commands:  
cdo -L -O -merge -apply,-selname,fld\_s04i201 [umnsa/\*.nc] merge.nc  
cdo daysum merge.nc daily.nc

Finally I will say that this type of work is at the limit of my abilities so I apologize in advance for any mistakes.

```py
import xarray as xr     
import numpy as np
import pandas as pd
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.util import add_cyclic_point
import xesmf as xe
import datetime
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import cmocean.cm as cmo
import os, datetime
from sys import exit, argv
import netCDF4 as nc4
import geopandas as gpd
import regionmask
import matplotlib.pyplot as plt
import cartopy.io.shapereader as shpreader
from cdo import Cdo
cdo = Cdo()
%matplotlib inline
#for a 24 hour segment
ds1 = xr.open_mfdataset('/g/data/v45/cc6171/experiments/ECL-2016_smallshift/netcdf/u-cs142/20160605T0000Z/aus2200/d0198/RA3/um/umnsa_spec_*.nc')
ds2 = xr.open_mfdataset('/g/data/v45/cc6171/experiments/ECL-2016_smallshift/netcdf/u-cs142/20160605T0600Z/aus2200/d0198/RA3/um/umnsa_spec_*.nc')
ds3 = xr.open_mfdataset('/g/data/v45/cc6171/experiments/ECL-2016_smallshift/netcdf/u-cs142/20160605T1200Z/aus2200/d0198/RA3/um/umnsa_spec_*.nc')
ds4 = xr.open_mfdataset('/g/data/v45/cc6171/experiments/ECL-2016_smallshift/netcdf/u-cs142/20160605T1800Z/aus2200/d0198/RA3/um/umnsa_spec_*.nc')
#sum the four 6 hourly periods
rain=ds1['fld_s04i201'].resample(time_0='1D').sum().values+ds2['fld_s04i201'].resample(time_0='1D').sum().values+ds3['fld_s04i201'].resample(time_0='1D').sum().values+ds4['fld_s04i201'].resample(time_0='1D').sum().values

```

```py

x = ds1.lat
y = ds1.lon
vmn = 1
vmx = 100
# interval for color bands
interval = 10     
tcks = np.arange(vmn,vmx+interval,interval)
lvl = np.arange(vmn,vmx+interval,interval)
colormap = cm.hsv_r 
fig,ax=plt.subplots(1,1)
cp = ax.contourf(y,x,rain[0,0,:,:],cmap=colormap,levels=lvl,vmin=vmn,vmax=vmx,extend='max')
cbar_ax = fig.add_axes([0.25, 0.01, 0.5, 0.04])
cbar = fig.colorbar(cp, cax=cbar_ax, orientation='horizontal', label= 'rain', ticks=tcks, shrink=0.8)
# plot it
plt.show()

```

and finally for the CDO calculated rainfall plot

```py
#trying the same plot but using the cdo calculated 24 hour rainfall dataset
ds = xr.open_dataset('/g/data/v45/cc6171/experiments/ECL-2016_smallshift/spec/daily.nc')
rain = ds.fld_s04i201
x = ds.lat
y = ds.lon
vmn = 1
vmx = 100
# interval for color bands
interval = 10     
tcks = np.arange(vmn,vmx+interval,interval)
lvl = np.arange(vmn,vmx+interval,interval)
colormap = cm.hsv_r 
combines this with the terrain above
fig,ax=plt.subplots(1,1)
cp = ax.contourf(y,x,rain[2,:,:],cmap=colormap,levels=lvl,vmin=vmn,vmax=vmx,extend='max')
cbar_ax = fig.add_axes([0.25, 0.01, 0.5, 0.04])
cbar = fig.colorbar(cp, cax=cbar_ax, orientation='horizontal', label= 'rain', ticks=tcks, shrink=0.8)
# plot it
plt.show()

```

---

<div class="post-metadata">

**Author:** ![Scott](https://sea2.discourse-cdn.com/flex020/user_avatar/forum.access-hive.org.au/scott/32/37_2.png) [@Scott](https://forum.access-hive.org.au/u/Scott)\
**Post date:** [14 May 2024 03:01 UTC](https://forum.access-hive.org.au/t/daily-total-rainfall-from-aus2200-output/2088/12 "2024-05-14T03:01:50Z")

</div>

That is curious that the fields seem to have two time axes. Could you copy one of the netcdf files to /scratch/public so I can take a look?

---

<div class="post-metadata">

**Author:** ![cchambers](https://avatars.discourse-cdn.com/v4/letter/c/3ec8ea/32.png) [@cchambers](https://forum.access-hive.org.au/u/cchambers)\
**Post date:** [14 May 2024 04:04 UTC](https://forum.access-hive.org.au/t/daily-total-rainfall-from-aus2200-output/2088/13 "2024-05-14T04:04:08Z")

</div>

I’ve put one of the 1-hourly files that contains the 10-minute accumulated rainfall variable fld\_s04i201 in  
/scratch/public/cc6171  
it’s file:  
umnsa\_spec\_20160603T0600.nc

---

<div class="post-metadata">

**Author:** ![Scott](https://sea2.discourse-cdn.com/flex020/user_avatar/forum.access-hive.org.au/scott/32/37_2.png) [@Scott](https://forum.access-hive.org.au/u/Scott)\
**Post date:** [14 May 2024 04:35 UTC](https://forum.access-hive.org.au/t/daily-total-rainfall-from-aus2200-output/2088/14 "2024-05-14T04:35:57Z")

</div>

Looking at the metadata for the times, we have

```plaintext
dimensions:
        time = UNLIMITED ; // (6 currently)
        bnds = 2 ;
        time_0 = 6 ;
        #...
variables:
        double time(time) ;
                time:axis = "T" ;
                time:units = "days since 1970-01-01 00:00" ;
                time:standard_name = "time" ;
                time:calendar = "proleptic_gregorian" ;
        double time_0(time_0) ;
                time_0:axis = "T" ;
                time_0:bounds = "time_0_bnds" ;
                time_0:units = "days since 1970-01-01 00:00" ;
                time_0:standard_name = "time" ;
                time_0:calendar = "proleptic_gregorian" ;
        double time_0_bnds(time_0, bnds) ;

```

so `time` is an instantaneous time, while `time_0` covers some time interval, the boundaries of which are in `time_0_bnds`.

With a single file open it looks like fields are just 3-dimensional, e.g.

```plaintext
        float fld_s04i201(time_0, lat, lon) ;
                fld_s04i201:_FillValue = 1.e+20f ;
                fld_s04i201:standard_name = "stratiform_rainfall_amount" ;
                fld_s04i201:long_name = "LARGE SCALE RAIN AMOUNT KG/M2/TS" ;
                fld_s04i201:units = "kg m-2" ;
                fld_s04i201:um_stash_source = "m01s04i201" ;
                fld_s04i201:missing_value = 1.e+20f ;
                fld_s04i201:cell_methods = "time_0: mean" ;
                fld_s04i201:grid_mapping = "latitude_longitude" ;

```

you can see in cell\_methods that it’s claiming that there’s been a mean along time (potentially an error in the post-processing script, this is supposed to be a sum).

Comparing against a field using the `time` dimension you can see that there’s no time processing

```plaintext
        float fld_s09i203(time, lat, lon) ;
                fld_s09i203:_FillValue = 1.e+20f ;
                fld_s09i203:standard_name = "low_type_cloud_area_fraction" ;
                fld_s09i203:long_name = "LOW CLOUD AMOUNT" ;
                fld_s09i203:units = "1" ;
                fld_s09i203:um_stash_source = "m01s09i203" ;
                fld_s09i203:missing_value = 1.e+20f ;
                fld_s09i203:grid_mapping = "latitude_longitude" ;

```

If I open just that file with Xarray then it’s getting the correct axes - when Aidan opened all the files he was getting (time\_0, time, lat, lon) as the axes

```python
In [2]: xarray.open_dataset('umnsa_spec_20160603T0600.nc').fld_s04i201
Out[2]:<xarray.DataArray 'fld_s04i201' (time_0: 6, lat: 2120, lon: 2600)>

```

I wonder if the time axis names are getting swapped in some of the files? so some files might have the averaged time marked as time\_0 and others have it marked as time - that would confuse a lot of tools.

---

<div class="post-metadata">

**Author:** ![cbengel](https://avatars.discourse-cdn.com/v4/letter/c/ce7236/32.png) [@cbengel](https://forum.access-hive.org.au/u/cbengel)\
**Post date:** [14 May 2024 22:49 UTC](https://forum.access-hive.org.au/t/daily-total-rainfall-from-aus2200-output/2088/15 "2024-05-14T22:49:51Z")

</div>

Hi Scott, Chris and Aidan,

I have noticed strange header information in fields that have accumulation time-profiles. I got around it in simpler experiments by outputting the stratiform rainfall at each timestep and manually calculating the rainfall totals that I want.

This gives a bit more control. I know that the AUS2200 grids are very large, but given that it is a single layer this approach may work for you too.

---

<div class="post-metadata">

**Author:** ![Scott](https://sea2.discourse-cdn.com/flex020/user_avatar/forum.access-hive.org.au/scott/32/37_2.png) [@Scott](https://forum.access-hive.org.au/u/Scott)\
**Post date:** [15 May 2024 00:21 UTC](https://forum.access-hive.org.au/t/daily-total-rainfall-from-aus2200-output/2088/16 "2024-05-15T00:21:50Z")

</div>

Are you also converting the output using `um2netcdf_iris`? Any idea where that script is coming from? I’ve not been able to locate it.

---

<div class="post-metadata">

**Author:** ![cbengel](https://avatars.discourse-cdn.com/v4/letter/c/ce7236/32.png) [@cbengel](https://forum.access-hive.org.au/u/cbengel)\
**Post date:** [15 May 2024 00:26 UTC](https://forum.access-hive.org.au/t/daily-total-rainfall-from-aus2200-output/2088/17 "2024-05-15T00:26:43Z")

</div>

Scott, the strange header behaviour for skipped-accumulation-time-profiles actually begins in the fields files. Yes, I wonder what a netcdf conversion tool would do with those …

---

<div class="post-metadata">

**Author:** ![Scott](https://sea2.discourse-cdn.com/flex020/user_avatar/forum.access-hive.org.au/scott/32/37_2.png) [@Scott](https://forum.access-hive.org.au/u/Scott)\
**Post date:** [15 May 2024 00:32 UTC](https://forum.access-hive.org.au/t/daily-total-rainfall-from-aus2200-output/2088/18 "2024-05-15T00:32:30Z")

</div>

Interesting, both means and accumulations are marked in the um files with LBPROC=128

---

<div class="post-metadata">

**Author:** ![cbengel](https://avatars.discourse-cdn.com/v4/letter/c/ce7236/32.png) [@cbengel](https://forum.access-hive.org.au/u/cbengel)\
**Post date:** [15 May 2024 01:13 UTC](https://forum.access-hive.org.au/t/daily-total-rainfall-from-aus2200-output/2088/19 "2024-05-15T01:13:36Z")

</div>

The assigned timing has issues as well.

I requested every 10-minutes (accumulations) over 6 hours and in the fields file the timing indicated every 5-minutes for 3 hours.

It was not clear if the header information was off or the time profile …

Hence the swap to the per-timestep information.

---

<div class="post-metadata">

**Author:** ![dale.roberts](https://avatars.discourse-cdn.com/v4/letter/d/848f3c/32.png) [@dale.roberts](https://forum.access-hive.org.au/u/dale.roberts)\
**Post date:** [15 May 2024 01:51 UTC](https://forum.access-hive.org.au/t/daily-total-rainfall-from-aus2200-output/2088/20 "2024-05-15T01:51:11Z")

</div>

Hi @cchambers @Aidan @Scott I’m back from leave. The ‘fields with two time axes’ in Aidan’s example is an artifact of `xarray.open_mfdataset()`. When it does its concatenation, it seems to try its best to put all fields on consistent axes, even when that makes no sense, like in this case. These files are massive, so you really want to use `preprocess` (e.g. [More efficent use of xarray.open\_mfdataset() — CLEX CMS Blog](https://coecms.github.io/posts/2023-06-22-mfdataset-preprocess.html)) to pick out the `fld_s04i201` variable before xarray gets the chance to merge anything. That’ll make the rest of the job (the resample and sum) a whole lot easier.

Edit: URL updated

[Next page](https://forum.access-hive.org.au/t/daily-total-rainfall-from-aus2200-output/2088.md?page=2)