source: nappy/trunk/cdms2na.py @ 339

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

Initial revision
  • 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 <infilename> <outfilename> [<options>]
14
15Where
16-----
17
18    infilename  - name of input file (NetCDF).
19    outfilename - name of output file (NASA Ames).
20    options     - list of options.
21"""
22
23# Imports from python standard library
24import sys, os, time, string, fileinput
25
26# Import cdat modules
27import cdms, cdtime
28
29# Import package modules
30from cdmsMap import cdmsMap
31import version
32import general
33from naCore import *
34
35def compareAxes(ax1, ax2):
36    """Takes 2 cmds axis objects returning 1 if they are essentially
37    the same and 0 if not."""
38    for axtype in ("time", "level", "latitude", "longitude"):
39        if cdms.axisMatches(ax1, axtype) and not cdms.axisMatches(ax2, axtype):
40            return 0
41
42    # Check ids
43    if ax1.id!=ax2.id: return 0
44    # Check lengths
45    if len(ax1)!=len(ax2): return 0
46    # Check values
47    if ax1._data_!=ax2._data_: return 0
48    # Check units
49    if ax1.units!=ax2.units: return 0
50    # OK, I think they are the same axis!
51    return 1
52
53def compareVariables(var1, var2):
54    """Compares two cdms variables to see if they are defined on identical
55    axes. Returns 1 if true, 0 if not."""
56    try:
57        for i in range(len(var1.getAxisList())):
58            ax1=var1.getAxis(i)
59            ax2=var2.getAxis(i)       
60            if compreAxes(ax1, ax2)==0:
61                return 0
62    except:
63        return 0
64    return 1
65
66def isAuxAndVar(avar, var):
67    """Compares two cdms variables and returns true if the first dimension of the
68    main variable is the only dimension the auxiliary variable is defined against."""
69    if len(avar.shape)>1:
70        return 0
71    auxax=avar.getAxis(0)
72    varax=var.getAxis(0)
73    return compareAxes(auxax, varax)
74
75def arrayToList(array, inlist):
76
77    """ Takes an n-dimensional Numeric array and converts it to an
78    n-dimensional list object.
79    """
80    dimlist=array.shape
81    if len(dimlist[1:])>0:
82            for i in range(dimlist[0]):
83                arrayToList(inlist[i], array[i])
84    else:
85            for i in range(dimlist[0]):
86                inlist.append(array[i])
87    return inlist
88
89def listOfListsCreator(inlist, dimlist):
90    """
91    Creates a list of lists with dimensions defined in dimlist (a list of integers).
92    """
93    if len(dimlist[1:])>0:
94        for i in range(dimlist[0]):
95            inlist.append([])
96            listOfListsCreator(inlist[i], array, dimlist[1:])
97    return inlist
98
99def getBestName(var): 
100    """Returns the most appropriate variable name for a NASA Ames header"""
101    if hasattr(var, "id"): name=var.id
102    if hasattr(var, "name"): name=var.name
103    if hasattr(var, "long_name"): name=var.long_name
104    if hasattr(var, "standard_name"): name=var.standard_name
105    if hasattr(var, "units"): name="%s (%s)" % (name, var.units)
106    return name
107
108def getMissingValue(var):
109    if hasattr(var, "missing_value"): return var.missing_value
110    if hasattr(var, "fill_value"): return var.fill_value
111    if hasattr(var, "_FillValue"): return var._FillValue
112    return 1.E40
113
114
115def fixHeaderLength(file):
116    lines=open(file).readlines()[:200]
117    count=0
118    for line in lines:
119        count=count+1
120        if line=="***Data Section begins on the next line***\n":
121            break
122    # Now replace 1st line NLHEAD with count
123    firstline=lines[0]
124    newline="%s    %s" % (count, firstline.split()[-1])
125    for line in fileinput.input(file, inplace=1):
126        if line==firstline:
127          # write the changed line
128          sys.stdout.write(newline+"\n")
129        else:
130          # write the line unchanged
131          sys.stdout.write(line)
132
133
134def cdms2na(infilename, outfilename, ffi="automatic", rules=None):
135    cdmsfile=cdms.open(infilename)
136    globals=cdmsfile.attributes
137    vars=[]
138    for var in cdmsfile.variables:
139        vars.append(cdmsfile(var))   
140    builder=CdmsToNABuilder(vars, globals)
141    #print builder.naDict["X"]
142    general.openNAFile(outfilename, 'w', builder.naDict)
143    nlhead=fixHeaderLength(outfilename)   
144    print "%s written successfully." % outfilename
145
146
147class CdmsToNABuilder(NACore):
148
149    def __init__(self, vars, global_attributes={}):
150        self.naDict={}
151        self.vars=vars
152        self.globals=global_attributes
153        (vars, auxVars)=self.analyseVariables()
154        self.naDict["NLHEAD"]="-999"
155        self.defineNAVars(vars)
156        self.defineNAAuxVars(auxVars)
157        self.defineNAGlobals()
158        self.defineNAComments()
159        self.defineGeneralHeader()
160        # Quick fudge
161        if self.naDict["FFI"]==1001: self.naDict["X"]=self.naDict["X"][0]
162
163    def analyseVariables(self):
164        vars4NA=[self.vars[0]]
165        auxVars4NA=[]
166        shape=self.vars[0].shape
167        ndims=len(shape)
168        self.naDict["NIV"]=ndims
169
170        # Work out which File Format Index is appropriate
171        if ndims in (2,3,4):
172            self.naDict["FFI"]=10+(ndims*1000)
173        elif ndims>4:
174            raise "Cannot write variables defined against greater than 4 axes in NASA Ames format."
175        else:
176            if len(auxVars4NA)>0 or (self.naDict.has_key("NAUXV") and self.naDict["NAUXV"]>0):
177                self.naDict["FFI"]=1010
178            else:
179                self.naDict["FFI"]=1001
180
181        axes=self.vars[0].getAxisList()
182        for var in self.vars[1:]:
183            if len(var.shape)!=ndims or var.shape!=shape: 
184                # Could it be an auxiliary variable
185                if len(var.shape)!=1:  continue
186                caxis=var.getAxis(0)
187                if compareAxes(axes[0], caxis)==0:  continue
188                # I think it is an auxiliary variable
189                auxVars4NA.append(var) 
190            else:
191                caxes=var.getAxisList()
192                for i in range(ndims):           
193                    if compareAxes(axes[i], caxes[i])==0:
194                        continue
195                # OK, I think they are compatible
196                vars4NA.append(var)
197
198        return (vars4NA, auxVars4NA)
199
200    def defineNAVars(self, vars):
201        self.naDict["NV"]=len(vars)
202        self.naDict["VNAME"]=[]
203        self.naDict["VMISS"]=[]
204        self.naDict["VSCAL"]=[]
205        self.naDict["V"]=[]
206        for var in vars:
207            name=getBestName(var)
208            self.naDict["VNAME"].append(name)
209            miss=getMissingValue(var)
210            if type(miss)!=float:  miss=miss[0]
211            self.naDict["VMISS"].append(miss)
212            #print self.naDict["VMISS"]
213            self.naDict["VSCAL"].append(1)
214            # AND THE ARRAY
215            # Populate the variable list 
216            ######## NOTE - might not have to do this #####
217            ######## It  might handle writing from a Numeric array ########
218            self.naDict["V"].append(var._data)
219            #listOfListsCreator(inlist, var.shape)
220            #arrayToList(var, inlist)
221
222            if not self.naDict.has_key("X"):
223                self.naDict["NXDEF"]=[]
224                self.naDict["NX"]=[]
225                # Create independent variable information
226                self.ax0=var.getAxis(0)
227                self.naDict["X"]=[list(self.ax0._data_)]
228                self.naDict["XNAME"]=[getBestName(self.ax0)]
229                if len(self.ax0)==1:
230                    self.naDict["DX"]=[0]
231                else:
232                    incr=self.ax0[1]-self.ax0[0]
233                    for i in range(1, len(self.ax0)):
234                        if (self.ax0[i]-self.ax0[i-1])!=incr:
235                            self.naDict["DX"]=[0]
236                            break
237                    else:
238                        self.naDict["DX"]=[incr]
239                # Now sort the rest of the axes
240                for axis in var.getAxisList()[1:]:
241                    self.getAxisDefinition(axis)
242
243    def defineNAAuxVars(self, auxVars):
244        self.naDict["NAUXV"]=len(auxVars)
245        self.naDict["ANAME"]=[]
246        self.naDict["AMISS"]=[]
247        self.naDict["ASCAL"]=[]
248        self.naDict["A"]=[]
249        for var in auxVars:
250            name=getBestName(var)
251            self.naDict["ANAME"].append(name)
252            miss=getMissingValue(var)
253            if type(miss)!=float:  miss=miss[0]
254            self.naDict["SMISS"].append(miss)
255            self.naDict["ASCAL"].append(1)
256            # AND THE ARRAY
257            # Populate the variable list 
258            ######## NOTE - might not have to do this #####
259            ######## It  might handle writing from a Numeric array ########
260            self.naDict["A"].append(var._data)
261            #listOfListsCreator(inlist, var.shape)
262            #arrayToList(var, inlist)     
263
264    def getAxisDefinition(self, axis):
265        self.naDict["NX"].append(len(axis))
266        self.naDict["XNAME"].append(getBestName(self.ax0))
267        incr=axis[1]-axis[0]
268        for i in range(1, len(axis)):
269            if (axis[i]-axis[i-1])!=incr:
270                self.naDict["DX"].append(0)
271                self.naDict["NXDEF"].append(len(axis))
272                self.naDict["X"].append(list(axis._data_))
273                break
274        else:
275            self.naDict["DX"].append(incr)
276            self.naDict["NXDEF"].append(3)
277            self.naDict["X"].append(axis[:3])
278        return
279
280    def defineNAGlobals(self):
281        globalmap=cdmsMap["toNA"]
282        self.extra_comments=[]
283        for key in self.globals.keys():
284            if type(self.globals[key]) not in (str, float, int): continue
285            if key in globalmap.keys():
286                if key=="history":
287                    timestring=time.strftime("%Y-%m-%d %H:%M:%S", time.localtime(time.time()))
288                    self.history="History:\n\t%s - Converted to NASA Ames format using nappy-%s.\n\t%s" % (timestring, version.version, self.globals[key])
289                    self.history=self.history.split("\n") 
290                   
291                elif key=="institution":
292                    self.naDict["ONAME"]=self.globals[key]
293                    self.naDict["ORG"]=self.globals[key]
294                elif key=="comment":
295                    self.ncom=[self.globals[key]]
296                elif key in ("Conventions", "references"):
297                    pass
298                else:
299                    self.naDict[globalmap[key]]=self.globals[key]
300            else:
301                self.extra_comments.append("%s:   %s" % (key, self.globals[key]))
302        return
303
304
305    def defineNAComments(self, normal_comments=[], special_comments=[]):
306        NCOM=["***NASA Ames Normal Comments follow***"]
307        if hasattr(self, "ncom"):  normal_comments=self.ncom+normal_comments
308        for ncom in normal_comments:
309            NCOM.append(ncom)
310        NCOM.append("")
311        NCOM.append("Additional Global Attributes defined in Cdmsfile and not translated elsewhere:")
312
313        if hasattr(self, "extra_comments"):
314            for excom in self.extra_comments:
315                NCOM.append(excom)
316
317        if hasattr(self, "history"):
318            for h in self.history:
319                NCOM.append(h)
320
321        NCOM.append("")
322        NCOM.append("***Data Section begins on the next line***")
323
324        SCOM=["***NASA Ames Special Comments follow***"]
325        for scom in special_comments:
326            SCOM.append(scom)
327
328        SCOM.append("Additional Variable Attributes defined in Cdmsfile and not translated elsewhere:")
329        used_var_atts=("name", "long_name", "standard_name", "id", "missing_value", "fill_value", "units")
330        for var in self.vars:
331            for scom in var.attributes.keys():
332                name=getBestName(var)
333                if type(scom) in (str, float, int) and scom not in used_var_atts:
334                    SCOM.append("\t%s: %s - %s" % (name, scom, var.attributes[scom]))
335
336        self.naDict["NCOM"]=NCOM
337        self.naDict["NNCOML"]=len(NCOM)
338        self.naDict["SCOM"]=SCOM
339        self.naDict["NSCOML"]=len(SCOM)
340        return
341
342    def defineGeneralHeader(self, header_items={}):
343        """Defines known header items and overwrites any with header_items key/value pairs."""
344        self.naDict["DATE"]=time.strftime("%Y %m %d", time.localtime(time.time())).split()
345        if self.ax0.isTime():
346            # Get first date in list
347            (unit, start_date)=re.match("(\w+)\s+?since\s+?(\d+-\d+-\d+)", self.ax0.units).groups()
348           
349            comptime=cdtime.s2c(start_date)
350            first_day=comptime.add(self.naDict["X"][0][0], getattr(cdtime, unit.capitalize()))
351            self.naDict["RDATE"]=string.replace(str(first_day).split(" ")[0], "-", " ").split()
352        else: 
353            self.naDict["RDATE"]=self.naDict["DATE"]
354        self.naDict["IVOL"]=1
355        self.naDict["NVOL"]=1
356        for key in header_items.keys():
357             self.naDict[key]=header_items[key]
358        return
359               
360
361
362if __name__=="__main__":
363
364    args=sys.argv[1:]
365    if len(args)!=2:
366        print helpMessage
367        print "Incorrect number of arguments used."
368        sys.exit()
369    else:
370        cdms2na(args[0], args[1]) 
Note: See TracBrowser for help on using the repository browser.