Changeset 27264


Ignore:
Timestamp:
09/05/22 02:05:50 (3 years ago)
Author:
bdef
Message:

CHG:Some netCDF IO modifications

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

Legend:

Unmodified
Added
Removed

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

    r26968 r27264  
    4141        DimValue(2)=DimSize(2).value;
    4242        % adding mesh related dimensions
    43         dimlist=[2,40,md.mesh.numberofelements,md.mesh.numberofvertices,size(md.mesh.elements,2)];
     43        dimlist=[2 40 md.mesh.numberofelements md.mesh.numberofvertices size(md.mesh.elements,2)];
    4444        dimnames=["DictDummy" "StringLength" "EltNum" "VertNum" "VertPerElt"];
     45        if isprop(md.mesh, 'edges'),
     46                dimlist(end+1)=md.mesh.numberofedges;
     47                dimnames(end+1)="EdgeNum";
     48        else
     49                dimlist(end+1)=0;
     50                dimnames(end+1)="EdgeNum";
     51        end
    4552        if verbose > 0,
    4653                disp('===Creating dimensions ===');
    4754        end
    4855        %define netcdf dimensions
    49         for i=1:5
     56        for i=1:length(dimlist)
    5057                % do not add the dimension if it exists already
    5158                if sum(dimlist(i) == DimValue) == 0
    … …  
    6168
    6269        for cl=1:length(issmclasses),
    63                 subclasses=fieldnames(md.(issmclasses{cl}))';
    64                 for sc=1:length(subclasses),
    65                         if sum(strcmp(class(md.(issmclasses{cl}).(subclasses{sc})), typelist)) == 0,
    66                                 issmclasses = [issmclasses class(md.(issmclasses{cl}).(subclasses{sc}))];
     70                if isempty(md.(issmclasses{cl})),
     71                        disp(sprintf("md.%s is empty and will be left as default",issmclasses{cl}));
     72                else
     73                        subclasses=fieldnames(md.(issmclasses{cl}))';
     74                        for sc=1:length(subclasses),
     75                                if sum(strcmp(class(md.(issmclasses{cl}).(subclasses{sc})), typelist)) == 0,
     76                                        issmclasses = [issmclasses class(md.(issmclasses{cl}).(subclasses{sc}))];
     77                                end
    6778                        end
    6879                end
    … …  
    8394                groupID=netcdf.defGrp(ncid,groups{i});
    8495                %In each group gather the fields of the class
     96                if isempty(md.(groups{i})),
     97                        disp(sprintf("WARNING: md.%s is empty, we skip it.",groups{i}))
     98                        continue
     99                end
    85100                fields=fieldnames(md.(groups{i}));
    86101                if isempty(fields),
    … …  
    169184                        elseif isa(Var,'struct')  % structures need special treatment
    170185                                if strcmp(groups{i}, 'results'),
    171                                         klasstring='results.results';
     186                                        klasstring=strcat(groups{i} ,'.', groups{i});
    172187                                        netcdf.putAtt(groupID,netcdf.getConstant('NC_GLOBAL'),'classtype',klasstring);
    173188                                        Listsize= length(md.(groups{i}).(fields{j}));
    … …  
    212227                                                end
    213228                                        end
     229                                elseif strcmp(groups{i}, 'toolkits'),
     230                                        klasstring=strcat(groups{i} ,'.', groups{i});
     231                                        netcdf.putAtt(groupID,netcdf.getConstant('NC_GLOBAL'),'classtype',klasstring);
     232                                        if verbose > 4,
     233                                                disp(sprintf("=}{=creating var for %s.%s",groups{i}, fields{j}));
     234                                        end
     235
     236                                        [DimSize,DimValue,varid]=CreateVar(ncid,Var,groupID,fields{j},DimSize,DimValue);
     237                                        if ~isempty(varid),
     238                                                FillVar(Var,groupID,varid);
     239                                        end
     240
    214241                                elseif isempty(fieldnames(md.(groups{i}).(fields{j}))) % this is an empty struct, jus treat it as normal
    215242                                        klass=class(md.(groups{i}));
    … …  
    399426        elseif isa(Var,'struct'),
    400427                %Start by getting the structure fields and size
    401                 locfields=fieldnames(Var)
     428                locfields=fieldnames(Var);
    402429                for i=1:length(locfields),
    403430                        for j=1:2,
    404431                                if j==1,
    405                                         CharVar=locfields{i};
     432                                        CharVar=locfields{i}';
     433                                        disp(size(CharVar))
    406434                                        if length(CharVar)==0
    407435                                                CharVar='emptystruct';
    408436                                        end
    409                                         startpoint=[i-1,0,0];
     437                                        startpoint=[0,0,i-1]
    410438                                else
    411439                                        if isa(Var.(locfields{i}),'char'),
    412                                                 CharVar=Var.(locfields{i});
     440                                                CharVar=Var.(locfields{i})';
    413441                                        else
    414                                                 CharVar=num2str(Var.(locfields{i}));
     442                                                CharVar=num2str(Var.(locfields{i}))';
    415443                                        end
    416444                                        if length(CharVar)==0
    417445                                                CharVar='emptystruct';
    418446                                        end
    419                                         startpoint=[i-1,1,0];
    420                                 end
    421 
    422                                 extent=[1,1,min(length(CharVar),40)];
     447                                        startpoint=[0,1,i-1]
     448                                end
     449
     450                                extent=[min(length(CharVar),40), 1, 1]
    423451                                if length(CharVar)>40,
    424452                                        netcdf.putVar(groupID,varid,startpoint,extent,CharVar(1:40));
    … …  
    458486        if dim>0,
    459487                for i=1:alldim,
    460                         if size(Var, i)>1 || i>dim,  %we skip dimensions with zero lenght but want to add dimensions from cells
     488                        if size(Var, i)>1 || i>dim || isa(Var, 'struct'),  %we skip dimensions with zero lenght but want to add dimensions from cells
    461489                                indsize=find(varsize(i)==DimValue);
    462490                                if length(indsize)>0
    … …  
    478506        % struct also need an extra dimension 2, but only if non empty
    479507        if isa(Var,'struct'),
    480                 dims=[dims DimSize(3).index DimSize(4).index];
     508                dims=[DimSize(4).index DimSize(3).index, dims];
    481509        end
    482510end

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

    r26968 r27264  
    9999    dimlist = [2, 40, md.mesh.numberofelements, md.mesh.numberofvertices, np.shape(md.mesh.elements)[1]]
    100100    dimnames = ['DictDummy', 'StringLength', 'EltNum', 'VertNum', 'VertPerElt']
     101    try:
     102        dimlist = dimlist + [md.mesh.numberofedges]
     103        dimnames = dimnames + ['EdgeNum']
     104    except AttributeError:
     105        #no edges on this mesh, we fix it at 0
     106        dimlist += [0]
     107        dimnames += ['EdgeNum']
    101108    if verbose > 0:
    102109        print('===Creating dimensions ===')
    … …  
    315322        if val_type.startswith('<U'):
    316323            val_type = 'stringarray'
    317             print(var)
    318324    except AttributeError:
    319325        val_type = type(var)

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

    r26968 r27264  
    5555
    5656    timeindex = False
     57    SteadySols = ['ThermalSolution', 'HydrologySolution', 'StressbalanceSolution']
    5758
    5859    for key, value in kwargs.items():
    … …  
    105106                    #here we have a more NC approach with time being a dimension
    106107                    listtype = split(r'\.', classtype[mod][0])[1]
    107                     print(listtype)
    108                     if len(NCFile.dimensions['Time']) == 1:
     108                    try:
     109                        soltype = str(getattr(curclass, 'SolutionType'))
     110                    except AttributeError:
     111                        #might be an older format try that instead :
     112                        soltype = str(getattr(curclass, 'sOLUTIONtYPE'))
     113                    if len(NCFile.dimensions['Time']) == 1 or soltype in SteadySols:
    109114                        nvdict['md'].__dict__[classtree[mod][0]].__dict__[classtree[mod][1]] = getattr(classtype[mod][1], listtype)()
    110115                        Tree = nvdict['md'].__dict__[classtree[mod][0]].__dict__[classtree[mod][1]]
    … …  
    148153                    Tree = nvdict['md'].__dict__[classtree[mod][0]].__dict__[defname][defindex - 1]
    149154                #}}}
     155                elif classtype[mod][0] == 'collections.OrderedDict':  #Treating multiple toolkits {{{
     156                    nvdict['md'].__dict__[classtree[mod][0]].__dict__[classtree[mod][1]] = getattr(classtype[mod][1], 'OrderedDict')
     157                    Tree = nvdict['md'].__dict__[classtree[mod][0]].__dict__[classtree[mod][1]]
    150158                else:
    151159                    if verbose > 0:
    … …  
    184192            else:
    185193                groupclass = [curclass]
    186             #==== We deal with Variables {{{
    187194            for groupindex, listclass in enumerate(groupclass):
     195                try:
     196                    soltype = str(getattr(listclass, 'SolutionType'))
     197                except AttributeError:
     198                    soltype = 'NoSol'
     199                #==== We deal with Variables {{{
    188200                for var in listclass.variables:
    189201                    if not resname or var == resname:
    … …  
    217229                                        else:
    218230                                            print('table dimension greater than 3 not implemented yet')
     231                                    elif soltype in SteadySols:
     232                                        Tree.__dict__[str(var)] = varval[:].data
    219233                                    else:  #old format had step sorted in difeerent group so last group is last time
    220234                                        Tree[0].__dict__[str(var)] = varval[:].data
    221235                                else:
    222236                                    if NewFormat:
    223                                         incomplete = 'Time' not in varval.dimensions
     237                                        incomplete = 'Time' not in varval.dimensions and soltype not in SteadySols
    224238                                        if incomplete:
    225239                                            try:
    … …  
    230244                                                #just one step, so no dimension, we just put it on the first solutionstep
    231245                                                timelist = [0]
     246                                        elif soltype in SteadySols:
     247                                            timelist = [0]
    232248                                        else:
    233249                                            timelist = np.arange(0, len(NCFile.dimensions['Time']))
    234                                         for t in timelist:
    235                                             if verbose > 5:
    236                                                 print("filing step {} for {}".format(t, var))
    237                                             if vardim == 0:
    238                                                 Tree[t].__dict__[str(var)] = varval[:].data
    239                                             elif vardim == 1:
    240                                                 stepval = ma.masked_array(varval[t].data, mask=np.where(np.isnan(varval[t]), 1, 0))
    241                                                 Tree[t].__dict__[str(var)] = ma.compressed(stepval)
    242                                             elif vardim == 2:
    243                                                 stepval = ma.masked_array(varval[t, :].data, mask=np.where(np.isnan(varval[t, :]), 1, 0))
    244                                                 Tree[t].__dict__[str(var)] = ma.compressed(stepval)
    245                                             elif vardim == 3:
    246                                                 stepval = ma.masked_array(varval[t, :, :].data, mask=np.where(np.isnan(varval[t, :, :]), 1, 0))
    247                                                 Tree[t].__dict__[str(var)] = ma.compressed(stepval).reshape((stepval.count(0)[0], stepval.count(1)[0]))
    248                                             else:
    249                                                 print('table dimension greater than 3 not implemented yet')
     250                                        if soltype in SteadySols:
     251                                            Tree.__dict__[str(var)] = varval[:].data
     252                                        else:
     253                                            for t in timelist:
     254                                                if verbose > 5:
     255                                                    print("filing step {} for {}".format(t, var))
     256                                                if vardim == 0:
     257                                                    Tree[t].__dict__[str(var)] = varval[:].data
     258                                                elif vardim == 1:
     259                                                    stepval = ma.masked_array(varval[t].data, mask=np.where(np.isnan(varval[t]), 1, 0))
     260                                                    Tree[t].__dict__[str(var)] = ma.compressed(stepval)
     261                                                elif vardim == 2:
     262                                                    stepval = ma.masked_array(varval[t, :].data, mask=np.where(np.isnan(varval[t, :]), 1, 0))
     263                                                    Tree[t].__dict__[str(var)] = ma.compressed(stepval)
     264                                                elif vardim == 3:
     265                                                    stepval = ma.masked_array(varval[t, :, :].data, mask=np.where(np.isnan(varval[t, :, :]), 1, 0))
     266                                                    Tree[t].__dict__[str(var)] = ma.compressed(stepval).reshape((stepval.count(0)[0], stepval.count(1)[0]))
     267                                                else:
     268                                                    print('table dimension greater than 3 not implemented yet')
    250269                                    else:
    251270                                        if verbose > 0:
    … …  
    315334                                            Tree.__dict__[str(var)] = varval[:, :].data
    316335                                elif vardim == 3:
    317                                     Tree.__dict__[str(var)] = varval[:, :, :].data
     336                                    if varval.dtype == "|S1":  #that is for matlab chararcter arrays
     337                                        #most likely that is a toolkit dictionar so should be treated as such
     338                                        #first we convert the character table to strings
     339                                        stringtable = []
     340                                        for i in range(np.shape(varval)[0]):
     341                                            stringtable.append([chartostring(varval[i, 0, :]), chartostring(varval[i, 1, :])])
     342                                        stringtable = np.asarray(stringtable, dtype=str)
     343                                        Tree.__dict__[str(var)] = OrderedDict([('toolkit', str(varval[np.where(stringtable[:, 0] == 'toolkit')[0][0], 1]))])
     344                                        strings1 = [str(arg[0]) for arg in stringtable if arg[0] != 'toolkits']
     345                                        strings2 = [str(arg[1]) for arg in stringtable if arg[0] != 'toolkits']
     346                                        Tree.__dict__[str(var)].update(list(zip(strings1, strings2)))
     347                                    else:
     348                                        Tree.__dict__[str(var)] = varval[:, :, :].data
    318349                                else:
    319350                                    print('table dimension greater than 3 not implemented yet')
    … …  
    387418                if class_dict[classe][0] not in ['dict', 'list', 'cell']:
    388419                    modulename = split(r'\.', class_dict[classe][0])[0]
    389                     if modulename == "giaivins":
    390                         print("WARNING: module {} does not exist anymore and is skipped".format(modulename))
    391                     else:
     420                    try:
    392421                        class_dict[classe].append(import_module(modulename))
    393422                        class_tree[classe] = [group, ]
     423                    except ModuleNotFoundError:
     424                        print("WARNING: module {} does not exist anymore and is skipped".format(modulename))
     425
    394426            except AttributeError:
    395427                print(('group {} is empty'.format(group)))
Note: See TracChangeset for help on using the changeset viewer.