#!/usr/bin/env python
#
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
#    http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing,
# software distributed under the License is distributed on an
# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied.  See the License for the
# specific language governing permissions and limitations
# under the License.

import cdms2
import cdtime
import cmor
import sys
import getopt
import factory
import numpy
from factory.formats import import_equation
from Toolbox.ESGFresources import *
from Toolbox.ESGFexcel import *
from Toolbox.CMORresources import CMORTable


# ************************************************************************
#                              process()                               
#                                                                      
#      Convert to obs4MIPS file format.                                
# ************************************************************************
def process( rc ):
    '''
    Convert netcdf/matlab/grads files into CMIP5 format.
    '''
    pdb.set_trace()
    # ----------------------------
    #  Loop yearly on file list.  
    # ----------------------------
    file_template = rc['file_template'].split(",");
    if( len(file_template) == 2 ):
        template_parameter = file_template[1]
        rc['file_template'] = file_template[0]
    else:
        template_parameter = 'years'

    for year in rc[template_parameter].split(","):
        if(year == ''):
            files= os.popen( "ls " + rc['file_template'] ).readlines()
        else:
            # ------------------------------------------------
            # Use string formating for path with same argument 
            # ------------------------------------------------
            try:
                tmplFile = rc['file_template'] % (year)
            except:
                tmplFile = rc['file_template'].format(year)
            if( not os.path.isfile( tmplFile) ) :
                print "**** Warning %s not found\n" % ( tmplFile )
                continue
            files= os.popen( "ls " + tmplFile).readlines()

        if( files == [] ):
            print "No file found: Check your resource file"
            return -1
        # ------------------------------------------------
        # Get the right handler to manage this file format
        # ------------------------------------------------
        Handler = factory.HandlerFormats(files[0].strip())

        # -----------------------------------
        # Take care of cmor initialization.
        # -----------------------------------
        cmor.setup(inpath=rc['inpath'],
                   netcdf_file_action = cmor.CMOR_REPLACE)
      
        cmor.dataset(experiment_id = rc['experiment_id'], 
                     institution   = rc['institution'  ],
                     calendar      = rc['calendar'     ],
                     institute_id  = rc['institute_id' ],
                     model_id      = rc['model_id'     ],
                     source        = rc['source'       ],
                     contact       = rc['contact'      ],
                     references    = rc['references'   ])
      
	
        # -----------------------------------------
        # add extra Global Attributes for obs4MIPs.
        # -----------------------------------------
        cmor.set_cur_dataset_attribute( 'instrument',     rc['instrument'    ])
        cmor.set_cur_dataset_attribute( 'mip_specs',      rc['mip_specs'     ])
        cmor.set_cur_dataset_attribute( 'data_structure', rc['data_structure'])
        cmor.set_cur_dataset_attribute( 'source_type',    rc['source_type'   ])
        cmor.set_cur_dataset_attribute( 'source_id',      rc['source_id'     ])
        cmor.set_cur_dataset_attribute( 'realm',          rc['realm'         ])
        cmor.set_cur_dataset_attribute( 'obs_project',    rc['obs_project'   ])
        cmor.set_cur_dataset_attribute( 'processing_version',
                                        rc['processing_version'] )
        cmor.set_cur_dataset_attribute( 'processing_level',
                                        rc['processing_level'] )

        cmor.load_table(rc['table'])

        # ---------------------------------------------------------------------
        # We loop on each file found, a new cmor file will be create on each
        # iteration.  If you want to aggregate, you need to use Grads ctl file
        # or NeCDF list of file.
        # ---------------------------------------------------------------------
        for file in files:
	
            fnm=file.strip()    # Get rid of \n
            aVariable  = eval(rc['original_var'])
            nbVariable = len(aVariable)

            # -----------------------------------------------------
            # ECMWF needs synoptic time 00z and 12z in he filename.
            # We extract it from the first file in the list.
            # -----------------------------------------------------
            if( rc['source_fn'] == 'SYNOPTIC' ):
                index = fnm.find('z.')
                rc['SYNOPTIC'] = fnm[index-2:index]

            # -----------------------
            # Work on all variables
            # -------------------------
            for j in arange(nbVariable):
                # ----------------------------------------------------
                # Fetch the variable directly or excute equation.
                # ----------------------------------------------------
                try:
                    variable=aVariable[j]
                    Handler.open(fnm, variable=variable)
                    rc['cvrt_original_var']   = aVariable[j]
                    print "Working on variable %s " % variable
                except:
                    if( aVariable[j] != 'equation' ) :
                        print "Variable %s can't open" % variable
                        continue
                    else:
                        print "Executing %s " % eval(rc['equation'])[j]
                
#                pdb.set_trace()
                rc['cvrt_original_units'] = eval(rc['original_units'])[j]
                rc['cvrt_cmor_var']       = eval(rc['cmor_var'])[j]
                rc['cvrt_equation']       = eval(rc['equation'])[j]
                rc['cvrt_level']          = eval(rc['level'])[j]
                
                data=Handler.getData()

                # ----------------------------------------------------------
                # Evaluate equation if needed. Usually used to change units
                # ----------------------------------------------------------
                if( rc['cvrt_equation'][0] == '@' ):
                    fncName = rc['cvrt_equation'][1:]
                    fnc = import_equation( "equations.%s" %  fncName )
                    data[:]= fnc(Handler)
                else:
                    data[:]=eval(rc['cvrt_equation'])
         
                # -------------------------------------------------------------
                # Save filled value in case data type is changed in createAxes
                # -------------------------------------------------------------
                fill_value = data.fill_value

                #  ---------------------------------------------
                #  Extract latitude/longitude
                #  ---------------------------------------------
                lonvals=Handler.getLongitude()
                latvals=Handler.getLatitude()
                # ---------------------
                # Create cmor time axis
                # ----------------------
                (rel_time, rel_time_bnds) = createTime(Handler, rc)

                # ---------------------------------------------------
                # Create cmor axes and add an axis to data if needed
                # ---------------------------------------------------
                (axes, data) = createAxes( rc, latvals, lonvals, data )
                                
                axis_ids = list()
                for axis in axes:
                    axis_id = cmor.axis(**axis)
                    axis_ids.append(axis_id)
            
                # ----------------------------------------------------------
                # Create cmor variable
                # Note: since this is in the loop, a new cmor file will be
                # create for each cmor write command.
                # ----------------------------------------------------------
                varid = cmor.variable(table_entry   = rc['cvrt_cmor_var'],
                                      axis_ids      = axis_ids,
                                      history       = '',
                                      missing_value = fill_value,
                                      original_name = rc['cvrt_original_var'],
                                      units         = rc['cvrt_original_units']
                                      )
                # -------------------------------
                # Write data for this time frame.
                # -------------------------------
                cmor.write(varid,data,\
                           time_vals=rel_time,time_bnds=rel_time_bnds)
                cmor.close(varid,file_name=True)

                # ---------------------------------------
                # Rename cmor files according to project.
                # ---------------------------------------
                if( movefiles(rc) ):
                    return -2 

        cmor.close()
    return 0


# ********************************************************************
#
#      createTime()
#
#   Define Time and Time bound axes for cmor
#
# ******************************************************************** 
def createTime(Handler, rc):
    '''
    InputtimeUnits: specified from resource file or from first file
    in a list of file.
    
    return relative time and time bounds using OutputTimeUnits from
    resource file.
    '''
    # ----------------------------------------------------
    # Retrieve time units from file if not provided in the
    # resource file.
    # ----------------------------------------------------
    InputTimeUnits = Handler.getTimeUnits(rc['InputTimeUnits'])
    
    #  --------------------------------------------------------
    #  Create time relative to January 1st 1900 to facilitate
    #  Threds software file handling.
    #  -------------------------------------------------------

    cur_time = Handler.getTime(InputTimeUnits)
    
    rel_time     =[cur_time[i].torel(rc['OutputTimeUnits']).value  
                   for i in range(len(cur_time))]
    
    if( len(rel_time) == 1 ) :
	deltarel = 1
    else:
       deltarel = rel_time[2] - rel_time[1]

    rel_time_bnds = rel_time[:]
    rel_time_bnds.append(rel_time[-1]+deltarel)
    return rel_time, rel_time_bnds

# ********************************************************************
# 
#  getCMIP5lev()
#
#  Extract CMIP5 mandatory level and recreate a new data array.
#  They are 16 mandatory levels.
# 
# ********************************************************************
def getCMIP5lev(data,rc):
    '''
    '''
    oTable               = CMORTable(rc['inpath'], rc['table'], "plevs")
    # ----------------------
    # Extract spefied levels
    # ----------------------
    if( 'levels' in oTable.dico.keys() ):
        #pdb.set_trace() 
        dataLevels = data.getLevel()[:]
        if( data.getLevel().units == "millibars" or
            data.getLevel().units == "hPa"       or
            data.getLevel().units == "mbar"    ):
            # --------------------------
            # Change units for to Pascal
            # ---------------------------
            LevelScaleFactor = 100
            dataLevels = data.getLevel()[:] * LevelScaleFactor

        # ----------------------------------------
        # No level selected, return all data array
        # ----------------------------------------
        if( len(rc['cvrt_level'].split(":")) == 1 ):
            levels =  [ float(item) for item in dataLevels ]
            lev=cdms2.createAxis( levels )
            lev.designateLevel()
            lev.units="pa"
            lev.long_name=data.getLevel().long_name
            #lev.id="lev"
            #lev=data.getAxis(1)
            #lev.__setattr__('_data_',dataLevels.astype(float))
            #lev.__setattr__('units',"Pa")
            #lev.units="hPa"
            data2=data.pressureRegrid(lev)
            return data2
        
        if( rc['cvrt_level'].split(':')[1] == "CMIP5" ):
            lev=cdms2.createAxis( [ float(item/LevelScaleFactor)
                                    for item in dataLevels
                                    if item in oTable.dico['levels' ] ] )

            lev.designateLevel()
            lev.units="pa"
            lev.long_name = data.getLevel().long_name
            data2=data.pressureRegrid(lev)
            lev[:]=lev[:]*LevelScaleFactor
            return data2
        else:
            # -----------------------
            # Assume a list of values
            # -----------------------
            levels = rc['cvrt_level'].split(':')[1].split(",")
            # --------------------------
            # Change units to Pascal
            # ---------------------------
            dataLevels = [ float(rc['cvrt_level'].split(":")[1].split(",")[i]) * \
                           LevelScaleFactor for i in range(len(levels)) ]
            # -----------------------------------
            # Match dataLevels with CMIP5 levels
            # Use file units
            # -----------------------------------
            lev=cdms2.createAxis( [ float(item/LevelScaleFactor)
                                    for item in dataLevels
                                    if item in oTable.dico['levels' ] ] )
            # -----------------------------------
            # Set axis metadata
            # -----------------------------------
            lev.units="pa"
            lev.long_name = data.getLevel().long_name
            lev.designateLevel()
            # -----------------------------------
            # Extract specified levels
            # -----------------------------------
            data2=data.pressureRegrid(lev)
            # -----------------------------------
            # Scale data back
            # -----------------------------------
            lev[:]=lev[:]*LevelScaleFactor
            return data2
            

    return data
# ********************************************************************
#
#      createAxes()
#
#   Define axes required by cmor and add z axis to data if needed
#
# ******************************************************************** 
def createAxes(rc, latvals, lonvals, data):
    #  ---------------------------------------------
    #  Create time/lat/lon axes using a dictionary
    #  ---------------------------------------------
    axes = [ 
        {'table_entry' : 'time',
         'units'       : rc['OutputTimeUnits']},
        
        {'table_entry' : 'latitude',
         'units'       : 'degrees_north',
         'coord_vals'  : latvals,
         'cell_bounds' : latvals.getBounds()},             
        
        {'table_entry' : 'longitude',
         'units'       : 'degrees_east',
         'coord_vals'  : lonvals,
         'cell_bounds' : lonvals.getBounds()},
        ]
    
    fill_value = data.fill_value

    if( rc['cvrt_level'] == 'height2m' ):
        axes.append({'table_entry' : 'height2m',
                     'units'       : 'm',
                     'coord_vals'  : [2.0] })
        data = numpy.array(data[:])
        data = data[:,:,:,numpy.newaxis]

    elif( rc['cvrt_level'] != '' ):
        data = getCMIP5lev( data, rc )
        levels=data.getLevel()[:]
        axes = numpy.insert(axes, 1,
                           {'table_entry' : 'plevs',
                            'units'       : 'Pa',
                            'coord_vals'  : levels })


    return axes, data

# ********************************************************************
#
#      usage()                                                          
#                                                                     
# ******************************************************************** 
def usage(message):
    '''
    Describe program synopsis.
    '''
    print
    print "*************************"
    print message
    print "*************************"
    print
    print
    print "obs4MIPS_process.py [-h] -r resource"
    print "   resource:   File containing Global attributes"
    print ""
    print "obs4MIPS will convert an input data file into CMIP5 format using "
    print "CMOR.  A directory path will be creating using CMOR by default or "
    print "using a template provided in the resource file."
    print
   
# ********************************************************************
#
#      main()                                                          
#                                                                     
# ******************************************************************** 
def main():
    '''
    '''
    pdb.set_trace()
    try:
        opts, args = getopt.getopt(sys.argv[1:], "hy:r:x:",
                                   ["help" ,"year=","resource=","excel="])
    except getopt.GetoptError, err:
        usage(str(err))# will print something like "option -a not recognized"
        return(2)

    # --------------------------
    # Verify passed arguments
    # --------------------------
    year     = -1
    resource = None
    excel    = None
    for o, a in opts:
        if o in ("-r", "--resource"):
            resource = a
        elif o in ("-x", "--excel"):
            excel = a
        elif o in ("-h", "--help"):
            usage()
            return(0)
        elif o in ("-y", "--year"):
            yr = a
        else:
            assert False, "unhandled option"

    # ------------------------------
    # Does the resource file exist?
    # ------------------------------
    if( ((resource == None ) or ( not os.path.isfile( resource ) )) and (( excel == None ) or ( not os.path.isfile( excel ) )) ):
        usage("bad Input Resource/Excel File")
        return 1

    # -----------------------
    # Read in "rc" file
    # -----------------------
    if( resource ):
       rc = ESGFresources( resource )
    if( excel ):
       rc = ESGFexcel( excel )

    # --------------------------------
    # Extract CMIP5 Table information
    # --------------------------------
    oTable               = CMORTable(rc['inpath'], rc['table'])
    if( not 'original_var' in rc.resources.keys() ):
        sys.exit(-1)
    rc['project_id']     = oTable[ 'project_id'     ]
    rc['product']        = oTable[ 'product'        ]
    rc['modeling_realm'] = oTable[ 'modeling_realm' ]
    rc['frequency']      = oTable[ 'frequency'      ]
    if( process(rc) ):
        return -1
            
    return 0

# ********************************************************************
#
#      Call main program and return exit code
#                                                                     
# ******************************************************************** 
if __name__ == '__main__':
    sys.exit(main())