import netCDF4 import os from earthdiagnostics import Utils, cdftools, cdo, TempFile from earthdiagnostics.basins import Basins from earthdiagnostics.diagnostic import Diagnostic from autosubmit.config.log import Log import shutil import threading class Siasiesiv(Diagnostic): """ Class containing all the diagnostics related to ice """ def __init__(self, data_manager, basin, startdate, member, chunk, lock): Diagnostic.__init__(self, data_manager) self.basin = basin if basin == Basins.Global: self._region = None else: self._region = basin.shortname self.startdate = startdate self.member = member self.chunk = chunk self.required_vars = ['sit', 'sic'] self.generated_vars = ['siextents', 'sivols', 'siareas', 'siextentn', 'sivoln', 'siarean'] self.lock = lock @staticmethod def generate_jobs(data_manager, startdates, members, chunks, options): basin = Basins.parse(options[1]) lock = threading.Lock() job_list = list() for startdate in startdates: for member in members: for chunk in range(1, chunks + 1): job_list.append(Siasiesiv(data_manager, basin, startdate, member, chunk, lock)) return job_list def compute(self): """ Compute the sea ice extent (1000km2), area (1000km2), volume (km3) and mean thickness (m) in both hemispheres or a specified region. Created in April 2012 Author : vguemas@ic3.cat Modified in June 2014 Author : neven.fuckar@ic3.cat Computation of the properties in various selected regions according to mask.regions.${NEMOVERSION}.nc (mask_regions.nc) is based on modification of mask.regions.ORCA1.noverticalinfo.Matt.nc from Matthieu Chevallier. :return: """ nco = Utils.nco sit_file = self.data_manager.get_file('seaIce', 'sit', self.startdate, self.member, self.chunk) sic_file = self.data_manager.get_file('seaIce', 'sic', self.startdate, self.member, self.chunk) nco.ncks(input=sic_file, output=sit_file, options='-A -v sic') ntime = int(cdo.ntime(input=sit_file)[0]) files = list() self.lock.acquire() if self.basin != Basins.Global: shutil.move('mask.nc', 'original_mask.nc') shutil.move('mask_regions.nc', 'mask.nc') Utils.rename_variable('mask.nc', self.basin.fullname, 'tmask') error = None try: for time in range(ntime): Log.info('Running time {0}', time) temp = TempFile.get() nco.ncks(input=sit_file, output=temp, options='-O -d time,{0}'.format(time)) cdftools.run('cdficediags', input=temp) Utils.move_file('icediags.nc', temp) files.append(temp) except Exception as ex: error = ex.message finally: if self.basin != Basins.Global: Utils.rename_variable('mask.nc', 'tmask', self.basin.fullname) shutil.move('mask.nc', 'mask_regions.nc') shutil.move('original_mask.nc', 'mask.nc') if error: raise Exception(error) self.lock.release() temp = TempFile.get() cdo.cat(input=' '.join(files), output=temp) nco.ncks(input=sit_file, output=temp, options='-A -v time') self.data_manager.send_file(self._extract_variable_and_rename(temp, 'SVolume', 'sivols', "10^3 km3"), 'seaIce', 'sivols', self.startdate, self.member, self.chunk, region=self.basin.fullname) self.data_manager.send_file(self._extract_variable_and_rename(temp, 'SArea', 'siareas', "10^6 km2"), 'seaIce', 'siareas', self.startdate, self.member, self.chunk, region=self.basin.fullname) self.data_manager.send_file(self._extract_variable_and_rename(temp, 'SExnsidc', 'siextents', "10^6 km2"), 'seaIce', 'siextents', self.startdate, self.member, self.chunk, region=self.basin.fullname) self.data_manager.send_file(self._extract_variable_and_rename(temp, 'NVolume', 'sivoln', "10^3 km3"), 'seaIce', 'sivoln', self.startdate, self.member, self.chunk, region=self.basin.fullname) self.data_manager.send_file(self._extract_variable_and_rename(temp, 'NArea', 'siarean', "10^6 km2"), 'seaIce', 'siarean', self.startdate, self.member, self.chunk, region=self.basin.fullname) self.data_manager.send_file(self._extract_variable_and_rename(temp, 'NExnsidc', 'siextentn', "10^6 km2"), 'seaIce', 'siextentn', self.startdate, self.member, self.chunk, region=self.basin.fullname) os.remove(temp) def _extract_variable_and_rename(self, input_file, variable, cmor_name, output_units): temp = TempFile.get() # Utils.nco.ncks(input=input_file, output=temp, options='-O -v {0},time,time_bnds'.format(variable)) input_handler = Utils.cdo.openCdf(input_file) os.remove(temp) handler = netCDF4.Dataset(temp, 'w') # Create dimensions handler.createDimension('time', None) handler.createDimension('bnds', 2) # Copy time variable time_var = input_handler.variables['time'] new_time = handler.createVariable('time', time_var.datatype, time_var.dimensions) new_time.setncatts({k: time_var.getncattr(k) for k in time_var.ncattrs()}) new_time[:] = time_var[:] original_bnds = input_handler.variables['time_bnds'] new_bnds = handler.createVariable('time_bnds', original_bnds.datatype, original_bnds.dimensions) new_bnds.setncatts({k: original_bnds.getncattr(k) for k in original_bnds.ncattrs()}) new_bnds[:] = original_bnds[:] original_variable = input_handler.variables[variable] values = original_variable[:, 0, 0] new_var = handler.createVariable(cmor_name, original_variable.datatype, 'time') new_var.setncatts({k: original_variable.getncattr(k) for k in original_variable.ncattrs()}) factor = self._get_conversion_factor(original_variable.units, output_units) values *= factor new_var[:] = values new_var.units = output_units new_var.short_name = cmor_name handler.close() return temp def _get_conversion_factor(self, input_units, output_units): units = input_units.split() if len(units) == 1: scale_unit = 1 unit = units[0] else: if '^' in units[0]: values = units[0].split('^') scale_unit = pow(int(values[0]), int(values[1])) else: scale_unit = float(units[0]) unit = units[1] units = output_units.split() if len(units) == 1: scale_new_unit = 1 new_unit = units[0] else: if '^' in units[0]: values = units[0].split('^') scale_new_unit = pow(int(values[0]), int(values[1])) else: scale_new_unit = float(units[0]) new_unit = units[1] factor = self._get_factor(new_unit, unit) invert = False if factor is None: factor = self._get_factor(unit, new_unit) invert = True if factor is None: raise Exception("Conversion from {0} to {1} not supported".format(units, output_units)) if invert: factor = scale_unit / float(scale_new_unit * factor) else: factor = (factor * scale_unit) / float(scale_new_unit) return factor def _get_factor(self, new_unit, unit): # Add only the conversions with a factor greater than 1 if unit == new_unit: return 1 if unit == 'km3': if new_unit == 'm3': return pow(1000, 3) elif unit == 'km2': if new_unit == 'm2': return pow(1000, 2) return None