Hi there,
Given several layers in a 3D NetCDF file (lon x lat x time), is there a way to compute a mean state of these layers … so a 2D ( lon x lat ) … without having to resort to loop?
Something like
grdmath FILE?var[a:b] MEAN
- OR -
grdmath FILE?var[list_indices] MEAN
Thanks
No.
Thanks for confirmation 
Loops are painful but only in shell(s)
A real case with Julia
julia> G = gmtread("MAOB/T2001_2019_SST_CLIMA.nc?z[0]");
julia> for k = 2:12 G += gmtread("MAOB/T2001_2019_SST_CLIMA.nc?z[@(k)]"); end
julia> G /= 12;
# Here some bug still be found-and-killed forced me to save grid in disk to visualize it
julia> gmtwrite("SST_avg.grd", G)
julia> imshow("SST_avg.grd", title="Average SST from 2001 to 2019", fmt=:png, colorbar=true)
Don’t forget grdmath’s stack operators; see grdmath -S for details on who they are.
So that would be something like
grdmath $file?var[a] $file?var[b] $file?var[c] ... -S MEAN = output.nc
can we pipe the list via cat list.txt | ?
sure, but like this: gmt grdmath $(cat list-of-grids) -S MEAN = output.nc
So …
cat << EOF > list.txt
grid?var[idt1]
grid?var[idt2]
grid?var[idtN]
EOF
grdmath $(cat list.txt) -S MEAN = output.nc
As simple as that ?
Sure, should be, caveat on keeping all those Unix special characters as they are, but otherwise should work.
There’s a third solution I routinely use for this, which is CDO’s operator timmean.
Given a stack of geographical regular latlon grids, one can
cdo timmean INFILE OUTFILE
GMT works pretty good in tandem with CDO
Hum, good to know, it can be useful in the case of continuous time-series.