In [15]:
import iris
import numpy
In [3]:
areacello_file = '/g/data/ua6/DRSv2/CMIP5/CSIRO-Mk3-6-0/historical/fx/ocean/r0i0p0/areacello/latest/areacello_fx_CSIRO-Mk3-6-0_historical_r0i0p0.nc'
data_file = '/g/data/ua6/DRSv2/CMIP5/CSIRO-Mk3-6-0/historical/mon/ocean/r1i1p1/hfds/latest/hfds_Omon_CSIRO-Mk3-6-0_historical_r1i1p1_185001-200512.nc'
In [4]:
area_cube = iris.load_cube(areacello_file)
In [6]:
print(area_cube)
In [7]:
data_cube = iris.load_cube(data_file)
In [8]:
print(data_cube)
In [12]:
if not data_cube.coord('latitude').has_bounds():
data_cube.coord('latitude').guess_bounds()
if not data_cube.coord('longitude').has_bounds():
data_cube.coord('longitude').guess_bounds()
area_weights = iris.analysis.cartography.area_weights(data_cube)
print(area_weights.shape)
print(type(area_weights))
In [16]:
area_weights = numpy.ma.masked_where(numpy.ma.getmask(data_cube.data), area_weights)
In [17]:
print(area_weights[0,50, 0:50])
In [18]:
print(area_cube.data[50, 0:50])
In [22]:
print(area_weights[0, 50, :].mean())
print(area_cube.data[50, 0:50].mean())
print(area_weights[0, 50, :].sum())
print(area_cube.data[50, :].sum())
In [23]:
def create_areacello_cube(cube):
"""Create an area cube."""
dim_coord_names = [coord.name() for coord in cube.dim_coords]
assert 'latitude' in dim_coord_names
assert 'longitude' in dim_coord_names
if not cube.coord('latitude').has_bounds():
cube.coord('latitude').guess_bounds()
if not cube.coord('longitude').has_bounds():
cube.coord('longitude').guess_bounds()
area_weights = iris.analysis.cartography.area_weights(cube)
out_dims = []
for index, name in enumerate(dim_coord_names):
out_dims.append((cube.coord(name), index))
area_cube = iris.cube.Cube(area_weights,
standard_name='cell_area',
long_name='Ocean Grid-Cell Area',
var_name='areacello',
units='m2',
dim_coords_and_dims=out_dims)
area_cube.data = numpy.ma.masked_where(numpy.ma.getmask(cube.data), area_cube.data)
if 'time' in dim_coord_names:
first_time = area_cube.coord('time').points[0]
area_cube = area_cube.extract(iris.Constraint(time=first_time))
return area_cube
In [24]:
new_area_cube = create_areacello_cube(data_cube)
In [25]:
print(new_area_cube.data[0, 50, :].mean())
print(area_cube.data[50, 0:50].mean())
print(new_area_cube.data[0, 50, :].sum())
print(area_cube.data[50, :].sum())
In [31]:
total_zonal_area = area_cube.collapsed('longitude', iris.analysis.SUM)
new_total_zonal_area = new_area_cube.collapsed('longitude', iris.analysis.SUM)
print(total_zonal_area.data[0:10])
print(new_total_zonal_area.data[0, 0:10])
In [ ]: