Commit 67fa2416 authored by Javier Vegas-Regidor's avatar Javier Vegas-Regidor
Browse files

Reworked days_over chain to adapt again to datasets without members

parent f3782d5d
......@@ -6,7 +6,7 @@ SCRATCH_DIR = /scratch/Earth/$USER
# Root path for the cmorized data to use
DATA_DIR = /esnas:/esarchive
# Specify if your data is from an experiment (exp), observation (obs) or reconstructions (recon)
DATA_TYPE = exp
DATA_TYPE = recon
# CMORization type to use. Important also for THREDDS as it affects variable name conventions.
# Options: SPECS (default), PRIMAVERA, CMIP6
DATA_CONVENTION = SPECS
......@@ -16,8 +16,10 @@ 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 = daysover,atmos,sfcWind,1981,2012,11
# DIAGS = OHC
#DIAGS = discretize,atmos,sfcWind,,0,40
#DIAGS = climpercent,atmos,sfcWind,2010,2012,11
DIAGS = daysover,atmos,sfcWind,2010,2012,11
#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.
FREQUENCY = 6hr
......@@ -26,7 +28,7 @@ CDFTOOLS_PATH = ~jvegas/CDFTOOLS/bin
# If true, copies the mesh files regardless of presence in scratch dir
RESTORE_MESHES = False
# Limits the maximum amount of threads used. Default: 0 (no limitation, one per virtual core available)z
#MAX_CORES = 1
MAX_CORES = 2
[CMOR]
# If true, recreates CMOR files regardless of presence. Default = False
......@@ -70,7 +72,7 @@ SERVER_URL = https://earth.bsc.es/thredds
[EXPERIMENT]
# Experiments parameters as defined in CMOR standard
INSTITUTE = ecmwf
MODEL = system4_m1
MODEL = erainterim
# Model version: Available versions
MODEL_VERSION = Ec3.2_O1L75
# Atmospheric output timestep in hours
......@@ -87,7 +89,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 = 20131101 20141101 20151101
STARTDATES = 20101101 20111101 20121101
# STARTDATES = 19840101 19850101
MEMBERS = 0
MEMBER_DIGITS = 1
......
# coding=utf-8
from earthdiagnostics.diagnostic import *
from earthdiagnostics.utils import Utils, TempFile
import os
from earthdiagnostics.modelingrealm import ModelingRealm, ModelingRealms
import numpy as np
import os
from earthdiagnostics.diagnostic import *
from earthdiagnostics.modelingrealm import ModelingRealm, ModelingRealms
from earthdiagnostics.utils import Utils, TempFile
class InterpolateCDO(Diagnostic):
......@@ -201,6 +202,7 @@ class InterpolateCDO(Diagnostic):
must_exist=False, rename_dimension=True)
handler = Utils.openCdf(variable_file)
var = handler.variables[self.variable]
units = var.units
coordinates = list()
for dim in var.dimensions:
if dim == 'i':
......@@ -225,6 +227,11 @@ class InterpolateCDO(Diagnostic):
temp = TempFile.get()
Utils.cdo.remap(','.join((self.grid.split('_')[0], self.weights)), input=variable_file, output=temp)
handler = Utils.openCdf(temp)
handler.variables[self.variable].units = units
handler.close()
self.regridded.set_local_file(temp)
......
# coding=utf-8
import iris
import iris.coord_categorisation
import iris.coords
import iris.exceptions
import numpy as np
import six
from bscearth.utils.log import Log
......@@ -7,11 +12,6 @@ from earthdiagnostics.diagnostic import Diagnostic, DiagnosticVariableOption, Di
from earthdiagnostics.frequency import Frequencies
from earthdiagnostics.utils import TempFile
from earthdiagnostics.variable_type import VariableType
import numpy as np
import iris
import iris.coord_categorisation
import iris.exceptions
import iris.coords
class ClimatologicalPercentile(Diagnostic):
......@@ -107,24 +107,12 @@ class ClimatologicalPercentile(Diagnostic):
def _save_results(self, percentile_values):
temp = TempFile.get()
percentile_coord = iris.coords.DimCoord(ClimatologicalPercentile.Percentiles, long_name='percentile')
results = iris.cube.CubeList()
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.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
iris.save(results.merge_cube(), temp, zlib=True)
iris.save(percentile_values.merge_cube(), temp, zlib=True)
self.percentiles_file.set_local_file(temp, rename_var='percent')
def _calculate_percentiles(self):
Log.debug('Calculating percentiles')
percentiles = {}
bins = self.distribution.coord('bin').points
......@@ -135,10 +123,18 @@ class ClimatologicalPercentile(Diagnostic):
index = np.searchsorted(cs, percentile_values)
return [bins[i] for i in index]
results = iris.cube.CubeList()
percentile_coord = iris.coords.DimCoord(ClimatologicalPercentile.Percentiles, long_name='percentile')
print(self.distribution)
for leadtime_slice in self.distribution.slices_over('leadtime'):
leadtime = leadtime_slice.coord('leadtime').points[0]
percentiles[leadtime] = np.apply_along_axis(calculate, 0, leadtime_slice.data)
return percentiles
result = iris.cube.Cube(np.apply_along_axis(calculate, 0, leadtime_slice.data), var_name='percent',
units=self.distribution.coord('bin').units)
result.add_dim_coord(percentile_coord, 0)
result.add_dim_coord(leadtime_slice.coord('latitude'), 1)
result.add_dim_coord(leadtime_slice.coord('longitude'), 2)
result.add_aux_coord(leadtime_slice.coord('leadtime'))
results.append(result)
return results
def _get_distribution(self):
for startdate, startdate_file in six.iteritems(self.leadtime_files):
......@@ -148,3 +144,5 @@ class ClimatologicalPercentile(Diagnostic):
self.distribution = data_cube
else:
self.distribution += data_cube
if len(self.distribution.coords('leadtime')) == 0:
self.distribution.add_aux_coord(iris.coords.AuxCoord(1, var_name='leadtime', units='months'))
# coding=utf-8
import os
import iris
import iris.analysis
import iris.coord_categorisation
......@@ -80,17 +82,17 @@ class DaysOverPercentile(Diagnostic):
self.variable_file = self.request_chunk(self.domain, self.variable, self.startdate, None, None)
def declare_data_generated(self):
var_over = self.variable + '_daysover_q{0}'
var_below = self.variable + '_daysbelow_q{0}'
var_over = self.variable + '_daysover_q{0}_{1.start_year}-{1.end_year}'
var_below = self.variable + '_daysbelow_q{0}_{1.start_year}-{1.end_year}'
self.days_over_file = {}
self.days_below_file = {}
for perc in ClimatologicalPercentile.Percentiles:
self.days_over_file[perc] = self.declare_chunk(self.domain, var_over.format(int(perc * 100)),
self.days_over_file[perc] = self.declare_chunk(self.domain, var_over.format(int(perc * 100), self),
self.startdate, None,
None, frequency=Frequencies.monthly,
vartype=VariableType.STATISTIC)
self.days_below_file[perc] = self.declare_chunk(self.domain, var_below.format(int(perc * 100)),
self.days_below_file[perc] = self.declare_chunk(self.domain, var_below.format(int(perc * 100), self),
self.startdate, None,
None, frequency=Frequencies.monthly,
vartype=VariableType.STATISTIC)
......@@ -99,7 +101,6 @@ class DaysOverPercentile(Diagnostic):
"""
Runs the diagnostic
"""
raise Exception('Pues me enfado y no respiro!!!')
iris.FUTURE.netcdf_promote = True
percentiles = iris.load_cube(self.percentiles_file.local_file)
......@@ -140,8 +141,8 @@ class DaysOverPercentile(Diagnostic):
results_over = {perc: iris.cube.CubeList() for perc in ClimatologicalPercentile.Percentiles}
results_below = {perc: iris.cube.CubeList() for perc in ClimatologicalPercentile.Percentiles}
var_daysover = 'daysover'
var_days_below = 'daysbelow'
var_daysover = 'days_over'
var_days_below = 'days_below'
long_name_days_over = 'Proportion of days over a given percentile for {0.start_year}-{0.end_year} ' \
'climatology'.format(self)
long_name_days_below = 'Proportion of days below a given percentile for {0.start_year}-{0.end_year} ' \
......@@ -177,23 +178,56 @@ class DaysOverPercentile(Diagnostic):
Log.debug('Saving percentiles startdate {0}', self.startdate)
for perc in ClimatologicalPercentile.Percentiles:
iris.FUTURE.netcdf_no_unlimited = True
temp = TempFile.get()
iris.save(results_over[perc].merge_cube(), temp, zlib=True, unlimited_dimensions=['time'])
Utils.rename_variables(temp, {'dim2': 'ensemble', 'dim1': 'ensemble'},
must_exist=False, rename_dimension=True)
self.days_over_file[perc].set_local_file(temp, rename_var='daysover')
temp = TempFile.get()
iris.save(results_below[perc].merge_cube(), temp, zlib=True, unlimited_dimensions=['time'])
Utils.rename_variables(temp, {'dim2': 'ensemble', 'dim1': 'ensemble'},
must_exist=False, rename_dimension=True)
self.days_below_file[perc].set_local_file(temp, rename_var='daysbelow')
self.days_over_file[perc].set_local_file(self.save_to_file(perc, results_over, var_daysover),
rename_var=var_daysover)
self.days_below_file[perc].set_local_file(self.save_to_file(perc, results_below, var_days_below),
rename_var=var_days_below)
del self.days_over_file
del self.days_below_file
del self.lat_coord
del self.lon_coord
def save_to_file(self, perc, results_over, var_daysover):
temp = TempFile.get()
iris.save(results_over[perc].merge_cube(), temp, zlib=True, unlimited_dimensions=['time'])
Utils.rename_variables(temp, {'dim2': 'ensemble', 'dim1': 'ensemble'},
must_exist=False, rename_dimension=True)
handler = Utils.openCdf(temp)
if 'time' not in handler.dimensions:
new_file = TempFile.get()
new_handler = Utils.openCdf(new_file, 'w')
new_handler.createDimension('time', 1)
for dimension in handler.dimensions:
Utils.copy_dimension(handler, new_handler, dimension)
for variable in handler.variables.keys():
if variable in (var_daysover, 'time', 'time_bnds'):
continue
Utils.copy_variable(handler, new_handler, variable)
old_var = handler.variables[var_daysover]
new_var = new_handler.createVariable(var_daysover, old_var.dtype, ('time',) + old_var.dimensions,
zlib=True, fill_value=1.0e20)
Utils.copy_attributes(new_var, old_var)
new_var[0, :] = old_var[:]
old_var = handler.variables['time']
new_var = new_handler.createVariable('time', old_var.dtype, ('time',))
Utils.copy_attributes(new_var, old_var)
new_var[0] = old_var[0]
old_var = handler.variables['time_bnds']
new_var = new_handler.createVariable('time_bnds', old_var.dtype, ('time',) + old_var.dimensions)
Utils.copy_attributes(new_var, old_var)
new_var[0, :] = old_var[:]
new_handler.close()
os.remove(temp)
temp = new_file
handler.close()
return temp
def create_results_cube(self, days_over, percentile, realization_coord, time_coord,
var_name, long_name):
result = iris.cube.Cube(days_over.astype(np.float32), var_name=var_name, long_name=long_name, units=1.0)
......
# coding=utf-8
from bscearth.utils.date import parse_date, add_months
from bscearth.utils.log import Log
import math
from earthdiagnostics.diagnostic import Diagnostic, DiagnosticVariableOption, DiagnosticDomainOption, \
DiagnosticIntOption, DiagnosticFloatOption
from earthdiagnostics.utils import Utils, TempFile
from earthdiagnostics.variable_type import VariableType
import numpy as np
import iris
from iris.cube import Cube
import iris.coord_categorisation
from iris.time import PartialDateTime
import iris.exceptions
import iris.coords
import math
import iris.exceptions
import iris.unit
import numpy as np
import psutil
import six
from bscearth.utils.date import parse_date, add_months, add_days
from bscearth.utils.log import Log
from iris.cube import Cube
from iris.time import PartialDateTime
from earthdiagnostics.diagnostic import Diagnostic, DiagnosticVariableOption, DiagnosticDomainOption, \
DiagnosticIntOption, DiagnosticFloatOption
from earthdiagnostics.utils import Utils, TempFile
from earthdiagnostics.variable_type import VariableType
class Discretize(Diagnostic):
......@@ -184,14 +186,25 @@ class Discretize(Diagnostic):
cubes = iris.cube.CubeList()
date = parse_date(self.startdate)
date = add_days(date, 14, self.data_manager.config.experiment.calendar)
for leadtime, distribution in six.iteritems(self.distribution):
leadtime_cube = Cube(distribution.astype(np.uint32), var_name=self.data_cube.var_name,
standard_name=self.data_cube.standard_name, units='1')
leadtime_cube.add_dim_coord(bins_coord, 0)
leadtime_cube.add_dim_coord(self.data_cube.coord('latitude'), 1)
leadtime_cube.add_dim_coord(self.data_cube.coord('longitude'), 2)
leadtime_cube.add_aux_coord(iris.coords.AuxCoord(np.array((leadtime,), np.int8), var_name='leadtime',
leadtime_cube.add_aux_coord(iris.coords.AuxCoord(leadtime,
var_name='leadtime',
units='months'))
lead_date = add_months(date, leadtime - 1, self.data_manager.config.experiment.calendar)
leadtime_cube.add_aux_coord(iris.coords.AuxCoord(iris.unit.date2num(lead_date,
unit='days since 1950-01-01',
calendar="standard"),
var_name='time',
units='days since 1950-01-01'))
cubes.append(leadtime_cube)
temp = TempFile.get()
iris.FUTURE.netcdf_no_unlimited = True
......@@ -201,23 +214,20 @@ class Discretize(Diagnostic):
def _get_distribution(self):
self.distribution = {}
Log.debug('Discretizing...')
for leadtime in set(self.data_cube.coord('leadtime').points):
Log.debug('Discretizing leadtime {0}', leadtime)
leadtime_cube = self.data_cube.extract(iris.Constraint(leadtime=leadtime))
for realization in self.data_cube.coord('realization').points:
Log.debug('Discretizing realization {0}', realization)
if 'realization' in leadtime_cube.coords():
for realization_cube in self.data_cube.slices_over('realization'):
Log.debug('Discretizing realization {0}', realization_cube.coord('realization').points[0])
self.print_memory_used()
if leadtime not in self.distribution:
self.distribution[leadtime] = self._calculate_distribution(realization_cube)
else:
self.distribution[leadtime] += self._calculate_distribution(realization_cube)
else:
self.print_memory_used()
try:
realization_cube = leadtime_cube.extract(iris.Constraint(realization=realization))
except iris.exceptions.CoordinateNotFoundError:
realization_cube = leadtime_cube
if realization_cube is None and realization == 0:
realization_cube = leadtime_cube
if leadtime not in self.distribution:
self.distribution[leadtime] = self._calculate_distribution(realization_cube)
else:
self.distribution[leadtime] += self._calculate_distribution(realization_cube)
self.distribution[leadtime] = self._calculate_distribution(leadtime_cube)
# noinspection PyTypeChecker
def _get_value_interval(self):
......
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