xarray RMSE of each XY cell from the median

Question

I've written a leave one out spatial interpolation routine using xarray where an interpolation is done n-1 times across the input data. Each iteration of the interpolation is stored in an xarray dataset under a "IterN" dimension. Here's an example:

grids
Out[2]: 
<xarray.Dataset>
Dimensions:   (IterN: 84, northing: 53, easting: 54)
Coordinates:
  * easting   (easting) float64 5.648e+05 5.653e+05 ... 5.907e+05 5.912e+05
  * northing  (northing) float64 4.333e+06 4.334e+06 ... 4.359e+06 4.359e+06
  * IterN     (IterN) int64 0 1 2 3 4 5 6 7 8 9 ... 75 76 77 78 79 80 81 82 83
Data variables:
    Data      (IterN, northing, easting) float64 2.065 2.018 ... 0.3037 0.2913

I can compute some statistics over the dataset using the in-built functions, for example median here:

median = grids.median(dim="IterN")
median
Out[5]: 
<xarray.Dataset>
Dimensions:   (northing: 53, easting: 54)
Coordinates:
  * easting   (easting) float64 5.648e+05 5.653e+05 ... 5.907e+05 5.912e+05
  * northing  (northing) float64 4.333e+06 4.334e+06 ... 4.359e+06 4.359e+06
Data variables:
    Data      (northing, easting) float64 2.111 2.061 2.011 ... 0.3104 0.2992

I would like to compute the root mean square error of the iteration sets against something like the median above. I can't see an inbuilt function for this, and I'm struggling to see how to implement this answer.

The desired output is an RMSE grid in the same XY dimensions as the other grids.

Any help would be greatly appreciated!

Answer 1

Does the following go in the direction you mean?

median = grids.median(dim="IterN")
error = grids - median
rmse = np.sqrt((error**2).mean("IterN"))