Thoughts on grdstack

Added by Paul 9 days ago

We are considering a new module called grdstack which will operate on a stack of 2-D grids (co-registered, perhaps). Action of interest are

  1. Computing statistics, e.g., mean, median, mode of all grids, standard deviation of grids
  2. Interpolate along the "time" axis, e.g., output equidistantly spaced grids or a 3-D netCDF cube.

For some of the basic operations where the time coordinate is not considered I wonder if it may be simpler to implement them in grdmath. E.g., maybe one could say

gmt grdmath *.grd SMEAN = mean.grd

where the SMEAN operator expects 2 or more grids (it would operate on all the items in the stack; hence the S prefix) and does the mean calculation. Or perhaps

gmt grdmath *.grd SMEANSTD = mean.grd = std.grd

which does both of the operations. Similarly for SMEDIAN, SMEDIANSCL, SMODE, AND SMODESCL. Could be extended to SMAX, SMIN and others as well. Then, grdstack could concern itself only with temporal interpolation operations.

Replies (2)

RE: Thoughts on grdstack - Added by Walter 7 days ago


This is a great idea. I asked Eric Leuliette about it, because he does a lot with this sort of thing.
Remko also did when he and I built the NIC 09 Ionosphere Climatology.

Here are comments from Eric L:
grdstack could potentially help us make some new products with the optimal interpolation work. I’ve been using the NetCDF Common Operators for some of the operations suggested for grdstack, and it would be better if a GMT-only script could be used instead.

For example, we could make climatology files by averaging over a stack of grids. Extending trend1d would also be useful to me. It would replace IDL or fortran code that I use instead, and it would be great to have all of the functionality in trend1d.
end quote.

What he means by trend1d (I added this suggestion when I forwarded your question to him) is, assuming each grid file is a geographical grid, and the sequence of files is a sequence of time samples, at each geographical grid point, fit a linear trend (possibly robustly with outlier detection) and create a grid with the amplitude of the linear trend coefficient (e.g. rate of sea level rise in mm/year from grids of sea height anomaly).

By "climatology" we mean things like min, max, average value at grid points.

It might also be helpful if there were some standard ways of filling "holes" where grids occasionally have NaNs. To this end, one might want to detect the case that (1) every grid has a NaN at certain geographical locations (which might mean the data are undefined/unmeasured at that location), versus (2) NaNs which show up once in awhile but not at every point in the time series. Then there could be options to interpolate or ignore the NaN in some way (such as fit the linear trend anyway, without a problem, but account for it in the trend slope uncertainty), or fill the hole with some sort of smooth surface spline.

We also have applications where a stack of grids might represent (lon,lat,time) of some process that varies in some complicated way, and we might have two such stacks, one representing measurements from satellite1 and the other from satellite2. For calibration purposes it might be nice to be able to fit a mean and time trend, and/or some low-order (very low order, like Degree 1) surface spherical harmonics, to the differences between stack1 and stack2. (Again, the surface spherical harmonics might need to handle NaNs over land, etc., so this is a least-squares type of estimator.)

Finally, I note that if averaging or fitting is involved it would be nice to have a switch to turn on/off that the grids are geographical, so that one could have spherical coordinate boundary conditions and also spherical patch or ellipsoidal patch area weighted averages and fits.

Just our ideas.


RE: Thoughts on grdstack - Added by Paul 7 days ago

Thanks Walter and Eric, very useful. FYI, over the weekend I added a new "reducing" mode to grdmath. This allows one of the reducing operators to apply to all items on the stack. For instance, if you want to compute the median value per pixel of a stack of grids called model_*.grd, you can now do

gmt grdmath model_*.grd -S MEAN = means.grd

The reducing operators are ADD, AND, MAD, LMSSCL, MAX, MEAN, MEDIAN, MIN, MODE, MUL, RMS, STD, SUB or XOR. So grdstack will focus on interpolations since plain stacking like above can be done in grdmath (i.e., there is no time involved). The -S option simply turns on the reducing mode which is in effect until the selected operator completes.