Commit 905a187a authored by Javier Vegas-Regidor's avatar Javier Vegas-Regidor
Browse files

Climpercent adapted to reading discretized variables

parent 7de0c4b8
......@@ -16,7 +16,7 @@ CON_FILES = /esnas/autosubmit/con_files/
# Diagnostics to run, space separated. You must provide for each one the name and the parameters (comma separated) or
# an alias defined in the ALIAS section (see more below). If you are using the diagnostics just to CMORize, leave it
# empty
DIAGS = discretize,atmos,sfcWind,500,0,40 climpercent,atmos,sfcWind,1984,1984,1
DIAGS = discretize,atmos,sfcWind,2000,0,40 climpercent,atmos,sfcWind,1984,1984,1
# DIAGS = OHC
# Frequency of the data you want to use by default. Some diagnostics do not use this value: i.e. monmean always stores
# its results at monthly frequency (obvious) and has a parameter to specify input's frequency.
......@@ -87,7 +87,7 @@ OCEAN_TIMESTEP = 6
# CHUNK_SIZE is the size of each data file, given in months
# CHUNKS is the number of chunks. You can specify less chunks than present on the experiment
EXPID = testing_erainterim
STARTDATES = 19840101
STARTDATES = 19911101 19921101
MEMBERS = 0
MEMBER_DIGITS = 1
CHUNK_SIZE = 1
......
......@@ -44,8 +44,9 @@ class ClimatologicalPercentile(Diagnostic):
self.start_year = start_year
self.end_year = end_year
self.forecast_month = forecast_month
self.cmor_var = data_manager.variable_list.get_variable(variable, silent=True)
self.distribution = None
self.leadtime_files = {}
def __eq__(self, other):
return self.domain == other.domain and self.variable == other.variable and \
......@@ -86,7 +87,6 @@ class ClimatologicalPercentile(Diagnostic):
return ['{0}{1:02}01'.format(year, self.forecast_month) for year in range(self.start_year, self.end_year+1)]
def request_data(self):
self.leadtime_files = {}
for startdate in self.requested_startdates():
if startdate not in self.leadtime_files:
self.leadtime_files[startdate] = {}
......@@ -95,7 +95,7 @@ class ClimatologicalPercentile(Diagnostic):
None, None, vartype=VariableType.STATISTIC)
def declare_data_generated(self):
var_name = '{0.variable}prct{0.start_year}{0.forecast_month}-{0.end_year}{0.forecast_month:2d}'.format(self)
var_name = '{0.variable}prct{0.start_year}{0.forecast_month}-{0.end_year}{0.forecast_month:02d}'.format(self)
self.percentiles_file = self.declare_chunk(self.domain, var_name, None, None, None,
frequency=Frequencies.climatology, vartype=VariableType.STATISTIC)
......@@ -104,6 +104,7 @@ class ClimatologicalPercentile(Diagnostic):
Runs the diagnostic
"""
iris.FUTURE.netcdf_promote = True
self._get_distribution()
percentile_values = self._calculate_percentiles()
self._save_results(percentile_values)
......@@ -111,11 +112,12 @@ class ClimatologicalPercentile(Diagnostic):
temp = TempFile.get()
percentile_coord = iris.coords.DimCoord(ClimatologicalPercentile.Percentiles, long_name='percentile')
results = iris.cube.CubeList()
for leadtime in percentile_values.keys():
result = iris.cube.Cube(percentile_values[leadtime], var_name='percent', units=self.units)
for leadtime, data in percentile_values.items():
result = iris.cube.Cube(data, var_name='percent',
units=self.distribution.coord('bin').units)
result.add_dim_coord(percentile_coord, 0)
result.add_dim_coord(self.lat_coord, 1)
result.add_dim_coord(self.lon_coord, 2)
result.add_dim_coord(self.distribution.coord('latitude'), 1)
result.add_dim_coord(self.distribution.coord('longitude'), 2)
result.add_aux_coord(iris.coords.AuxCoord(np.int8(leadtime), long_name='leadtime'))
results.append(result)
iris.FUTURE.netcdf_no_unlimited = True
......@@ -125,24 +127,27 @@ class ClimatologicalPercentile(Diagnostic):
def _calculate_percentiles(self):
Log.debug('Calculating percentiles')
percentiles = iris.cube.CubeList()
percentiles = {}
bins = self.distribution.coord('bin').points
def calculate(point_distribution):
cs = np.cumsum(point_distribution)
total = cs[-1]
percentile_values = ClimatologicalPercentile.Percentiles * total
index = np.searchsorted(cs, percentile_values)
return [(self._bins[i + 1] + self._bins[i]) / 2 for i in index]
return [bins[i] for i in index]
for leadtime_slice in self.distribution.slices_over('leadtime'):
percentiles.append(np.apply_along_axis(calculate, 0, leadtime_slice.data))
return percentiles.merge_cube()
leadtime = leadtime_slice.coord('leadtime').points[0]
percentiles[leadtime]=np.apply_along_axis(calculate, 0, leadtime_slice.data)
return percentiles
def _get_distribution(self):
for startdate in self.leadtime_files:
for startdate, startdate_file in self.leadtime_files.iteritems():
Log.info('Getting data for startdate {0}', startdate)
data_cube = iris.load_cube(startdate.local_file)
data_cube = iris.load_cube(startdate_file.local_file)
if self.distribution is None:
self.distribution = data_cube
else:
self.distribution += data_cube
#!/usr/bin/env bash
#SBATCH -n 1
#SBATCH -w gustafson
#SBATCH -n 8
#SBATCH --time 2:00:00
#SBATCH --error=job.%J.err
#SBATCH --output=job.%J.out
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment