Misleading (?) averaging in the COP gallery example

[Matthew Mizielinski]

The COP-1d-plot example in the iris gallery reads in a series of cubes and then collapses them over latitude and longitude to get a time series using the iris.analysis.MEAN aggregator.

From what I understand of the aggregators this gives equal weights to every grid box rather than weighting by area. In some cases this will not make much difference to results (I have tested in several cases), but it does mislead the user as to what is going on. In other cases it can significantly alter results.

A simple example of what can go wrong;

import iris
import numpy as np
import iris.analysis.cartography

#build a simple cube with data = sin(latitude)**2
latitude = iris.coords.DimCoord(range(-90, 90, 5), standard_name='latitude',  units='degrees')
longitude = iris.coords.DimCoord(range(0, 360, 5), standard_name='longitude', units='degrees')

s = np.array([np.sin(latitude.points * np.pi/180.) ** 2,]*len(longitude.points)).transpose()

mycube = iris.cube.Cube(s, dim_coords_and_dims=[(latitude, 0), (longitude, 1)])

#cube must know about shape of globe (Earth radius supplied here) for area weighting
cs = iris.coord_systems.GeogCS(6371229)
for c in ['latitude', 'longitude']:
   mycube.coord(c).guess_bounds()
   mycube.coord(c).coord_system=cs

cube_weights = iris.analysis.cartography.area_weights(mycube)

print "grid box mean:      ", mycube.collapsed(['latitude','longitude'],iris.analysis.MEAN).data
print "area weighted mean: ", mycube.collapsed(['latitude','longitude'],iris.analysis.MEAN,weights=cube_weights).data

which gives the output

grid box mean:       [ 0.5]
area weighted mean:  [ 0.33290945]

i.e a significant difference (admittedly in a contrived situation to illustrate the point).

The documentation for iris.analysis.MEAN has also chosen the most simple example possible, zonal mean on a lat-lon grid, where area weights would not alter the result (although I would argue that it would be good practise to include them). A meridional mean example should be added here to emphasise the need for the area weights.

I've spoken to new users of iris who find this behaviour problematic, probably because some other analysis tools assume that MEAN == weighted MEAN when working with spatial dimensions and the example and documentation appears to be misleading.

I haven't considered what might happen when working with lat-lon-height data, but I suspect similar appropriate weights would be needed to account for different thicknesses of vertical levels.

Could the documentation and in particular the example be improved to reflect the methods required to extract meaningful data from ?