source: nappy/trunk/cdms2na.py @ 351

Subversion URL: http://proj.badc.rl.ac.uk/svn/ndg/nappy/trunk/cdms2na.py@351
Revision 351, 20.9 KB checked in by badc, 16 years ago (diff)

* empty log message *

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