Skip to content

Commit

Permalink
Add anomalies computation
Browse files Browse the repository at this point in the history
  • Loading branch information
Javier Vegas-Regidor committed Aug 20, 2019
1 parent a358175 commit bfe1456
Show file tree
Hide file tree
Showing 2 changed files with 152 additions and 24 deletions.
88 changes: 77 additions & 11 deletions esmvalcore/preprocessor/_time.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,7 @@

import cf_units
import iris
import iris.util
import iris.coord_categorisation
import numpy as np

Expand Down Expand Up @@ -411,26 +412,91 @@ def climate_statistics(cube, operator='mean', period='full'):
cube = cube.collapsed('time', operator_method)
return cube

clim_coord = 'clim_coord'
if period in ['daily', 'day']:
iris.coord_categorisation.add_day_of_year(
cube, 'time', name=clim_coord
)
iris.coord_categorisation.add_day_of_year(cube, 'time')
clim_coord = 'day_of_year'
elif period in ['monthly', 'month', 'mon']:
iris.coord_categorisation.add_month_number(
cube, 'time', name=clim_coord
)
iris.coord_categorisation.add_month_number(cube, 'time')
clim_coord = 'month_number'
elif period in ['seasonal', 'season']:
iris.coord_categorisation.add_season(
cube, 'time', name=clim_coord
)
iris.coord_categorisation.add_season_number(cube, 'time')
iris.coord_categorisation.add_season(cube, 'time')
clim_coord = 'season_number'
else:
raise ValueError(
'Climate_statistics does not support period %s' % period
)
operator = get_iris_analysis_operation(operator)
cube = cube.aggregated_by(clim_coord, operator)
cube.remove_coord(clim_coord)
cube.remove_coord('time')
iris.util.promote_aux_coord_to_dim_coord(cube, clim_coord)
return cube


def anomalies(cube, period):
"""
Compute anomalies using a mean with the specified granularity.
Computes anomalies based on daily, monthly, seasonal or yearly means for
the full available period
Parameters
----------
cube: iris.cube.Cube
input cube.
period: str, optional
Period to compute the statistic over.
Available periods: 'full', 'season', 'seasonal', 'monthly', 'month',
'mon', 'daily', 'day'
Returns
-------
iris.cube.Cube
Monthly statistics cube
"""

reference = climate_statistics(cube, period=period)
if period in ['full']:
return cube - reference
elif period in ['daily', 'day']:
if not cube.coords('day_of_year'):
iris.coord_categorisation.add_day_of_year(cube, 'time')
if not reference.coords('day_of_year'):
iris.coord_categorisation.add_day_of_year(reference, 'time')
cube_coord = cube.coord('day_of_year')
ref_coord = reference.coord('day_of_year')
elif period in ['monthly', 'month', 'mon']:
if not cube.coords('month_number'):
iris.coord_categorisation.add_month_number(cube, 'time')
if not reference.coords('month_number'):
iris.coord_categorisation.add_month_number(reference, 'time')
cube_coord = cube.coord('month_number')
ref_coord = reference.coord('month_number')
elif period in ['seasonal', 'season']:
if not cube.coords('season_number'):
iris.coord_categorisation.add_season_number(cube, 'time')
if not reference.coords('season_number'):
iris.coord_categorisation.add_season_number(reference, 'time')
cube_coord = cube.coord('season_number')
ref_coord = reference.coord('season_number')
else:
raise ValueError('Period %s not supported')

data = cube.core_data()
cube_time = cube.coord('time')
ref = {}
for ref_slice in reference.slices_over(ref_coord):
ref[ref_slice.coord(ref_coord).points[0]] = np.ravel(
ref_slice.core_data())

for i in range(cube_time.shape[0]):
time = cube_time.points[i]
indexes = cube_time.points == time
indexes = np.ones_like(data, dtype=bool) * indexes
data[indexes] = data[indexes] - ref[cube_coord.points[i]]

cube = cube.copy(data)
return cube


Expand Down
88 changes: 75 additions & 13 deletions tests/unit/preprocessor/_time/test_time.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,7 +18,7 @@
regrid_time,
decadal_statistics, annual_statistics, seasonal_statistics,
monthly_statistics, daily_statistics,
climate_statistics
climate_statistics, anomalies
)


Expand Down Expand Up @@ -53,7 +53,6 @@ def setUp(self):
def test_get_january(self):
"""Test january extraction"""
sliced = extract_month(self.cube, 1)
print(sliced)
assert_array_equal(
np.array([1, 1]),
sliced.coord('month_number').points)
Expand All @@ -76,7 +75,6 @@ def setUp(self):
def test_extract_time(self):
"""Test extract_time."""
sliced = extract_time(self.cube, 1950, 1, 1, 1950, 12, 31)
print(sliced)
assert_array_equal(
np.arange(1, 13, 1),
sliced.coord('month_number').points)
Expand All @@ -91,15 +89,12 @@ def test_extract_time_one_time(self):
cube = _create_sample_cube()
cube = cube.collapsed('time', iris.analysis.MEAN)
sliced = extract_time(cube, 1950, 1, 1, 1952, 12, 31)
print(sliced)
assert_array_equal(np.array([360.]), sliced.coord('time').points)

def test_extract_time_no_time(self):
"""Test extract_time with no time step."""
cube = _create_sample_cube()[0]
sliced = extract_time(cube, 1950, 1, 1, 1950, 12, 31)
print('sliced', sliced, sliced.shape)
print('cube', cube, cube.shape)
assert cube == sliced


Expand All @@ -113,39 +108,34 @@ def setUp(self):
def test_get_djf(self):
"""Test function for winter"""
sliced = extract_season(self.cube, 'djf')
print(sliced)
assert_array_equal(
np.array([1, 2, 12, 1, 2, 12]),
sliced.coord('month_number').points)

def test_get_djf_caps(self):
"""Test function works when season specified in caps"""
sliced = extract_season(self.cube, 'DJF')
print(sliced)
assert_array_equal(
np.array([1, 2, 12, 1, 2, 12]),
sliced.coord('month_number').points)

def test_get_mam(self):
"""Test function for spring"""
sliced = extract_season(self.cube, 'mam')
print(sliced)
assert_array_equal(
np.array([3, 4, 5, 3, 4, 5]),
sliced.coord('month_number').points)

def test_get_jja(self):
"""Test function for summer"""
sliced = extract_season(self.cube, 'jja')
print(sliced)
assert_array_equal(
np.array([6, 7, 8, 6, 7, 8]),
sliced.coord('month_number').points)

def test_get_son(self):
"""Test function for summer"""
sliced = extract_season(self.cube, 'son')
print(sliced)
assert_array_equal(
np.array([9, 10, 11, 9, 10, 11]),
sliced.coord('month_number').points)
Expand Down Expand Up @@ -559,7 +549,6 @@ def setUp(self):
),
0,
)
print(self.cube_1)
add_auxiliary_coordinate([self.cube_1, self.cube_2])

def test_regrid_time_6hour(self):
Expand Down Expand Up @@ -607,7 +596,6 @@ def setUp(self):
),
0,
)
print(self.cube_1)
add_auxiliary_coordinate([self.cube_1, self.cube_2])

def test_regrid_time_3hour(self):
Expand Down Expand Up @@ -725,5 +713,79 @@ def get_decade(coord, value):
assert_array_equal(result.coord('time').points, expected_time)


def make_map_data(number_years=2):
"""Make a cube with time, lat and lon dimensions."""
times = np.arange(0.5, number_years * 360)
bounds = np.stack(((times - 0.5), (times + 0.5)), 1)
time = iris.coords.DimCoord(
times,
bounds=bounds,
standard_name='time',
units=Unit('days since 1950-01-01', calendar='360_day'))
lat = iris.coords.DimCoord(
range(2),
standard_name='latitude',
)
lon = iris.coords.DimCoord(
range(2),
standard_name='longitude',
)
data = np.array([[[0], [1], ], [[1], [0], ]]) * times
cube = iris.cube.Cube(
data,
dim_coords_and_dims=[(lon, 0), (lat, 1), (time, 2)]
)
return cube


@pytest.mark.parametrize('period', ['full', 'day', 'month', 'season'])
def test_anomalies(period):
cube = make_map_data(number_years=2)
result = anomalies(cube, period)
if period == 'full':
anom = np.arange(-359.5, 360, 1)
zeros = np.zeros_like(anom)
assert_array_equal(
result.data,
np.array([[zeros, anom], [anom, zeros]])
)
elif period == 'day':
anom = np.concatenate((np.ones(360) * -180, np.ones(360) * 180))
zeros = np.zeros_like(anom)
assert_array_equal(
result.data,
np.array([[zeros, anom], [anom, zeros]])
)
elif period == 'month':
anom1 = np.concatenate([np.arange(-194.5, -165) for x in range(12)])
anom2 = np.concatenate([np.arange(165.5, 195) for x in range(12)])
anom = np.concatenate((anom1, anom2))
zeros = np.zeros_like(anom)
print(result.data[0, 1])
assert_array_equal(
result.data,
np.array([[zeros, anom], [anom, zeros]])
)
elif period == 'season':
anom = np.concatenate((
np.arange(-314.5, -255),
np.arange(-224.5, -135),
np.arange(-224.5, -135),
np.arange(-224.5, -135),
np.arange(15.5, 105),
np.arange(135.5, 225),
np.arange(135.5, 225),
np.arange(135.5, 225),
np.arange(375.5, 405),
))
zeros = np.zeros_like(anom)
print(result.data[0, 1])
assert_array_equal(
result.data,
np.array([[zeros, anom], [anom, zeros]])
)
assert_array_equal(result.coord('time').points, cube.coord('time').points)


if __name__ == '__main__':
unittest.main()

0 comments on commit bfe1456

Please sign in to comment.