source: tags/ORCHIDEE_OL/Utilitaire/Regrid_restart/ajoute_bords.py @ 6

#!/bin/env cdat
# -*- coding: iso-8859-15 -*-

# python system
import os, sys, gc, getopt, time
#import string
import re
import Numeric


# CDAT/CDMS
import cdms, MV
#, cdutil
from cdms.axis import TransientAxis
from cdms.coord import TransientVirtualAxis
from dataproc.utils import cdat

def usage():
    """Display usage information."""
    
    man = """
    This script with corrected the mask of a restart.
    It will get file_to_modify datas and the mask of file_to_add file
    and make a new file_out. The missing datas in file_to_modify will be simply
    replaced by the same points in file_to_add file.
    
Usage : 
    prompt> cdat ajoute_bords.py file_to_add file_to_modify file_out [variable]
    Add some points defines in file_to_add for all variables (or just one) of file_to_modify to file_out.
    
    examples :
    > ajoute_bords.py CM1414A_20001231_sechiba_rest.nc regrid_CF_CM5PIRC11_29901231_sechiba_rest.nc grille_regrid_CF_CM5PIRC11_29901231_sechiba_rest.nc
    > ajoute_bords.py CM1414A_20001231_stomate_rest.nc regrid_CF_CM5PIRC11_29901231_stomate_rest.nc grille_regrid_CF_CM5PIRC11_29901231_stomate_rest.nc
    > ajoute_bords.py CM5PIRC8_26601231_stomate_rest.nc regrid_CF_OOL-CTRLPI-1_9_5_SPIN_v2ORC_43_01001231_stomate_rest.nc grille_regrid_CF_OOL-CTRLPI-1_9_5_SPIN_v2ORC_43_01001231_stomate_rest.nc

"""
    print man


options = [
    "help",
    "version" ]



# Extract arguments from the command line
try:
    myopts, myargs = getopt.getopt(sys.argv[1:], 'Vh', longopts = options)
except getopt.GetoptError:
    print 'getopt error'
    usage()
    sys.exit(1)

if (len(myargs) < 3):
    print 'I need at least three files.'
    usage()
    sys.exit(1)
else:
    FileIn = myargs[0]
    print "file_to_add : ",FileIn
    if (os.path.exists(FileIn)==False):
        print 'file ',FileIn,' doesn''t exist !' 
        sys.exit(1)
    else:
        fin=cdms.open(FileIn, 'r')

    FileOut = myargs[1]
    print "file_to_modify : ",FileOut
    if (os.path.exists(FileOut)==False):
        print 'file ',FileOut,' doesn''t exist !' 
        sys.exit(1)
    else:
        fout=cdms.open(FileOut, 'r')

    FileNew = myargs[2]
    print "file_out : ",FileNew
    fnew=cdms.open(FileNew, 'w')

    if (len(myargs) > 3):
        varlist = [myargs[3]]
        print 'Only one variable :',varlist
    else:
        varlist = fin.listvariables()
        #varlist = [varlist[0]]
        print "Liste des variables de "+FileIn,varlist
        
# Affichage des variables : fixe la limite max
MV.set_print_limit(1700)

var0=fin(varlist[0])
missing = var0.missing_value
missing_ = var0.missing_value - 1.
print missing, missing_

for var in varlist:
    print "-----"
    print "Variable : ",var,"------"

    if (var in ["nav_lon","nav_lat","nav_lev","time_steps","time",
                "lat","lon","lev","day_counter",
                "dt_days", "routingcounter", "veget_year", "date" ]):
        var1=fin(var)
        fnew.write(var1,attributes=cdat.getAttributes(var1), axes=var1.getAxisList(), id=var1.id)
        print "remplace."
        continue
    if (var in ["bounds_lat","bounds_lon","routingcounter","routingarea","routetogrid","routetobasin",
                "basinid","topoindex","fastres","slowres","streamres","lakeres","previous_outflow",
                "lakeinflow","returnflow","riverflow","coastalflow","hydrographs","runoff_route",
                "drainage_route","evapot_route","precip_route","humrel_route","totnobio_route","vegtot_route"]):
        print "supprime."
        continue

    var0 = fin(var)
    var0m = MV.getmaskarray(var0)

    #var0_notmissing=(var0 < missing_)
    var0_notmissing=(var0m == 0)

    var1 = fout(var)

    var1m = MV.getmaskarray(var1)
    #var1_missing=(var1 > missing_)
    var1_missing=(var1m == 1)
#    var1or = MV.logical_or(var1,var1m)

    #test = MV.logical_and(var0_notmissing,var1or)
    test = MV.logical_and(var0_notmissing,var1_missing)
#     if ( var in ["temp_sol"] ):
#         fnew.write(var0m,attributes=cdat.getAttributes(var1), axes=var1.getAxisList(), id="var0m")
#         fnew.write(var0_notmissing,attributes=cdat.getAttributes(var1), axes=var1.getAxisList(), id="var0_notmissing")
#         fnew.write(var1_missing,attributes=cdat.getAttributes(var1), axes=var1.getAxisList(), id="var1_missing")
#         fnew.write(var1m,attributes=cdat.getAttributes(var1), axes=var1.getAxisList(), id="var1m")
#         fnew.write(var1or,attributes=cdat.getAttributes(var1), axes=var1.getAxisList(), id="var1or")
#         fnew.write(test,attributes=cdat.getAttributes(var1), axes=var1.getAxisList(), id="test")

    var1b = MV.where( test , var0, var1)
    
#     var1w = cdms.createVariable(var1b,
#                                fill_value = missing,
#                                attributes = cdat.getAttributes(var1),
#                                id = var,
#                                axes = var1.getAxisList())
#     print var1w

#     var1w = cdms.createVariable(var1b,
#                                 attributes = cdat.getAttributes(var1),
#                                 id = var,
#                                 axes = [var1.getAxis(0), var1.getAxis(1), var1.getAxis(2)])
#     print var1w,cdat.getAttributes(var1w),

    fnew.write(var1b,attributes=cdat.getAttributes(var1), axes=var1.getAxisList(), id=var1.id)

fnew.close()