Changeset 21235


Ignore:
Timestamp:
09/28/16 07:45:35 (9 years ago)
Author:
bdef
Message:

NEW:adding netcd export import ability in python

Location:
issm/trunk-jpl/src/m
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/m/contrib/defleurian/netCDF/export_netCDF.py

    r21130 r21235  
    11from netCDF4 import Dataset, stringtochar
    2 import numpy
     2import numpy as np
    33import time
    44import collections
    … …  
    2626        Duration=md.timestepping.final_time-md.timestepping.start_time
    2727        if Duration>0 and md.timestepping.time_step*md.settings.output_frequency>0:
    28                 StepNum=Duration/(md.timestepping.time_step*md.settings.output_frequency)
     28                StepNum=int(Duration/(md.timestepping.time_step*md.settings.output_frequency))
    2929        else:
    3030                StepNum=1
    … …  
    3232        Dimension1=NCData.createDimension('Dimension1',md.mesh.numberofelements)
    3333        Dimension2=NCData.createDimension('Dimension2',md.mesh.numberofvertices)
    34         Dimension3=NCData.createDimension('Dimension3',numpy.shape(md.mesh.elements)[1])
     34        Dimension3=NCData.createDimension('Dimension3',np.shape(md.mesh.elements)[1])
    3535        Dimension4=NCData.createDimension('Dimension4',StepNum)
    3636        Dimension5=NCData.createDimension('Dimension5',2)
    … …  
    4040                                                 len(Dimension3):'Dimension3',
    4141                                                 len(Dimension4):'Dimension4',
    42                                                  len(Dimension5):'Dimension6'}
     42                                                 len(Dimension5):'Dimension5'}
    4343
    4444        #get all model classes and create respective groups
    … …  
    5454                                NCgroup.__setattr__('classtype', "results")
    5555                                Subgroup=NCgroup.createGroup(str(supfield))
    56                                 Subgroup.__setattr__('classtype',str(supfield))
     56                                #Subgroup.__setattr__('classtype',md.results.__dict__[supfield].__class__.__name__)
     57                                Subgroup.__setattr__('classtype',md.results.__class__.__name__)
    5758                                if type(md.results.__dict__[supfield])==list:#the solution have several timestep
    5859                                        #get last timesteps and output frequency
    59                                         last_step = numpy.size(md.results.__dict__[supfield])
     60                                        last_step = np.size(md.results.__dict__[supfield])
    6061                                        step_freq = md.settings.output_frequency
    6162                                        #grab first time step
    … …  
    7071                                        for field in subfields:
    7172                                                if str(field)!='errlog' and str(field)!='outlog' and str(field)!='SolutionType':
    72                                                         print 'Treating '+str(group)+'.'+str(supfield)+'.'+str(field)
    7373                                                        Var=md.results.__dict__[supfield].__dict__[field]
    7474                                                        DimDict=CreateVar(NCData,Var,field,NCgroup,DimDict,False)
    7575                                else:
    76                                         print 'Result format not suported'
     76                                        print ('Result format not suported')
    7777                else:
    7878                       
    7979                        for field in fields:
    80                                 print 'Treating ' +str(group)+'.'+str(field)
    8180                                NCgroup.__setattr__('classtype', md.__dict__[group].__class__.__name__)
    8281                                Var=md.__dict__[group].__dict__[field]
    … …  
    9695                val_shape=dict.keys(var)
    9796        except TypeError:
    98                 val_shape=numpy.shape(var)
     97                val_shape=np.shape(var)
    9998
    10099
    … …  
    104103                                                        'int64':'i8'}
    105104               
    106         val_dim=numpy.shape(val_shape)[0]
     105        val_dim=np.shape(val_shape)[0]
    107106        #Now define and fill up variable
    108107        #treating scalar string or bool as atribute
    109         if val_type==str or val_type==bool:
     108        if val_type==str or val_type==unicode or val_type==bool:
    110109                Group.__setattr__(str(field), str(var))
    111 
    112110        #treating list as string table
    113         #matlab does not recognise strings so we have to settle down with char arrays
    114111        elif val_type==list:
    115112                dimensions,DimDict=GetDim(NCData,var,val_shape,DimDict,val_dim,istime)
    … …  
    119116                                ncvar[elt] = var[elt]
    120117                        except IndexError:
    121                                 ncvar[0]= " "
    122                                 #treating bool tables as string tables
     118                                ncvar= []
     119        #treating bool tables as string tables
    123120        elif val_type=='bool':
    124121                dimensions,DimDict=GetDim(NCData,var,val_shape,DimDict,val_dim,istime)
    … …  
    126123                for elt in range(0,val_shape[0]):
    127124                        ncvar[elt] = str(var[elt])
    128                         #treating dictionaries as string tables of dim 2
     125        #treating dictionaries as tables of strings
    129126        elif val_type==collections.OrderedDict:
    130127                dimensions,DimDict=GetDim(NCData,var,val_shape,DimDict,val_dim,istime)
    … …  
    133130                        ncvar[elt,0]=dict.keys(var)[elt]
    134131                        ncvar[elt,1]=str(dict.values(var)[elt]) #converting to str to avoid potential problems
    135                         #Now dealing with numeric variables
     132        #Now dealing with numeric variables
    136133        else:
    137134                dimensions,DimDict=GetDim(NCData,var,val_shape,DimDict,val_dim,istime)
    138135                ncvar = Group.createVariable(str(field),TypeDict[val_type],dimensions,zlib=True)
    139                
    140136                if istime:
    141137                        last=step_args[0]
    142                         freq=step_args[1]
    143138                        md=step_args[2]
    144139                        supfield=step_args[3]
    145140                        vartab=var
    146                         for time in range(freq-1,last,freq):
     141                        for time in range(0,last,1):
    147142                                if time!=0:
    148143                                        timevar=md.results.__dict__[supfield].__getitem__(time).__dict__[field]
    149                                         print 'Treating results.'+str(supfield)+'.'+str(field)+' for time '+str(time)
    150                                         vartab=numpy.column_stack((vartab,timevar))
    151                                 #print numpy.shape(vartab)
     144                                        vartab=np.column_stack((vartab,timevar))
    152145                        try:
    153146                                ncvar[:,:]=vartab[:,:]
    … …  
    156149                else:
    157150                        try:
    158                                 nan_val=numpy.isnan(var)
     151                                nan_val=np.isnan(var)
    159152                                if nan_val.all():
    160153                                        ncvar [:] = 'NaN'
    … …  
    171164        #grab dimension
    172165        for dim in range(0,i): #loop on the dimensions
    173                 if type(shape[0])==int: 
     166                if type(shape[0])==int:
    174167                        try:
    175168                                output=output+[str(DimDict[shape[dim]])] #test if the dimension allready exist
    176169                        except KeyError: #if not create it
    177                                 if (shape[dim])>1:
     170                                if (shape[dim])>0:
    178171                                        index=len(DimDict)+1
    179172                                        NewDim=NCData.createDimension('Dimension'+str(index),(shape[dim]))
    180173                                        DimDict[len(NewDim)]='Dimension'+str(index)
    181174                                        output=output+[str(DimDict[shape[dim]])]
    182                                         #print 'Defining dimension ' +'Dimension'+str(index)
    183175                elif type(shape[0])==str:#dealling with a dictionnary
    184176                        try:
    185                                 output=[str(DimDict[numpy.shape(shape)[0]])]+['DictDim']
     177                                #dimension5 is 2 to treat with dict
     178                                output=[str(DimDict[np.shape(shape)[0]])]+['Dimension5']
    186179                        except KeyError:
    187180                                index=len(DimDict)+1
    188                                 NewDim=NCData.createDimension('Dimension'+str(index),numpy.shape(shape)[0])
     181                                NewDim=NCData.createDimension('Dimension'+str(index),np.shape(shape)[0])
    189182                                DimDict[len(NewDim)]='Dimension'+str(index)
    190                                 output=[str(DimDict[numpy.shape(dict.keys(var))[0]])]+['Dimension5']
    191                                 #print 'Defining dimension ' +'Dimension'+str(index)
     183                                output=[str(DimDict[np.shape(dict.keys(var))[0]])]+['Dimension5']
    192184                        break
    193185        if istime:

  • issm/trunk-jpl/src/m/io/loadmodel.py

    r17497 r21235  
    11from loadvars import loadvars
    22from whichdb import whichdb
     3from netCDF4 import Dataset
    34
    45def loadmodel(path):
    … …  
    1718                pass
    1819        else:
    19                 raise IOError("loadmodel error message: file '%s' does not exist" % path)
     20                try:
     21                        NCFile=Dataset(path,mode='r')
     22                        NCFile.close()
     23                        pass
     24                except RuntimeError:
     25                        raise IOError("loadmodel error message: file '%s' does not exist" % path)
     26                #       try:
     27        #recover model on file and name it md
     28        struc=loadvars(path)
     29        name=[key for key in struc.iterkeys()]
     30        if len(name)>1:
     31                raise IOError("loadmodel error message: file '%s' contains several variables. Only one model should be present." % path)
    2032
    21         try:
    22                 #recover model on file and name it md
    23                 struc=loadvars(path)
    24 
    25                 name=[key for key in struc.iterkeys()]
    26                 if len(name)>1:
    27                         raise IOError("loadmodel error message: file '%s' contains several variables. Only one model should be present." % path)
    28 
    29                 md=struc[name[0]]
    30                 return md
    31 
    32         except Exception as me:
    33                 print me
    34                 raise IOError("could not load model '%s'" % path)
    35 
     33        md=struc[name[0]]
     34        return md

  • issm/trunk-jpl/src/m/io/loadvars.py

    r15494 r21235  
    11import shelve
    22import os.path
     3import numpy as np
     4from netCDF4 import Dataset
     5from netCDF4 import chartostring
     6from os import path
    37from whichdb import whichdb
     8from model import *
    49
    510def loadvars(*args):
    … …  
    5156        if whichdb(filename):
    5257                print "Loading variables from file '%s'." % filename
    53         else:
    54                 raise IOError("File '%s' not found." % filename)
     58               
     59                my_shelf = shelve.open(filename,'r') # 'r' for read-only
     60                if nvdict:
     61                        for name in nvdict.iterkeys():
     62                                try:
     63                                        nvdict[name] = my_shelf[name]
     64                                        print "Variable '%s' loaded." % name
     65                                except KeyError:
     66                                        value = None
     67                                        print "Variable '%s' not found." % name
    5568
    56         my_shelf = shelve.open(filename,'r') # 'r' for read-only
    57 
    58         if nvdict:
    59                 for name in nvdict.iterkeys():
    60                         try:
     69                else:
     70                        for name in my_shelf.iterkeys():
    6171                                nvdict[name] = my_shelf[name]
    6272                                print "Variable '%s' loaded." % name
    63                         except KeyError:
    64                                 value = None
    65                                 print "Variable '%s' not found." % name
     73
     74                my_shelf.close()
    6675
    6776        else:
    68                 for name in my_shelf.iterkeys():
    69                         nvdict[name] = my_shelf[name]
    70                         print "Variable '%s' loaded." % name
     77                try:
     78                        NCFile=Dataset(filename,mode='r')
     79                        NCFile.close()
     80                        classtype,classtree=netCDFread(filename)
     81                        nvdict['md']=model()
     82                        module=map(__import__,dict.values(classtype))
     83                        for i,mod in enumerate(dict.keys(classtype)):
     84#                               print('treating md.{}'.format(mod))
     85                                if np.size(classtree[mod])>1:
     86                                        if classtree[mod][0]=='results':
     87                                                #treating results (Dimension4 is time)
     88                                                resdim=len(NCFile.dimensions['Dimension4'])
     89                                                curclass=NCFile.groups[classtree[mod][0]].groups[classtree[mod][1]]
     90                                                nvdict['md'].__dict__[classtree[mod][0]].__dict__[classtree[mod][1]] = [getattr(module[i],classtype[mod])()]
     91                                                if resdim>1:
     92                                                        for t in range(1,resdim):
     93                                                                nvdict['md'].__dict__[classtree[mod][0]].__dict__[classtree[mod][1]].append(getattr(module[i],classtype[mod])())
     94                                                Tree=nvdict['md'].__dict__[classtree[mod][0]].__dict__[classtree[mod][1]][:]
     95                                        else:
     96                                                curclass=NCFile.groups[classtree[mod][0]].groups[classtree[mod][1]]
     97                                                nvdict['md'].__dict__[classtree[mod][0]].__dict__[classtree[mod][1]] = getattr(module[i],classtype[mod])()
     98                                                Tree=nvdict['md'].__dict__[classtree[mod][0]].__dict__[classtree[mod][1]]
     99                                else:
     100                                        curclass=NCFile.groups[classtree[mod][0]]
     101                                        nvdict['md'].__dict__[mod] = getattr(module[i],classtype[mod])()
     102                                        Tree=nvdict['md'].__dict__[classtree[mod][0]]
     103                                for var in curclass.variables:
     104                                        #print('    treating {}'.format(var))
     105                                        varval=curclass.variables[str(var)]
     106                                        vardim=varval.ndim
     107                                        try:
     108                                                val_type=str(varval.dtype)
     109                                        except AttributeError:
     110                                                val_type=type(varval)
     111                                        if vardim==0:
     112                                                try:
     113                                                        Tree.__dict__[str(var)]=varval.getValue()
     114                                                        if varval.getValue()=='True':
     115                                                                Tree.__dict__[str(var)]=True
     116                                                        elif varval.getValue()=='False':
     117                                                                Tree.__dict__[str(var)]=False
     118                                                except IndexError:
     119                                                        Tree.__dict__[str(var)]=[]
     120                                        elif vardim==1:
     121                                                if varval.dtype==str:
     122                                                        if varval.shape==1:
     123                                                                Tree.__dict__[str(var)]=varval[0]
     124                                                        if 'True' in varval[:] or 'False' in varval[:]:
     125                                                                Tree.__dict__[str(var)]=np.asarray(varval[:],dtype=bool)
     126                                                else:
     127                                                        if classtree[mod][0]=='results' and resdim>1:
     128                                                                for t in range(0,resdim):
     129                                                                        Tree[t].__dict__[str(var)]=varval[t]
     130                                                        else:
     131                                                                Tree.__dict__[str(var)]=varval[:]
     132                                        elif vardim==2:
     133                                                #dealling with dict
     134                                                if varval.dtype==str:
     135                                                        Tree.__dict__[str(var)]=dict(zip(varval[:,0], varval[:,1]))
     136                                                else:
     137                                                        if classtree[mod][0]=='results' and resdim>1:
     138                                                                for t in range(0,resdim):
     139                                                                        Tree[t].__dict__[str(var)]=varval[:,t]
     140                                                        else:
     141                                                                Tree.__dict__[str(var)]=varval[:,:]
     142                                        elif vardim==3:
     143                                                if classtree[mod][0]=='results' and resdim>1:
     144                                                        for t in range(0,resdim):
     145                                                                Tree[t].__dict__[str(var)]=varval[:,:,t]
     146                                                else:
     147                                                        Tree.__dict__[str(var)]=varval[:,:,:]
     148                                        else:
     149                                                print 'table dimension greater than 3 not implemented yet'
     150                                for attr in curclass.ncattrs():
     151                                        #print('    treating {}'.format(attr))
     152                                        if classtree[mod][0]!='results': #no attributes in results
     153                                                Tree.__dict__[str(attr)]=curclass.getncattr(attr)
     154                                                if curclass.getncattr(attr)=='True':
     155                                                        Tree.__dict__[str(attr)]=True
     156                                                elif curclass.getncattr(attr)=='False':
     157                                                        Tree.__dict__[str(attr)]=False
    71158
    72         my_shelf.close()
     159                except RuntimeError:
     160                        raise IOError("File '%s' not found." % filename)
    73161
    74162        if   len(args) >= 2 and isinstance(args[1],(str,unicode)):    # (value)
    … …  
    83171                return nvdict
    84172
     173
     174def netCDFread(filename):
     175        def walktree(data):
     176                keys = data.groups.keys()
     177                for key in keys:
     178                        yield [str(key)]
     179                        for children in walktree(data.groups[str(key)]):
     180                                child=[str(key)]
     181                                child.append(str(children[0]))
     182                                yield child
     183        print ('Opening {} for reading '.format(filename))
     184        NCData=Dataset(filename, 'r')
     185        class_dict={}
     186        class_tree={}
     187       
     188        for children in walktree(NCData):
     189                classe=str(children[0])
     190                if np.size(children)>1:
     191                        for name in children[1:]:
     192                                classe=classe+'.'+name
     193                        class_dict[classe]=str(getattr(NCData.groups[children[0]].groups[children[1]],'classtype'))
     194                else:
     195                        class_dict[classe]=str(getattr(NCData.groups[classe],'classtype'))
     196                class_tree[classe]=children
     197
     198        return class_dict,class_tree
Note: See TracChangeset for help on using the changeset viewer.