Repository navigation
cf.Field.collapse() fails when computing a volume mean where the volume cell measure is time-dependant #952
Description
Activity
I checked current
main, and the likely narrow spot is in the cell-measure axis reduction path:cf/weights.py:1533-1537buildsclm_axesby keeping only cell-measure axes whose data shape is greater than 1.- When
clm_axes != clm_axes0,cf/weights.py:1551-1553computesiaxesfrom the retained axes and passes them toclm.squeeze(...). - For a time-dependent volume measure with a singleton time axis, that looks like it may try to squeeze the retained non-singleton axes rather than the singleton axis that was removed from
clm_axes. - I did not see an open PR that appears to cover this issue.
A focused regression could use a volume cell measure with shape like
(1, z, y, x)and assert the singleton time axis is removed while the spatial axes remain in the weight key/data.Thanks for the report, @dlrhodson - I'm looking at your reproducer ...
OK - it's bug, which I've fixed:
In [2]: import cf ...: print(cf.__version__) ...: data=cf.read("test1.nc")[0] ...: print(data) ...: mean=data.collapse('volume: mean',weights='volume') 3.20.0 Field: sea_water_potential_temperature (ncvar%thetao) ----------------------------------------------------- Data : sea_water_potential_temperature(time(1), depth(5), Y(5), X(5)) degree_C Cell methods : time(1): mean Dimension coords: time(1) = [1989-06-01 00:00:00] 360_day : depth(5) = [0.5057600140571594, ..., 5.140361309051514] m : Y(5) = [0, ..., 4] : X(5) = [0, ..., 4] Auxiliary coords: latitude(Y(5), X(5)) = [[-84.21070861816406, ..., -83.79289245605469]] degrees_north : longitude(Y(5), X(5)) = [[72.5, ..., 76.5]] degrees_east Cell measures : measure:area(Y(5), X(5)) = [[16.0, ..., 16.0]] m2 : measure:volume(time(1), depth(5), Y(5), X(5)) = [[[[16.382505416870117, ..., 21.442960739135742]]]] m3 In [3]: print(mean) Field: sea_water_potential_temperature (ncvar%thetao) ----------------------------------------------------- Data : sea_water_potential_temperature(time(1), depth(1), Y(1), X(1)) degree_C Cell methods : time(1): mean depth(1): Y(1): X(1): mean Dimension coords: time(1) = [1989-06-01 00:00:00] 360_day : depth(1) = [2.91261887550354] m : Y(1) = [2.0] : X(1) = [2.0] In [4]: data.weights('volume', methods=True) Out[4]: {(1, 2, 3): 'volume cell measure'} In [5]: data.weights('volume') Out[5]: <CF Field: long_name=weights(depth(5), Y(5), X(5)) 1>
PR to follow, next release due in about 2 weeks (or you can pick it from main as soon as the PR is merged, etc)
Or in verbose mode:
In [5]: import cf ...: print(cf.__version__) ...: data=cf.read("test1.nc")[0] ...: mean=data.collapse('volume: mean',weights='volume', verbose=-1) 3.20.0 Can't create latitude and longitude coordinates: overwrite=False and latitude and/or longitude coordinates already exist: <CF AuxiliaryCoordinate: latitude(5, 5) degrees_north>, <CF AuxiliaryCoordinate: longitude(5, 5) degrees_east> all_methods, all_axes, all_within, all_over = ['mean'] [['domainaxis3', 'domainaxis2', 'domainaxis1']] [None] [None] axes = ['domainaxis3', 'domainaxis2', 'domainaxis1'] method = mean collapse_axes_all_sizes = Constructs: {'domainaxis1': <CF DomainAxis: size(5)>, 'domainaxis2': <CF DomainAxis: size(5)>, 'domainaxis3': <CF DomainAxis: size(5)>} collapse_axes = {'domainaxis1': <CF DomainAxis: size(5)>, 'domainaxis2': <CF DomainAxis: size(5)>, 'domainaxis3': <CF DomainAxis: size(5)>} collapse_axes_sizes = [5, 5, 5] Input weights = 'volume' Output weights = {(1, 2, 3): <CF Data(5, 5, 5): [[[16.382505416870117, ..., 21.442960739135742]]] 1>} Before collapse of data: iaxes, d_kwargs = [1, 2, 3] {'weights': {(1, 2, 3): <CF Data(5, 5, 5): [[[16.382505416870117, ..., 21.442960739135742]]] 1>}} f.shape = (1, 5, 5, 5) f.dtype = float32 After collapse of data: f.shape = (1, 1, 1, 1) f.dtype = float64 collapse_axes = {'domainaxis1': <CF DomainAxis: size(5)>, 'domainaxis2': <CF DomainAxis: size(5)>, 'domainaxis3': <CF DomainAxis: size(5)>} Removing cell_measure Changing axis size to 1: domainaxis1 Removing cell_measure Removing auxiliary_coordinate 'auxiliarycoordinate0' Removing auxiliary_coordinate 'auxiliarycoordinate1' Changing axis size to 1: domainaxis2 Changing axis size to 1: domainaxis3 Update cell methods: Original cell methods = {'cellmethod0': <CF CellMethod: domainaxis0: mean>} method = 'mean' within = None over = None Modified cell methods = Constructs: {'cellmethod0': <CF CellMethod: domainaxis0: mean>, 'cellmethod1': <CF CellMethod: domainaxis1: domainaxis2: domainaxis3: mean>} In [6]: print(mean) Field: sea_water_potential_temperature (ncvar%thetao) ----------------------------------------------------- Data : sea_water_potential_temperature(time(1), depth(1), Y(1), X(1)) degree_C Cell methods : time(1): mean depth(1): Y(1): X(1): mean Dimension coords: time(1) = [1989-06-01 00:00:00] 360_day : depth(1) = [2.91261887550354] m : Y(1) = [2.0] : X(1) = [2.0]
Reacted by Sadie L. Bartholomew and Dan Hodson@dlrhodson - is it OK if we add your test file to the cf-python test suite (and so will stored in this repo)?
Brilliant - thanks! Yes, please do add the test file to the test suite
Reacted by Sadie L. Bartholomew and David Hassell
Hi All,
If I try to compute a volume mean using cf.Field.collapse() where the cell measure is time dependant (sometimes true for ocean fields), cf.Field.collapse() throws an error.
This bug occurs due to a line in the
Weights.cell_meaure()function. The line:iaxes = [clm_axes0.index(axis) for axis in clm_axes]is intended to select redundant axes to be squeezed, but instead selects the non-redundant axes. This part of the code is only ever reached if the cell_measure has at least one axis of size 1. This occurs if cell_measure is time dependant and we are trying to collapse a field with a single time point.
What this line should do is to select all the axes in clm_axes0 that are not in clm_axes. This line:
iaxes = [clm_axes0.index(axis) for axis in (set(clm_axes0)-set(clm_axes))]does this. Replacing with this line fixes the bug.
Attached is a zip file with a demonstration of the bug.
cf_bug.zip