source: nappy/tags/nappy-beta-0-1/cdms2na.py @ 349

Subversion URL: http://proj.badc.rl.ac.uk/svn/ndg/nappy/tags/nappy-beta-0-1/cdms2na.py@3343
Revision 349, 20.8 KB checked in by selatham, 16 years ago (diff)

Inserted license information.

  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
Line 
1#!/usr/bin/env python
2#    Copyright (C) 2004 CCLRC & NERC( Natural Environment Research Council ).
3#    This software may be distributed under the terms of the
4#    Q Public License, version 1.0 or later. http://ndg.nerc.ac.uk/public_docs/QPublic_license.txt
5
6helpMessage="""
7
8cdms2na.py
9==========
10
11Converts a CdmsFile object into one a NASA Ames file.
12
13Usage
14=====
15
16    cdms2na.py [<options>] -i <infilename> -o <outfilename>
17
18Where
19-----
20
21    infilename  - name of input file (NetCDF).
22    outfilename - name of output file (NASA Ames).
23    options     - list of options [NOT YET IMPLEMENTED].
24   
25"""
26
27# Imports from python standard library
28import sys, os, time, string, fileinput, re
29sys.path.append("..")
30
31# Import cdat modules
32if sys.platform.find("win")>-1:
33    pass
34else:
35    import cdms, cdtime
36    # Define cdms variables
37    cdms.setAutoBounds("off")
38
39# Import nappy modules
40import cdmsMap
41import version
42import general
43import localRules
44from naCore import *
45from naError import *
46
47def compareAxes(ax1, ax2):
48    """Takes 2 cmds axis objects returning 1 if they are essentially
49    the same and 0 if not."""
50    for axtype in ("time", "level", "latitude", "longitude"):
51        if cdms.axisMatches(ax1, axtype) and not cdms.axisMatches(ax2, axtype):
52            return 0
53
54    # Check ids
55    if ax1.id!=ax2.id: return 0
56    # Check lengths
57    if len(ax1)!=len(ax2): return 0
58    # Check values
59    if ax1._data_!=ax2._data_: return 0
60    # Check units
61    if ax1.units!=ax2.units: return 0
62    # OK, I think they are the same axis!
63    return 1
64
65def compareVariables(var1, var2):
66    """Compares two cdms variables to see if they are defined on identical
67    axes. Returns 1 if true, 0 if not."""
68    try:
69        for i in range(len(var1.getAxisList())):
70            ax1=var1.getAxis(i)
71            ax2=var2.getAxis(i)       
72            if compareAxes(ax1, ax2)==0:
73                return 0
74    except:
75        return 0
76    return 1
77
78def isAuxAndVar(avar, var):
79    """Compares two cdms variables and returns true if the first dimension of the
80    main variable is the only dimension the auxiliary variable is defined against."""
81    if len(avar.shape)>1:
82        return 0
83    auxax=avar.getAxis(0)
84    varax=var.getAxis(0)
85    return compareAxes(auxax, varax)
86
87def arrayToList(array, inlist):
88
89    """ Takes an n-dimensional Numeric array and converts it to an
90    n-dimensional list object.
91    """
92    dimlist=array.shape
93    if len(dimlist[1:])>0:
94            for i in range(dimlist[0]):
95                arrayToList(inlist[i], array[i])
96    else:
97            for i in range(dimlist[0]):
98                inlist.append(array[i])
99    return inlist
100
101def listOfListsCreator(inlist, dimlist):
102    """
103    Creates a list of lists with dimensions defined in dimlist (a list of integers).
104    """
105    if len(dimlist[1:])>0:
106        for i in range(dimlist[0]):
107            inlist.append([])
108            listOfListsCreator(inlist[i], array, dimlist[1:])
109    return inlist
110
111def getBestName(var): 
112    """
113    Returns the most appropriate variable name for a NASA Ames header
114    """
115    if hasattr(var, "id"): name=var.id
116    if hasattr(var, "name"): name=var.name
117    if hasattr(var, "long_name"): name=var.long_name
118    if hasattr(var, "standard_name"): name=var.standard_name
119    if hasattr(var, "units") and not re.match("^\s+$", var.units): name="%s (%s)" % (name, var.units)   
120    # Do a check to see units are not duplicated
121    match=re.match("(.*\(%s\)\s*)\(%s\)(.*)$" % (var.units, var.units), name)
122    if match:
123        name=match.groups()[0]+match.groups()[1]
124       
125    if name[-2:]=="()": name=name[:-2]
126    return name
127
128def getMissingValue(var):
129    if hasattr(var, "missing_value"): return var.missing_value
130    if hasattr(var, "fill_value"): return var.fill_value
131    if hasattr(var, "_FillValue"): return var._FillValue
132    return 1.E20
133
134
135def fixHeaderLength(file):
136    lines=open(file).readlines()[:200]
137    count=0
138    for line in lines:
139        count=count+1
140        if line=="###Data Section begins on the next line###\n":
141            break
142    # Now replace 1st line NLHEAD with count
143    firstline=lines[0]
144    newline="%s    %s" % (count, firstline.split()[-1])
145    for line in fileinput.input(file, inplace=1):
146        if line==firstline:
147          # write the changed line
148          sys.stdout.write(newline+"\n")
149        else:
150          # write the line unchanged
151          sys.stdout.write(line)
152
153
154def cdms2na(infilename, outfilename, variable=None, ffi="automatic", rules=None):
155    print "Reading data from: %s" % infilename
156    cdmsfile=cdms.open(infilename)
157    globals=cdmsfile.attributes
158    if not variable:
159        vars=[]
160        for var in cdmsfile.variables:
161            vars.append(cdmsfile(var))   
162    else:
163        vars=[cdmsfile(variable)]
164
165    # Re-order variables if they have the attribute 'nasa_ames_var_number'
166    orderedVars=[None]*1000
167    otherVars=[]
168    for var in vars:
169        varMetadata=cdmsfile[var]
170        if hasattr(varMetadata, "nasa_ames_var_number"):
171            num=varMetadata.nasa_ames_var_number
172            orderedVars[num]=var
173        else:
174            otherVars.append(var)
175   
176    vars=[]
177    for var in orderedVars:
178        if var!=None:
179            vars.append(var)
180           
181    vars=vars+otherVars
182   
183    builder=CdmsToNABuilder(vars, globals)
184    #print builder.naDict["X"]
185   
186    print "\nWriting output NASA Ames file: %s" % outfilename
187    general.openNAFile(outfilename, 'w', builder.naDict)
188    nlhead=fixHeaderLength(outfilename)   
189    print "\nNASA Ames file written successfully: %s" % outfilename
190
191
192class CdmsToNABuilder(NACore):
193    """
194    Class to build a NASA Ames File object from a set of
195    CDMS variables and global attributes (optional).
196    """
197   
198    def __init__(self, vars, global_attributes={}):
199        """
200        Sets up instance variables and calls appropriate methods to
201        generate sections of NASA Ames file object.
202        """
203        self.naDict={}
204        self.vars=vars
205        self.globals=global_attributes
206        (self.orderedVars, auxVars)=self.analyseVariables()
207        self.naDict["NLHEAD"]="-999"
208        self.defineNAVars(self.orderedVars)
209        self.defineNAAuxVars(auxVars)
210        self.defineNAGlobals()
211        self.defineNAComments()
212        self.defineGeneralHeader()
213        # Quick fudge
214        if self.naDict["FFI"]==1001: self.naDict["X"]=self.naDict["X"][0]
215
216    def analyseVariables(self):
217        """
218        Method to examine the content of CDMS variables to return
219        a tuple of two lists containing variables and auxiliary variables
220        for the NASA Ames file object.
221        """
222        # Get largest ranked variable as the one we use as standard
223        rank=0
224        count=0
225        for var in self.vars:
226            count=count+1
227            if var.rank()>rank:
228                rank=var.rank()
229                bestVar=var
230                bestVarIndex=count
231            elif var.rank()==rank:
232                if len(var.flat)>len(bestVar.flat):
233                    bestVar=var
234                    bestVarIndex=count
235               
236        vars4NA=[bestVar]
237        auxVars4NA=[]
238        shape=bestVar.shape
239        ndims=len(shape)
240        self.naDict["NIV"]=ndims
241
242        # Work out which File Format Index is appropriate
243        if ndims in (2,3,4):
244            self.naDict["FFI"]=10+(ndims*1000)
245        elif ndims>4:
246            raise "Cannot write variables defined against greater than 4 axes in NASA Ames format."
247        else:
248            if len(auxVars4NA)>0 or (self.naDict.has_key("NAUXV") and self.naDict["NAUXV"]>0):
249                self.naDict["FFI"]=1010
250            else:
251                self.naDict["FFI"]=1001
252        print self.naDict["FFI"]
253        axes=bestVar.getAxisList()
254       
255        # Get other variable info
256        for var in self.vars[:bestVarIndex]+self.vars[bestVarIndex:]:
257            if len(var.shape)!=ndims or var.shape!=shape: 
258                # Could it be an auxiliary variable
259                if len(var.shape)!=1:  continue
260                caxis=var.getAxis(0)
261                if compareAxes(axes[0], caxis)==0:  continue
262                # I think it is an auxiliary variable
263                auxVars4NA.append(var) 
264            else:
265                caxes=var.getAxisList()
266                for i in range(ndims):           
267                    if compareAxes(axes[i], caxes[i])==0:
268                        continue
269                # OK, I think they are compatible
270                vars4NA.append(var)
271               
272        # Re-order if they previously came from NASA Ames files (i.e. including
273        # the attribute 'nasa_ames_var_number')
274        orderedVars=[None]*1000
275        otherVars=[]
276        for var in vars4NA:
277            if hasattr(var, "nasa_ames_var_number"):
278                orderedVars[var.nasa_ames_var_number[0]]=var
279            else:
280                otherVars.append(var)
281        # Remake vars4NA now in order
282        vars4NA=[]
283        for var in orderedVars:
284            if var!=None: vars4NA.append(var)
285        vars4NA=vars4NA+otherVars
286
287        # Now re-order the Auxiliary variables if they previously came from NASA
288        # Ames files (i.e. including the attribute 'nasa_ames_aux_var_number')
289
290        orderedAuxVars=[None]*1000
291        otherAuxVars=[]
292        for var in auxVars4NA:
293            if hasattr(var, "nasa_ames_aux_var_number"):
294                orderedAuxVars[var.nasa_ames_aux_var_number[0]]=var
295            else:
296                otherAuxVars.append(var)
297        # Remake auxVars4NA now in order
298        auxVars4NA=[]
299        for var in orderedAuxVars:
300            if var!=None: auxVars4NA.append(var)
301        auxVars4NA=auxVars4NA+otherAuxVars     
302        return (vars4NA, auxVars4NA)
303
304    def defineNAVars(self, vars):
305        """
306        Method to define NASA Ames file object variables and their
307        associated metadata.
308        """
309        self.naDict["NV"]=len(vars)
310        self.naDict["VNAME"]=[]
311        self.naDict["VMISS"]=[]
312        self.naDict["VSCAL"]=[]
313        self.naDict["V"]=[]
314        for var in vars:
315            name=getBestName(var)
316            self.naDict["VNAME"].append(name)
317            miss=getMissingValue(var)
318            if type(miss)!=float:  miss=miss[0]
319            self.naDict["VMISS"].append(miss)
320            #print self.naDict["VMISS"]
321            self.naDict["VSCAL"].append(1)
322            # AND THE ARRAY
323            # Populate the variable list 
324            ######## NOTE - might not have to do this #####
325            ######## It  might handle writing from a Numeric array ########
326            self.naDict["V"].append(var._data)
327            #listOfListsCreator(inlist, var.shape)
328            #arrayToList(var, inlist)
329
330            if not self.naDict.has_key("X"):
331                self.naDict["NXDEF"]=[]
332                self.naDict["NX"]=[]
333                # Create independent variable information
334                self.ax0=var.getAxis(0)
335                self.naDict["X"]=[list(self.ax0._data_)]
336                self.naDict["XNAME"]=[getBestName(self.ax0)]
337                if len(self.ax0)==1:
338                    self.naDict["DX"]=[0]
339                else:
340                    incr=self.ax0[1]-self.ax0[0]
341                    for i in range(1, len(self.ax0)):
342                        if (self.ax0[i]-self.ax0[i-1])!=incr:
343                            self.naDict["DX"]=[0]
344                            break
345                    else:
346                        self.naDict["DX"]=[incr]
347                # Now sort the rest of the axes
348                for axis in var.getAxisList()[1:]:
349                    self.getAxisDefinition(axis)
350
351    def defineNAAuxVars(self, auxVars):
352        """
353        Method to define NASA Ames file object auxiliary variables and their
354        associated metadata.
355        """
356        self.naDict["NAUXV"]=len(auxVars)
357        self.naDict["ANAME"]=[]
358        self.naDict["AMISS"]=[]
359        self.naDict["ASCAL"]=[]
360        self.naDict["A"]=[]
361        for var in auxVars:
362            name=getBestName(var)
363            self.naDict["ANAME"].append(name)
364            miss=getMissingValue(var)
365            if type(miss)!=float:  miss=miss[0]
366            self.naDict["AMISS"].append(miss)
367            self.naDict["ASCAL"].append(1)
368            # AND THE ARRAY
369            # Populate the variable list 
370            ######## NOTE - might not have to do this #####
371            ######## It  might handle writing from a Numeric array ########
372            self.naDict["A"].append(var._data)
373            #listOfListsCreator(inlist, var.shape)
374            #arrayToList(var, inlist)     
375
376    def getAxisDefinition(self, axis):
377        """
378        Method to create the appropriate NASA Ames file object
379        items associated with an axis (independent variable in
380        NASA Ames).
381        """
382        self.naDict["NX"].append(len(axis))
383        self.naDict["XNAME"].append(getBestName(axis))
384        # If only one item in axis values
385        if len(axis)<3:
386            self.naDict["DX"].append(0)
387            self.naDict["NXDEF"].append(len(axis))
388            self.naDict["X"].append(list(axis._data_))       
389            return
390   
391        incr=axis[1]-axis[0]
392        for i in range(1, len(axis)):
393            if (axis[i]-axis[i-1])!=incr:
394                self.naDict["DX"].append(0)
395                self.naDict["NXDEF"].append(len(axis))
396                self.naDict["X"].append(list(axis._data_))
397                break
398        else:
399            self.naDict["DX"].append(incr)
400            self.naDict["NXDEF"].append(3)
401            self.naDict["X"].append(axis[:3])
402        return
403
404    def defineNAGlobals(self):
405        """
406        Maps CDMS (NetCDF) global attributes into NASA Ames Header fields.
407        """
408        # Get the global mapping dictionary
409        globalmap=cdmsMap.toNA
410        # Check if we should add to it with locally set rules
411        locGlobs=localRules.localGlobalAttributes
412        for att in locGlobs.keys():
413            if not globalmap.has_key(att):
414                globalmap[key]=locGlobs[key]
415
416        self.extra_comments=[[],[],[]]  # Normal comments, special comments, other comments
417        for key in self.globals.keys():
418            if key!="first_valid_date_of_data" and type(self.globals[key]) not in (str, float, int): continue
419            if key in globalmap.keys():
420                if key=="history":
421                    timestring=time.strftime("%Y-%m-%d %H:%M:%S", time.localtime(time.time()))
422                    self.history="History:\n\t%s - Converted to NASA Ames format using nappy-%s.\n\t%s" % (timestring, version.version, self.globals[key])
423                    self.history=self.history.split("\n") 
424                   
425                elif key=="institution":
426                    # If fields came from NA then extract appropriate fields.
427                    match=re.match(r"(.*)\s+\(ONAME from NASA Ames file\);\s+(.*)\s+\(ORG from NASA Ames file\)\.", self.globals[key])
428                    if match:
429                        self.naDict["ONAME"]=match.groups()[0]
430                        self.naDict["ORG"]=match.groups()[1]
431                    else:
432                        self.naDict["ONAME"]=self.globals[key]
433                        self.naDict["ORG"]=self.globals[key]               
434                                   
435                elif key=="comment":
436                    # Need to work out if they are actually comments from NASA Ames in the first place
437                    #self.ncom=[self.globals[key]]
438                    comLines=self.globals[key].split("\n")
439                    normComms=[]
440                    normCommFlag=None
441                    specComms=[]
442                    specCommFlag=None
443                    for line in comLines:
444                        if line.find("###NASA Ames Special Comments follow###")>-1:
445                            specCommFlag=1
446                        elif line.find("###NASA Ames Special Comments end###")>-1:
447                            specCommFlag=None
448                        elif line.find("###NASA Ames Normal Comments follow###")>-1:
449                            normCommFlag=1
450                        elif line.find("###NASA Ames Normal Comments end###")>-1:
451                            normCommFlag=None   
452                        elif specCommFlag==1:
453                            specComms.append(line)
454                        elif normCommFlag==1:
455                            normComms.append(line)
456                        elif line.find("###Data Section begins on the next line###")>-1:
457                            pass
458                        else:
459                            normComms.append(line)         
460                   
461                    self.extra_comments=[specComms, normComms, []]                 
462                                   
463                elif key=="first_valid_date_of_data":
464                    self.naDict["DATE"]=self.globals[key]
465               
466                elif key in ("Conventions", "references"):
467                    pass
468                else:
469                    self.naDict[globalmap[key]]=self.globals[key]
470            else:
471                self.extra_comments[2].append("%s:   %s" % (key, self.globals[key]))
472        return
473
474
475    def defineNAComments(self, normal_comments=[], special_comments=[]):
476        """
477        Defines the Special and Normal comments sections in the NASA Ames file
478        object - including information gathered from the defineNAGlobals method.
479        """
480       
481        if hasattr(self, "ncom"):  normal_comments=self.ncom+normal_comments
482        NCOM=[]
483        for ncom in normal_comments:
484            NCOM.append(ncom)
485        if len(NCOM)>0:   NCOM.append("")
486       
487        if len(self.extra_comments[2])>0:
488            for excom in extra_comments[2]:
489                NCOM.append(excom)
490       
491        if len(self.extra_comments[1])>0: 
492            NCOM.append("Additional Global Attributes defined in the source file and not translated elsewhere:")
493            for excom in self.extra_comments[1]:
494                NCOM.append(excom)
495
496        if hasattr(self, "history"):
497            for h in self.history:
498                NCOM.append(h)
499       
500        if len(NCOM)>0:
501            NCOM.insert(0, "###NASA Ames Normal Comments follow###")
502            NCOM.append("")
503            NCOM.append("###NASA Ames Normal Comments end###")
504        NCOM.append("###Data Section begins on the next line###")
505
506        specCommentsFlag=None
507        SCOM=[]
508        special_comments=self.extra_comments[0]
509        if len(special_comments)>0: 
510            SCOM=["###NASA Ames Special Comments follow###"]
511            specCommentsFlag=1
512        for scom in special_comments:
513            SCOM.append(scom)
514
515
516        used_var_atts=("name", "long_name", "standard_name", "id", 
517                "missing_value", "fill_value", "units", 
518                "nasa_ames_var_number", "nasa_ames_aux_var_number")
519        varCommentsFlag=None
520
521        for var in self.orderedVars:
522            varflag="unused"
523            name=getBestName(var)
524            for scom in var.attributes.keys():
525                if type(scom) in (str, float, int) and scom not in used_var_atts:
526                    if varflag=="unused":
527                        if varCommentsFlag==None:
528                            varCommentsFlag=1
529                            if specCommentsFlag==None:
530                                SCOM=["###NASA Ames Special Comments follow###"]
531                            SCOM.append("Additional Variable Attributes defined in the source file and not translated elsewhere:")
532                            SCOM.append("###Variable attributes from source (NetCDF) file follow###")
533                        varflag=="using" 
534                        SCOM.append("\tVariable (%s): %s" % (var.id, name))
535                    SCOM.append("\t\t%s = %s" % (scom, var.attributes[scom]))
536        if varCommentsFlag==1:  SCOM.append("###Variable attributes from source (NetCDF) file end###")
537        if specCommentsFlag==1:
538            SCOM.append("###NASA Ames Special Comments end###")
539
540        """used_var_atts=("name", "long_name", "standard_name", "id", "missing_value", "fill_value", "units")
541        for var in self.vars:
542            for scom in var.attributes.keys():
543                name=getBestName(var)
544                if type(scom) in (str, float, int) and scom not in used_var_atts:
545                    SCOM.append("\t%s: %s - %s" % (name, scom, var.attributes[scom]))"""
546
547        # Strip out empty lines (or returns)
548        NCOM_cleaned=[]
549        SCOM_cleaned=[]
550        for c in NCOM:
551            if c.strip() not in ("", " ", "  "): NCOM_cleaned.append(c)
552        for c in SCOM:
553            if c.strip() not in ("", " ", "  "): SCOM_cleaned.append(c)
554
555        self.naDict["NCOM"]=NCOM_cleaned
556        self.naDict["NNCOML"]=len(self.naDict["NCOM"])
557        self.naDict["SCOM"]=SCOM_cleaned
558        self.naDict["NSCOML"]=len(self.naDict["SCOM"])
559        return
560
561    def defineGeneralHeader(self, header_items={}):
562        """
563        Defines known header items and overwrites any with header_items
564        key/value pairs.
565        """
566        # Check if DATE field previously known in NASA Ames file
567        time_now=time.strftime("%Y %m %d", time.localtime(time.time())).split()
568        if not self.naDict.has_key("RDATE"):
569            self.naDict["RDATE"]=time_now
570       
571        if self.ax0.isTime():
572            # Get first date in list
573            try:
574                (unit, start_date)=re.match("(\w+)\s+?since\s+?(\d+-\d+-\d+)", self.ax0.units).groups()           
575                comptime=cdtime.s2c(start_date)
576                first_day=comptime.add(self.naDict["X"][0][0], getattr(cdtime, unit.capitalize()))
577                self.naDict["DATE"]=string.replace(str(first_day).split(" ")[0], "-", " ").split()
578            except:
579                print "Nappy Warning: Could not get the first date in the file. You will need to manually edit the output file."
580                self.naDict["DATE"]="<DATE_UNKNOWN>"
581        else: 
582            if not self.naDict.has_key("DATE"):
583                print "Nappy Warning: Could not get the first date in the file. You will need to manually edit the output file."
584                self.naDict["DATE"]="<DATE_UNKNOWN>"
585        self.naDict["IVOL"]=1
586        self.naDict["NVOL"]=1
587        for key in header_items.keys():
588             self.naDict[key]=header_items[key]
589        return
590               
591
592
593if __name__=="__main__":
594
595    args=sys.argv[1:]
596    if len(args)<4:
597        print helpMessage
598        print "Incorrect number of arguments used."
599        sys.exit()
600       
601    for arg in args:
602        if arg=="-i":
603            infile=args[args.index(arg)+1]
604        elif arg=="-o":
605            outfile=args[args.index(arg)+1]
606
607    cdms2na(infile, outfile) 
Note: See TracBrowser for help on using the repository browser.