Merge ERA5-Land and ERA5 time series

Question

I am merging some ERA5-Land and ERA5 (netCDF) time series:

1. extract individual ERA5 and ERA5-Land maps and apply scale and offset factors -- this by looping the following

for BAND in $(seq 1 8784) ;do
    gdal_translate -unscale -b $BAND era5.nc era5_band_${BAND}.nc
done

and likewise

for BAND in $(seq 1 8784) ;do
    gdal_translate -unscale -b $BAND era5land.nc era5land_band_${BAND}.nc
done

This, however, does not preserve the NETCDF_DIM_time stamp of each map.

2. merge corresponding maps in one netCDF by copying ERA5-Land over ERA5 data and resampling to the spatial resolution of ERA5-Land data -- this by looping over:

for BAND in $(seq 1 8784) ;do
    gdal_merge.py -ps 0.1 0.1 -o era5land_band_${BAND}_patched.nc era5_band_${BAND}.nc era5land_band_${BAND}.nc 
done

This produces thousands of files (as many as the time steps in the initial time series). A small subset list:

era5land_band_1_patched.nc
era5land_band_2_patched.nc
era5land_band_3_patched.nc
era5land_band_4_patched.nc
era5land_band_5_patched.nc
era5land_band_6_patched.nc
era5land_band_7_patched.nc
era5land_band_8_patched.nc
era5land_band_9_patched.nc

Finally, I would like to stack them back to a single netCDF file. I went on and updated the variable name and long_name using hints from

ncrename -v Band1,t2m era5land_band_1_patched.nc 
ncatted -O -a long_name,t2m,o,c,'2 metre temperature' era5land_band_1_patched.nc

Given both source data ERA5 and ERA5-Land, have the same temporal resolution and time-stamps, how is it recommended building a time series with the merged maps and the original time-stamps?

I would like to stay practical with GDAL and Bash (also due to the availability of GNU Parallel). Are there better alternatives (simplicity first, speed second)?

Answer 1

Quick UPDATE of the workflow!

I am merging some ERA5-Land and ERA5 (netCDF) time series using GDAL and Bash as follows:

1. apply scale and offset factors and extract individual ERA5 and ERA5-Land maps and -- this by looping the following

for BAND in $(seq 1 8784) ;do
    gdal_translate -ot Float32 -unscale -b $BAND era5.nc era5_band_${BAND}.nc
done

and likewise

for BAND in $(seq 1 8784) ;do
    gdal_translate -ot Float32 -unscale -b $BAND era5land.nc era5land_band_${BAND}.nc
done

This, however, does not preserve the NETCDF_DIM_time stamp of each map.

2. merge corresponding maps in one netCDF by copying ERA5-Land over ERA5 data and resample, at the same time, to the spatial resolution of ERA5-Land data -- this by looping over:

for BAND in $(seq 1 8784)
    gdal_merge.py \
      -ot Float32 \
      -ps 0.1 0.1 \
      -o era5land_band_${BAND}_patched.nc \
      era5_band_${BAND}.nc \
      era5land_band_${BAND}.nc 
done

This produces thousands of files (as many as the time steps in the initial time series). A small subset list:

era5land_band_1_patched.nc
era5land_band_2_patched.nc
era5land_band_3_patched.nc
era5land_band_4_patched.nc
era5land_band_5_patched.nc
era5land_band_6_patched.nc
era5land_band_7_patched.nc
era5land_band_8_patched.nc
era5land_band_9_patched.nc

Given both source data ERA5 and ERA5-Land have the same temporal resolution and time-stamps,

6. Compute effective scale and offset values based on all input maps via the following formula (Python syntax below)

scale_factor = (data_maximum - data_minimum) / (2 ** bits - 1)
add_offset = data_minimum

where minimum and maximum values are the overall minimum and maximum from all maps

7. Scale

# define $OUT_MAP, something like
OUT_MAP=$(basename $IN_MAP .nc)_scaled.nc
for IN_MAP in $(/usr/bin/ls era5land_band_*_patched.nc) ;do
  EXPRESSION="(A - $ADD_OFFSET) / $SCALE_FACTOR"
  gdal_calc.py \
      --format NETCDF \
      -A $IN_MAP \
      --outfile=$OUT_MAP \
      --calc="$EXPRESSION" \
      --type=Int16 \
      --quiet
  ncatted -O -h -a add_offset,Band1,c,f,$ADD_OFFSET $OUT_MAP
  ncatted -O -h -a scale_factor,Band1,c,f,$SCALE_FACTOR $OUT_MAP

done

also, create a function with the above and use parallel!

8. Time-stamp and add-metadata for each individual patched map

# define the $OUT, something like:
HOURS_SINCE=$(gdalinfo "$cds_era5_land_2m_temperature_${YEAR}.nc" |grep -m 3 'NETCDF_DIM_time' |tail -n 1 |cut -d '=' -f2)
for BAND in $(/usr/bin/ls era5land_band_*_patched_scaled.nc) ;do
  IN_TIME=$(basename $BAND .nc)_timedimension.nc
  OUT=$(basename $BAND .nc)_timestamped.nc
  ncecat -u time $1 $IN_TIME  # add time dimension
  ncap2 \
     -h \
     -s "time[time]=array($HOURS_SINCE,1,\$time)" \
     -s 'time@units="hours since 1900-01-01 00:00:00.0"' \
     -O $IN_TIME $OUT
  ncrename -v Band1,t2m $OUT
  ncatted -a long_name,t2m,o,c,'2 metre temperature' $OUT
  ncatted -a units,t2m,o,c,K $OUT
  HOURS_SINCE=$(($HOURS_SINCE + 1))
done

Get HOURS_SINCE from the original file, increase by 1 after each time-stamping iteration, for example:

  9. Finally, merge maps to a time series in a single netCDF file

cdo mergetime era5_and_land_t2m_${YEAR}band*patched_scaled_timestamped.nc era5_and_land_t2m${YEAR}.nc

The workflow doesn't seem elegant (mixture of GDAL, NCO and CDO tools).

- Are there better alternatives (simplicity first, speed second)?
- With CDO, an input and an output is always expected (?). How to edit netCDF attributes of a file using CDO tools in-place?
- If best to use CDO/NCO, how to repeat this exercise using exclusively CDO/NCO?


## Useful

Don't create .xml files when computing statistics!

```bash
GDAL_PAM_ENABLED=NO
export GDAL_PAM_ENABLED

and some helper function

function get_number_of_bands() {
    gdalinfo $1 \
    |grep NETCDF_DIM_time_DEF \
    |cut -d'=' -f2 \
    |tr -d '{' \
    |cut -d',' -f1 ;
}