Python Cdo.sub方法代码示例

# 需要导入模块: from cdo import Cdo [as 别名]
# 或者: from cdo.Cdo import sub [as 别名]
def normalize(resource, grouping='year', region='AUT', start_date="1971-01-01", end_date="2010-12-31", out_dir=None):
    """
    noramlize netcdf file for region and timeperiod

    :param resource: netcdf filename
    :param start_date: string with start date
    :param end_date: string with end date
    :param out_dir: output directory for result file (netcdf)

    :return: normalized netcdf file
    """
    rd = ocgis.RequestDataset(resource)
    variable = rd.variable
    from dateutil import parser as date_parser
    time_range=[ date_parser.parse(start_date) , date_parser.parse(end_date) ]
    calc = [{'func':'mean','name':'ref_' + variable }] 

    filename = drs_filename(resource)
    filename = filename.replace("EUR", region)
    from os.path import join
    output = join(out_dir, filename)
    prefix = "ref_%s" % filename
    prefix = prefix.replace('.nc', '')
    
    try: 
        reference = ocgis.OcgOperations(
            dataset=rd,
            geom=COUNTRY_SHP,
            dir_output=out_dir,
            output_format="nc",
            select_ugid=select_ugid(region),
            prefix=prefix,
            add_auxiliary_files=False,
            calc=calc,
            calc_grouping=calc_grouping(grouping),
            time_range=time_range  ).execute()

        from tempfile import mkstemp
        from cdo import Cdo   
        cdo = Cdo()
        _,out_resource = mkstemp(prefix="out_resource_", dir=out_dir)
        cdo.fldmean(input = resource , output = out_resource)
        _,out_ref = mkstemp(prefix="out_ref_", dir=out_dir)
        cdo.fldmean(input = reference , output = out_ref )
        cdo.sub(input = "%s %s" % (out_resource, out_ref) , output = output)
    except:
        msg = 'normalize failed for file : %s ' % filename
        logger.exception(msg)
        raise CalculationException(msg)
    return output    

# 需要导入模块: from cdo import Cdo [as 别名]
# 或者: from cdo.Cdo import sub [as 别名]
def get_anomalies(nc_file, frac=0.2, reference=None):
  '''
  anomalisation of data subsets for weather classification. 
  Anomalisation is done by substrcting a smoothed anual cycle 

  :parm nc_file: input netCDF file 
  :param frac: Number between 0-1 for stregth of smoothing 
               (0 = close to the original data, 1=flat line)
               default=0.2
  :param reference: Period to calulate anual cycle 
  :return string: path to output netCDF file
  '''
  try: 
    variable = utils.get_variable(nc_file)
    calc = [{'func': 'mean', 'name': variable}]
    calc_grouping = calc_grouping = ['day','year']
    nc_anual_cycle = call(nc_file, calc=calc, calc_grouping=calc_grouping, time_range=reference)
    logger.info('anual cycle calculated')  
  except Exception as e:
    msg = 'failed to calcualte anual cycle %s' % e
    logger.error(msg)
    raise Exception(msg)

  ### spline for smoothing
  import statsmodels.api as sm
  from numpy import tile, empty, linspace
  from netCDF4 import Dataset
  from cdo import Cdo
  cdo = Cdo()
  
  try:
    # variable = utils.get_variable(nc_file)
    ds = Dataset(nc_anual_cycle, mode='a')
    vals = ds.variables[variable]
    vals_sm = empty(vals.shape)
    ts = vals.shape[0]
    x = linspace(1, ts*3 , num=ts*3 , endpoint=True)

    for lat in range(vals.shape[1]):
      for lon in range(vals.shape[2]):
        try:
          y = tile(vals[:,lat,lon], 3)
          # ys = smooth(y, window_size=91, order=2, deriv=0, rate=1)[ts:ts*2]
          ys = sm.nonparametric.lowess(y, x, frac=frac )[ts:ts*2,1]
          vals_sm[:,lat,lon] = ys
        except Exception as e:
          msg = 'failed for lat %s lon %s  %s ' % (lat,lon,e)
          logger.debug('failed for lat %s lon %s  %s ' % (lat,lon,e))
          raise Exception(msg)
      print 'done for %s - %s ' % (lat, lon)    
    vals[:,:,:] = vals_sm[:,:,:]
    ds.close()
    logger.info('smothing of anual cycle done')
  except Exception as e:
    msg = 'failed smothing of anual cycle %s ' % e
    logger.error(msg)
    raise Exception(msg)
  try:   
    ip , nc_anomal = mkstemp(dir='.',suffix='.nc')
    nc_anomal = cdo.sub(input=[nc_file, nc_anual_cycle], output= nc_anomal )
    logger.info('anomalisation done: %s ' % nc_anomal)
  except Exception as e:
    msg = 'failed substraction of anual cycle %s ' % e
    logger.error(msg)
    raise Exception(msg)    
  return nc_anomal

# 需要导入模块: from cdo import Cdo [as 别名]
# 或者: from cdo.Cdo import sub [as 别名]
def modelUncertaintyWorker(resource):
  """retuns the result
  
  :param resource: list of pathes to netCDF files
  """
  from cdo import Cdo
  cdo = Cdo()
  from flyingpigeon.utils import check_timestepps

  # check resource for consistency
  resource_qc =  check_timestepps(resource)  
  
  try: 
    # ensemble mean and magnitude
    nc_ensmean = cdo.ensmean(input = resource_qc, output = 'nc_ensmean.nc')    
    logger.info('ensmean calculations done')
  except Exception as e: 
    logger.error('ensmean calculations failed: %s ' % e )

  try:     
    #nc_delta = nc_ensmean(lastpt) - nc_ensmean(firstpt)    
    nc_laststep = cdo.seltimestep('-1', input = nc_ensmean, output = 'nc_laststep.nc')
    nc_firststep = cdo.seltimestep(1, input = nc_ensmean, output = 'nc_firststep.nc')
    nc_delta = cdo.sub(input = [nc_laststep, nc_firststep], output = 'nc_delta.nc')
    nc_absdelta = cdo.abs(input = nc_delta, output = 'nc_absdelta.nc')
    logger.info('delta calculation done')
  except Exception as e: 
    logger.error('delta calculation failed: %s ' % e )  

  try: 
    # ensemble std    
    nc_ensstd = cdo.ensstd(input = resource_qc, output = 'nc_ensstd.nc')
    logger.info('std calculation done')
  except Exception as e: 
    logger.error('std calculation failed: %s ' % e )

  try:
    # compute mask: if nc_absdelta > nc_enssstd
    nc_level = cdo.mulc(1, input = nc_ensstd, output = 'nc_level.nc')
    nc_binmask = cdo.gt(input = [nc_absdelta, nc_level], output = 'nc_binmask.nc')
    logger.info('calculated mask')
  except Exception as e: 
    logger.error('mask calculation failed: %s ' % e )

  
  #if > const
  #cdo.gtc('arg', input=infile, ontput=outfile)

  # mean + sigma as a mask 
  # sigma1 = cdo.fldstd(input = nc_ensmean, output = 'nc_sigma1.nc')
  # mask 
  
  # mean + 2* sigma as a mask
  
  # merge to on result netCDF
  # cdo.merge(input=[file1, file2], output='result.nc')

  result = nc_ensmean    #ensemble mean
  result2 = nc_ensstd    #ensemble std
  result3 = nc_absdelta  #magnitude of model change  
  result4 = nc_binmask   #absdelta > std
  
  return result, result2, result3, result4

# 需要导入模块: from cdo import Cdo [as 别名]
# 或者: from cdo.Cdo import sub [as 别名]

#.........这里部分代码省略.........
    if timeslice == None: 
      timeslice = int((end - start) / 3)
      if timeslice == 0: 
        timeslice = 1
    else: 
      timeslice = int(timeslice)
    start1 = start
    start2 = start1 + timeslice - 1 
    end1 = end - timeslice + 1
    end2 = end
    logger.info('timeslice and periodes set')
  except Exception as e:
    msg = 'failed to set the periodes'
    logger.exception(msg)
    raise Exception(msg)

  try:
    files = []
    for i, mf in enumerate(mergefiles):
      files.append(cdo.selyear('{0}/{1}'.format(start1,end2), input=[mf.replace(' ','\ ')] , output='file_{0}_.nc'.format(i) )) #python version
    logger.info('timeseries selected from defined start to end year')
  except Exception as e:
    msg = 'seltime and mergetime failed'
    logger.exception(msg)
    raise Exception(msg)    

  try: 
    # ensemble mean 
    nc_ensmean = cdo.ensmean(input=files , output='nc_ensmean.nc')
    logger.info('ensemble mean calculation done')
  except Exception as e:
    msg = 'ensemble mean failed'
    logger.exception(msg)
    raise Exception(msg)
  
  try: 
    # ensemble std 
    nc_ensstd  = cdo.ensstd(input=files , output='nc_ensstd.nc')
    logger.info('ensemble std and calculation done')
  except Exception as e:
    msg = 'ensemble std or failed'
    logger.exception(msg)
    raise Exception(msg)
  
  # get the get the signal as difference from the beginning (first years) and end period (last years), :
  try:
    selyearstart = cdo.selyear('%s/%s' % (start1,start2), input = nc_ensmean, output = 'selyearstart.nc' ) 
    selyearend = cdo.selyear('%s/%s' % (end1,end2), input = nc_ensmean, output = 'selyearend.nc' )
    meanyearst = cdo.timmean(input = selyearstart, output= 'meanyearst.nc')
    meanyearend = cdo.timmean(input = selyearend, output= 'meanyearend.nc')
    signal = cdo.sub(input=[meanyearend, meanyearst], output = 'signal.nc')
    logger.info('Signal calculation done')
  except Exception as e:
    msg = 'calculation of signal failed'
    logger.exception(msg)
    raise Exception(msg)
  
  # get the intermodel standard deviation (mean over whole period)
  try:
    #std_selyear = cdo.selyear('%s/%s' % (end1,end2), input=nc_ensstd, output='std_selyear.nc')
    #std = cdo.timmean(input = std_selyear, output = 'std.nc')
    
    std = cdo.timmean(input = nc_ensstd, output = 'std.nc')
    std2 = cdo.mulc('2', input = std, output = 'std2.nc')
    logger.info('calculation of internal model std for time period done')
  except Exception as e:
    msg = 'calculation of internal model std failed'
    logger.exception(msg) 
    raise Exception(msg)
  try:
    absolut = cdo.abs(input=signal, output='absolut_signal.nc')
    high_agreement_mask = cdo.gt(input=[absolut,std2],  output= 'large_change_with_high_model_agreement.nc')
    low_agreement_mask = cdo.lt(input=[absolut,std], output= 'small_signal_or_low_agreement_of_models.nc')
    logger.info('high and low mask done')
  except Exception as e:
    msg = 'calculation of robustness mask failed'
    logger.exception(msg)
    raise Exception(msg)
  
  try: 
    if variable == None: 
      variable = get_variable(signal)
    logger.info('variable to be plotted: %s' % variable)
    
    if title == None: 
      title='Change of %s (difference of mean %s-%s to %s-%s)' % (variable, end1, end2, start1, start2)  
    
    graphic = None
    graphic = map_ensembleRobustness(signal, high_agreement_mask, low_agreement_mask, 
              variable=variable, 
              cmap=cmap,
              title = title)
    
    logger.info('graphic generated')
  except Exception as e:
    msg('graphic generation failed: %s' % e)
    logger.debug(msg)
    raise Exception(msg)

  return signal, low_agreement_mask, high_agreement_mask, graphic, text_src #