source: nappy/trunk/bin/nc2na @ 352

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

Replaced sys.path lines to point to likely nappy directories. Fixed issue if no time axis defined.
Ag Stephens

  • 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.