Skip to first unread message
Joshua Dorrington
May 29, 2019, 5:27:42 AM5/29/19
to SciTools (iris, cartopy, cf_units, etc.) - https://github.com/scitools
I have two iris cubes that won't concatenate when I believe they should, F0 and F1:
>> print(F0)
Dimension coordinates:
realisation x -
forecast_period - x
Auxiliary coordinates:
hour - x
lead_time - x
month_number - x
season - x
season_year - x
Scalar coordinates:
latitude: 46.5 degrees, bound=(42.0, 51.0) degrees
longitude: 1.5 degrees, bound=(-5.0, 8.0) degrees
Attributes:
Conventions: CF-1.5
Cell methods:
mean: latitude, longitude
>> print(F1)
Dimension coordinates:
realisation x -
forecast_period - x
Auxiliary coordinates:
hour - x
lead_time - x
month_number - x
season - x
season_year - x
Scalar coordinates:
latitude: 46.5 degrees, bound=(42.0, 51.0) degrees
longitude: 1.5 degrees, bound=(-5.0, 8.0) degrees
Attributes:
Conventions: CF-1.5
Cell methods:
mean: latitude, longitude
F0 and F1 were loaded from two different files, representing forecasts at different lead times. The forecast_period coordinate contains the true datetime, while the auxilliary lead_time coordinate contains an offset from the initialisation date:
DimCoord([1999-12-01 12:00:00, 1999-12-02 12:00:00, 1999-12-04 12:00:00,
1999-12-05 12:00:00, 1999-12-06 12:00:00], standard_name='forecast_period', calendar='gregorian', long_name='time_since_forecast_initialisation', var_name='L')
>>print(F1.coord("forecast_period")[:5])
DimCoord([1999-12-03 12:00:00, 1999-12-10 12:00:00, 1999-12-17 12:00:00,
1999-12-24 12:00:00, 2000-01-07 12:00:00], standard_name='forecast_period', calendar='gregorian', long_name='time_since_forecast_initialisation', var_name='L')
>>print(F0.coord("lead_time")[:5])
AuxCoord(array([12.5, 13.5, 8.5, 9.5, 10.5]), standard_name=None, units=Unit('1'), var_name='lead_time')
>>print(F1.coord("lead_time")[:5])
AuxCoord(array([14.5, 14.5, 14.5, 14.5, 14.5]), standard_name=None, units=Unit('1'), var_name='lead_time')
I have run the following checks to try and identify the clash, with no success:
assert F0.dtype==F1.dtype
assert F0.units==F1.units
assert F0.metadata==F1.metadata
assert F0.attributes==F1.attributes
for coord1,coord2 in zip(F0.coords(),F1.coords()):
assert coord1.standard_name==coord2.standard_name
assert coord1.var_name==coord2.var_name
assert coord1.long_name==coord2.long_name
assert coord1.dtype==coord2.dtype
assert coord1.units==coord2.units
assert coord1.attributes==coord2.attributes
assert coord1.has_bounds()==coord2.has_bounds()
if coord1.has_bounds():
assert coord1.bounds.dtype==coord2.bounds.dtype
assert coord1.is_compatible(coord2)
#There should be one and only one concatenatable coordinate
concat_coords=0
for coord1,coord2 in zip(F0.dim_coords,F1.dim_coords):
assert coord1.circular==coord2.circular
overlap=np.isin(coord1.points,coord2.points).any()
same=np.isin(coord1.points,coord2.points).all()
if not (overlap or same):
concat_coords+=1
assert concat_coords==1
return True
check_cubes_compatible(F0,F1)
True
iris.cube.CubeList([F0,F1]).concatenate_cube()
Traceback (most recent call last):
File "<ipython-input-183-b4d9aa82cafd>", line 1, in <module>
iris.cube.CubeList([F0,F1]).concatenate_cube()
File "/home/josh/anaconda3/envs/main/lib/python3.7/site-packages/iris/cube.py", line 511, in concatenate_cube
raise iris.exceptions.ConcatenateError(msgs)
ConcatenateError: failed to concatenate into a single cube.
An unexpected problem prevented concatenation.
Expected only a single cube, found 2.
Simon Peatman
May 29, 2019, 6:20:27 AM5/29/19
to SciTools (iris, cartopy, cf_units, etc.) - https://github.com/scitools
Joshua,
I notice your forecast_period
coordinates overlap each other and I think this is the reason the Cubes
won't concatenate. You're asking Iris to interweave the two Cubes,
which the concatenate method isn't able to do as far as I know.
Unfortunately, this is one of those times when Iris gives a very unhelpful error message!
Simon
Joshua Dorrington
May 29, 2019, 6:31:21 AM5/29/19
to SciTools (iris, cartopy, cf_units, etc.) - https://github.com/scitools
Hi Simon,
Thanks for taking a look at this.
The forecast_periods overlap as you say, but there are no actually repeated dates. I thought iris would be able to deal with that. Is that not the case?
Joshua Dorrington
May 29, 2019, 6:42:54 AM5/29/19
to SciTools (iris, cartopy, cf_units, etc.) - https://github.com/scitools
Thanks to Simon's comment about iris being unable to interweave coordinates, I have found a hack, although it isn't very pretty:
for f in F0.slices_over("forecast_period"):
CL.append(f)
for f in F1.slices_over("forecast_period"):
CL.append(f)
CL.sort(key=lambda cube:cube.coord("forecast_period").points[0])
F=CL.merge_cube()
F.transpose()
print(F)
2 metre temperature / (K) (realisation: 10; forecast_period: 1607)
Dimension coordinates:
realisation x -
forecast_period - x
Auxiliary coordinates:
lead_time - x
month_number -
season_year - x
Scalar coordinates:
hour: 12
latitude: 46.5 degrees, bound=(42.0, 51.0) degrees
longitude: 1.5 degrees, bound=(-5.0, 8.0)
season: djf
Attributes:
Simon Peatman
May 29, 2019, 6:48:11 AM5/29/19
to SciTools (iris, cartopy, cf_units, etc.) - https://github.com/scitools
I was about to suggest using slices_over(), but once you've separated the Cubes along the concatenation coordinate like this then you don't need to sort them yourself. This should work:
iris.cube.CubeList(list(F0.slices_over('forecast_period')) + list(F1.slices_over('forecast_period'))).merge_cube()
Simon
Joshua Dorrington
May 29, 2019, 6:55:04 AM5/29/19
to SciTools (iris, cartopy, cf_units, etc.) - https://github.com/scitools
Ah great, thanks Simon!
Reply all
Reply to author
Forward
0 new messages