source: nappy/trunk/bin/nc2na @ 347

Subversion URL: http://proj.badc.rl.ac.uk/svn/ndg/nappy/trunk/bin/nc2na@347
Revision 347, 20.6 KB checked in by selatham, 15 years ago (diff)

Incorporates changes by selatham and Ag Stephens. as of 5/11/04.

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