If I have a bunch of small tiles for a domain, but they don't necessarily cover the whole area. Is there a way of merging them together to make a single cube in Iris?
Short answer: no
Long answer: see this notebook :)
import iris
import iris.cube
import iris.coords
fname = iris.sample_data_path('E1_north_america.nc')
full_cube = iris.load_cube(fname)[0:2]
full_cube
Air Temperature (K) | time | latitude | longitude |
---|---|---|---|
Shape | 2 | 37 | 49 |
Dimension coordinates | |||
time | x | - | - |
latitude | - | x | - |
longitude | - | - | x |
Auxiliary coordinates | |||
forecast_period | x | - | - |
Scalar coordinates | |||
forecast_reference_time | 1859-09-01 06:00:00 | ||
height | 1.5 m | ||
Attributes | |||
Conventions | CF-1.5 | ||
Model scenario | E1 | ||
STASH | m01s03i236 | ||
source | Data from Met Office Unified Model 6.05 | ||
Cell methods | |||
mean | time (6 hour) |
# Setup the problem so that we have a bunch of chunks, but not necessarily complete ones...
chunks = [full_cube[:, 0:10, 0:20], full_cube[:, 2:10, 22:],
full_cube[:, 12:28, 0:10], full_cube[:, 10:30, 10:],
full_cube[:, 30:, 2:-3]]
import iris.plot as iplt
import matplotlib.pyplot as plt
Let's take a look at these individual chunks:
norm = plt.Normalize(vmin=full_cube.data.min(), vmax=full_cube.data.max())
for cube in chunks:
iplt.pcolormesh(cube[0], norm=norm)
plt.gca().coastlines()
plt.show()
/Users/pelson/dev/scitools/iris/lib/iris/coords.py:1000: UserWarning: Coordinate 'longitude' is not bounded, guessing contiguous bounds. 'contiguous bounds.'.format(self.name())) /Users/pelson/dev/scitools/iris/lib/iris/coords.py:1000: UserWarning: Coordinate 'latitude' is not bounded, guessing contiguous bounds. 'contiguous bounds.'.format(self.name()))
The problem boils down to us wanting this data in the form of a single cube...
One approach is to regrid each smaller cube into a single big one, and join those together:
all_lons = sorted(set().union(*[cube.coord('longitude').points for cube in chunks]))
all_lats = sorted(set().union(*[cube.coord('latitude').points for cube in chunks]))
print('Lon range: ', min(all_lons), max(all_lons))
print('Lat range: ', min(all_lats), max(all_lats))
Lon range: 225.0 315.0 Lat range: 15.0 60.0
import iris.analysis
sample_points = [('longitude', all_lons), ('latitude', all_lats)]
scheme = iris.analysis.Nearest(extrapolation_mode='mask')
resampled = [cube.interpolate(sample_points, scheme) for cube in chunks]
resampled[0]
Air Temperature (K) | time | latitude | longitude |
---|---|---|---|
Shape | 2 | 37 | 49 |
Dimension coordinates | |||
time | x | - | - |
latitude | - | x | - |
longitude | - | - | x |
Auxiliary coordinates | |||
forecast_period | x | - | - |
Scalar coordinates | |||
forecast_reference_time | 1859-09-01 06:00:00 | ||
height | 1.5 m | ||
Attributes | |||
Conventions | CF-1.5 | ||
Model scenario | E1 | ||
STASH | m01s03i236 | ||
source | Data from Met Office Unified Model 6.05 | ||
Cell methods | |||
mean | time (6 hour) |
iplt.pcolormesh(resampled[0][0])
plt.gca().coastlines()
plt.show()
/Users/pelson/dev/scitools/iris/lib/iris/coords.py:1000: UserWarning: Coordinate 'longitude' is not bounded, guessing contiguous bounds. 'contiguous bounds.'.format(self.name())) /Users/pelson/dev/scitools/iris/lib/iris/coords.py:1000: UserWarning: Coordinate 'latitude' is not bounded, guessing contiguous bounds. 'contiguous bounds.'.format(self.name()))
import numpy as np
mega_cube = resampled[0].copy()
for cube in resampled:
# Keep bolting on the non-masked parts of the resampled cubes
# until we have processed them all.
mega_cube.data = np.ma.where(
cube.data.mask, mega_cube.data, cube.data)
iplt.pcolormesh(mega_cube[0])
plt.gca().coastlines()
plt.show()
/Users/pelson/dev/scitools/iris/lib/iris/coords.py:1000: UserWarning: Coordinate 'longitude' is not bounded, guessing contiguous bounds. 'contiguous bounds.'.format(self.name())) /Users/pelson/dev/scitools/iris/lib/iris/coords.py:1000: UserWarning: Coordinate 'latitude' is not bounded, guessing contiguous bounds. 'contiguous bounds.'.format(self.name()))