source: nappy/tags/NDG0-1/cdms2na.py @ 353

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

Two bug fixes. One to deal with regular expression error if units
contained a regex special character. The other to stop the NetCDF-style
attributes including the variable name between each.

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