I do this by adding a new field, and getting the new field to work with the grid data. This works for Enzo, it's probably similar in other codes. Follow the examples on the website for adding new fields, my example is a little out of date.
#this is to make sure your output is the right shape
new_field = na.zeros(data['Density'].shape)
#watch these indices, I just made this up on the spot
new_field[1:,:,:] = 0.5*(data['Density'][:-1,:,:] + data['Density'][1:,:,:])
add_field('XAverageDensity', function = _XAverageDensity, validators=[ValidateSpatial(1,['Density'])])
The last bit about validators forces this function to work on the 3d data cube that yt gets off disk, rather than the flattened arrays that are better for variable resolution. The 1 indicates how many extra ghost zones yt gets, and ['Density'] can be a whole list of fields.