Changeset 26761


Ignore:
Timestamp:
01/10/22 02:58:43 (3 years ago)
Author:
bdef
Message:

CHG:update of ma netCDF export and related changes

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

Legend:

Unmodified
Added
Removed

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

    r25502 r26761  
    11function export_netCDF(md,filename)
    2 
    3 %Now going on Real treatment
     2%verbosity of the code, 0 is no messages, 5 is chatty
     3        verbose = 5;
    44        if exist(filename),
    5                 disp(sprintf('File %s allready exist', filename));
    6                 prompt = 'Give a new name or "delete" to replace: ';
    7                 newname = input(prompt,'s');
    8                 if strcmp(newname,'delete')
    9                         delete(filename)
    10                 else
    11                         disp(sprintf('New file name is %s ', newname));
    12                         filename=newname
    13           end
    14   end
     5                delete(filename)
     6                % disp(sprintf('File %s allready exist', filename));
     7                % prompt = 'Give a new name or "delete" to replace: ';
     8                % newname = input(prompt,'s');
     9                % if strcmp(newname,'delete')
     10                %       delete(filename)
     11                % else
     12                %       disp(sprintf('New file name is %s ', newname));
     13                %       filename=newname
     14                % end
     15        end
    1516        %open file and write description
    1617        mode = netcdf.getConstant('NC_NETCDF4');
    17         mode = bitor(mode,netcdf.getConstant('NC_NOCLOBBER'));%NOCLOBBER to avoid overwrite
     18        mode = bitor(mode,netcdf.getConstant('NC_NOCLOBBER')); %NOCLOBBER to avoid overwrite
    1819        ncid = netcdf.create(filename,mode);
    1920        netcdf.putAtt(ncid,netcdf.getConstant('NC_GLOBAL'),'description',['Results for run ' md.miscellaneous.name]);
    … …  
    2122
    2223        %gather geometry and timestepping as dimensions
    23         resfields=fieldnames(md.results);
    24         Duration=size(eval(['md.results. ' resfields{1} ]),2);
     24        if isempty(fieldnames(md.results)),
     25                %results as no field so no time is present
     26                Duration = 0;
     27        else
     28                resfields = fieldnames(md.results);
     29                Duration = size(eval(['md.results. ' resfields{1} ]),2);
     30        end
    2531        if Duration>0,
    26                 StepNum=Duration;
     32                StepNum = Duration;
    2733        else
    2834                StepNum=1;
    29   end
    30 
    31    dimlist=[2,md.mesh.numberofelements,md.mesh.numberofvertices,size(md.mesh.elements,2)];
    32 
    33         %define netcdf dimensions
    34         DimSize(1).index=netcdf.defDim(ncid,'Time',StepNum);
     35        end
     36        DimSize(1).index=netcdf.defDim(ncid,'Time',StepNum);  %time is the first dimension
    3537        [DimSize(1).name,DimSize(1).value]=netcdf.inqDim(ncid,DimSize(1).index);
    3638        DimValue(1)=DimSize(1).value;
    37         DimSize(2).index=netcdf.defDim(ncid,'UnLim',netcdf.getConstant('NC_UNLIMITED'));
     39        DimSize(2).index=netcdf.defDim(ncid,'UnLim',netcdf.getConstant('NC_UNLIMITED')); % we add an unlimited dimension if needed
    3840        [DimSize(2).name,DimSize(2).value]=netcdf.inqDim(ncid,DimSize(2).index);
    3941        DimValue(2)=DimSize(2).value;
     42        % adding mesh related dimensions
     43        dimlist=[2,40,md.mesh.numberofelements,md.mesh.numberofvertices,size(md.mesh.elements,2)];
     44        dimnames=["DictDummy" "StringLength" "EltNum" "VertNum" "VertPerElt"];
     45        if verbose > 0,
     46                disp('===Creating dimensions ===');
     47        end
     48        %define netcdf dimensions
    4049        for i=1:5
     50                % do not add the dimension if it exists already
    4151                if sum(dimlist(i) == DimValue) == 0
    42                         DimSize(i+2).index=netcdf.defDim(ncid,['Dimension' num2str(i+2)],dimlist(i));
     52                        DimSize(i+2).index=netcdf.defDim(ncid,dimnames(i),dimlist(i));
    4353                        [DimSize(i+2).name,DimSize(i+2).value]=netcdf.inqDim(ncid,DimSize(i+2).index);
    4454                        DimValue(i+2)=DimSize(i+2).value;
    4555                end
    4656        end
    47 
    48         typelist=[{'numeric'} {'logical'} {'string'} {'char'} {'cell'}];
    49 
     57        issmclasses = fieldnames(md)';
     58        typelist={'half', 'single','double','int8','int16'...
     59                  ,'int32','int64','uint8','uint16','uint32'...
     60                  ,'uint64','logical','char','string'};  %all malab types that are 0D
     61
     62        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}))];
     67                        end
     68                end
     69        end
    5070        %get all model classes and create respective groups
    5171        groups=fieldnames(md);
     72        if verbose > 0,
     73                disp('===Creating and populating groups===');
     74        end
    5275        for i=1:length(groups),
    53                 disp(sprintf('group name in tree %s ',groups{i}));
     76                if verbose >1,
     77                        disp(sprintf('===Now treating %s===',groups{i}));
     78                end
     79                if strcmp(groups{i}, 'qmu'),
     80                        disp('qmu is skipped until it is more stable');
     81                        continue
     82                end
    5483                groupID=netcdf.defGrp(ncid,groups{i});
    5584                %In each group gather the fields of the class
    56                 groupfields=fieldnames(md.(groups{i}));
    57                 for j=1:length(groupfields),
    58                         Var=md.(groups{i}).(groupfields{j});
     85                fields=fieldnames(md.(groups{i}));
     86                if isempty(fields),
     87                        disp(sprintf("WARNING: md.%s as no fields, we skip it.",groups{i}))
     88                        continue
     89                end
     90                %looping on fields in each group
     91                for j=1:length(fields),
     92                        Var=md.(groups{i}).(fields{j});
     93                        %treatment for lists
    5994                        if isa(Var,'cell')
    60                                 Stdlist=false;
     95                                Stdlist=false;  %first assume it is not a standard list
    6196                                if length(Var) == 0
    62                                         Stdlist=true;
     97                                        Stdlist=true;  %It is empty and so standard (for us)
    6398                                else
    6499                                        for k=1:length(typelist)
    65100                                                if isa(Var{1},typelist{k})
    66                                                         Stdlist=true;
     101                                                        Stdlist=true;  %if the list is of a known type (to matlab) if not it is probably some exotic ISSM stuff
    67102                                                end
    68103                                        end
    69104                                end
    70 
    71                                 netcdf.putAtt(groupID,netcdf.getConstant('NC_GLOBAL'),'classtype',class(md.(groups{i})));
    72                                 if(Stdlist)
    73                                         disp(sprintf('=====Field name in tree %s ',groupfields{j}));
    74                                         [DimSize,DimValue]=DefCreateVar(ncid,Var,groupID,groupfields{j},DimSize,DimValue);
    75                                 else
    76                                         listsize=length(Var);
    77                                         subgroupID=netcdf.defGrp(groupID,groupfields{j});
    78                                         netcdf.putAtt(subgroupID,netcdf.getConstant('NC_GLOBAL'),'classtype',class(Var));
    79                                         for l=1:listsize
    80                                                 if isprop(Var{l},'name')
    81                                                         lname=Var{l}.name;
    82                                                 elseif isprop(Var{l},'step')
    83                                                         lname=Var{l}.step
    84                                                 else
    85                                                         lname=[class(Var{l}) int2str(l)];
    86                                                 end
    87                                                 listgroupID=netcdf.defGrp(subgroupID,lname);
    88                                                 netcdf.putAtt(listgroupID,netcdf.getConstant('NC_GLOBAL'),'classtype',class(Var{l}));
    89                                                 subfields=fieldnames(Var{l});
    90                                                 for m=1:length(subfields)
    91                                                         if ~strcmp(subfields{m},'outlog')
    92                                                                 [DimSize,DimValue]=DefCreateVar(ncid,Var{l}.(subfields{m}),listgroupID,subfields{m},DimSize,DimValue);
     105                                %print the issm class as a classtype attribute
     106                                klass = class(md.(groups{i}));
     107                                klasstring = strcat(klass, '.',klass);
     108                                netcdf.putAtt(groupID,netcdf.getConstant('NC_GLOBAL'),'classtype',klasstring);
     109                                if(Stdlist)  % this is a standard or empty list just proceed
     110                                        if verbose > 4,
     111                                                disp(sprintf("=££=creating var for %s.%s with classtype : %s",groups{i}, fields{j}, klasstring))
     112                                        end
     113                                        if ~isempty(Var) && isa(Var{1}, 'char'),  % we have a char array, pad it to a given length
     114                                                Var=char(Var)';
     115                                        end
     116                                        [DimSize,DimValue,varid]=CreateVar(ncid,Var,groupID,fields{j},DimSize,DimValue);
     117                                        if ~isempty(varid),
     118                                                FillVar(Var,groupID,varid);
     119                                        end
     120
     121                                else % this is a list of fields, specific treatment needed (perhaps)
     122                                        if verbose > 4,
     123                                                disp(sprintf("=??=we have a list of fields for %s.%s with classtype : %s",groups{i}, fields{j}, klasstring));
     124                                        end
     125                                        if strcmp(groups{i}, 'outputdefinition'),
     126                                                listsize=length(Var);
     127                                                for k=1:listsize,
     128                                                        subgroupname=md.(groups{i}).(fields{j}){k}.definitionstring;
     129                                                        subgroupID=netcdf.defGrp(groupID,subgroupname);
     130                                                        klass=class(md.(groups{i}).(fields{j}){k});
     131                                                        klasstring = strcat(klass, '.',klass);
     132                                                        netcdf.putAtt(subgroupID,netcdf.getConstant('NC_GLOBAL'),'classtype',klasstring);
     133                                                        subfields=fieldnames(md.(groups{i}).(fields{j}){k});
     134                                                        for l=1:length(subfields)
     135                                                                if verbose > 4,
     136                                                                        disp(sprintf("=--=creating var for %s.%s[%i].%s",groups{i}, fields{j}, k, subfields{l}));
     137                                                                end
     138                                                                Var = md.(groups{i}).(fields{j}){k}.(subfields{l});
     139                                                                [DimSize,DimValue,varid]=CreateVar(ncid,Var,subgroupID,subfields{l},DimSize,DimValue);
     140                                                                if ~isempty(varid),
     141                                                                        FillVar(Var,subgroupID,varid);
     142                                                                end
    93143                                                        end
    94144                                                end
    95                                         end
    96                                 end
    97                         elseif isa(Var,'struct') && ~strcmp(groupfields{j},'bamg')
    98                                 classtype=class(md.(groups{i}));
    99                                 if strcmp(classtype,'struct')
    100                                         classtype=groups{i};
    101                                 end
    102                                 netcdf.putAtt(groupID,netcdf.getConstant('NC_GLOBAL'),'classtype',classtype);
    103                                 if length(Var)>1
    104                                         listsize=length(Var);
    105                                         subgroupID=netcdf.defGrp(groupID,groupfields{j});
    106                                         classtype=class(Var);
    107                                         if strcmp(classtype,'struct')
    108                                                 classtype=groups{i};
    109                                         end
    110                                         netcdf.putAtt(subgroupID,netcdf.getConstant('NC_GLOBAL'),'classtype',classtype);
    111                                         for l=1:listsize
    112                                                 if isfield(Var(l),'step')
    113                                                         lname=[int2str(Var(l).step)];
    114                                                 else
    115                                                         lname=[class(Var(l)) int2str(l)];
     145                                        else
     146                                                disp(sprintf("WARNING: unknown treatment for md.%s",groups{i}));
     147                                        end
     148                                end
     149                        elseif sum(strcmp(class(Var), typelist))==1, %this is a standard matlab class with no subgrouping
     150                                if verbose > 4,
     151                                        disp(sprintf("====creating var for %s.%s", groups{i}, fields{j}))
     152                                end
     153                                klass=class(md.(groups{i}));
     154                                klasstring = strcat(klass, '.',klass);
     155                                netcdf.putAtt(groupID,netcdf.getConstant('NC_GLOBAL'),'classtype',klasstring);
     156                                if sum(numel(Var) == size(Var)) == 0,  %this is a 2D array or more (and not a vector with dimension 2 = 1)
     157                                        Var = Var';
     158                                end
     159
     160                                [DimSize,DimValue,varid]=CreateVar(ncid,Var,groupID,fields{j},DimSize,DimValue);
     161                                if ~isempty(varid),
     162                                        FillVar(Var,groupID,varid);
     163                                end
     164
     165
     166                        elseif isa(Var,'struct')  % structures need special treatment
     167                                if strcmp(groups{i}, 'results'),
     168                                        klasstring='results.results';
     169                                        netcdf.putAtt(groupID,netcdf.getConstant('NC_GLOBAL'),'classtype',klasstring);
     170                                        subgroupname=fields{j};
     171                                        subgroupID=netcdf.defGrp(groupID,subgroupname);
     172                                        klasstring = 'results.solutionstep';
     173                                        netcdf.putAtt(subgroupID,netcdf.getConstant('NC_GLOBAL'),'classtype',klasstring);
     174                                        subfields=fieldnames(md.(groups{i}).(fields{j}));
     175                                        if isempty(subfields),
     176                                                disp(sprintf("WARNING: md.%s.%s as no subfields, we skip it.",groups{i}, fields{j}));
     177                                                continue
     178                                        end
     179                                        for k=1:length(subfields),
     180                                                if verbose > 4,
     181                                                        disp(sprintf("=@@=creating var for %s.%s.%s",groups{i}, fields{j}, subfields{k}));
    116182                                                end
    117                                                 classtype=class(Var(l));
    118                                                 if strcmp(classtype,'struct')
    119                                                         classtype=groups{i};
     183                                                Var = md.(groups{i}).(fields{j}).(subfields{k});
     184                                                [DimSize,DimValue,varid]=CreateVar(ncid,Var,subgroupID,subfields{k},DimSize,DimValue);
     185                                                if ~isempty(varid),
     186                                                        FillVar(Var,subgroupID,varid);
    120187                                                end
    121                                                 listgroupID=netcdf.defGrp(subgroupID,lname);
    122                                                 netcdf.putAtt(listgroupID,netcdf.getConstant('NC_GLOBAL'),'classtype',classtype);
    123                                                 subfields=fieldnames(Var(l));
    124                                                 for m=1:length(subfields)
    125                                                         if ~strcmp(subfields{m},'outlog')
    126                                                                 [DimSize,DimValue]=DefCreateVar(ncid,Var(l).(subfields{m}),listgroupID,subfields{m},DimSize,DimValue);
     188                                        end
     189                                elseif isempty(fieldnames(md.(groups{i}).(fields{j}))) % this is an empty struct, jus treat it as normal
     190                                        klass=class(md.(groups{i}));
     191                                        klasstring = strcat(klass, '.',klass);
     192                                        netcdf.putAtt(groupID,netcdf.getConstant('NC_GLOBAL'),'classtype',klasstring);
     193                                        if verbose > 4,
     194                                                disp(sprintf("=[]=creating var for %s.%s",groups{i}, fields{j}));
     195                                        end
     196
     197                                        [DimSize,DimValue,varid]=CreateVar(ncid,Var,groupID,fields{j},DimSize,DimValue);
     198                                        if ~isempty(varid),
     199                                                FillVar(Var,groupID,varid);
     200                                        end
     201
     202                                else
     203                                        disp(sprintf("WARNING, md.%s.%s is not treated as it does not fall in one of the existing cases with class '%s'.",groups{i}, fields{j}, class(md.(groups{i}).(fields{j}))))
     204                                end
     205                        % elseif strcmp(fields{j},'bamg')
     206                        %       klass=class(md.(groups{i}));
     207                        %       % if strcmp(klass,'struct')
     208                        %       %       klass=groups{i};
     209                        %       % end
     210                        %       netcdf.putAtt(groupID,netcdf.getConstant('NC_GLOBAL'),'classtype',klass);
     211                        %       if length(Var)>1
     212                        %               listsize=length(Var);
     213                        %               subgroupID=netcdf.defGrp(groupID,fields{j});
     214                        %               klass=class(md.(groups{i}).(fields{j}));
     215                        %               netcdf.putAtt(subgroupID,netcdf.getConstant('NC_GLOBAL'),'classtype',klass);
     216                        %               for k=1:listsize
     217                        %                       if isfield(Var(k),'step')
     218                        %                               lname=[int2str(Var(k).step)];
     219                        %                       else
     220                        %                               lname=[class(Var(k)) int2str(k)];
     221                        %                       end
     222
     223                        %                       listgroupID=netcdf.defGrp(subgroupID,lname);
     224                        %                       klass=class(md.(groups{i}).(fields{j}){k});
     225                        %                       netcdf.putAtt(listgroupID,netcdf.getConstant('NC_GLOBAL'),'classtype',klass);
     226                        %                       subfields=fieldnames(Var(k));
     227                        %                       for l=1:length(subfields)
     228                        %                               if ~strcmp(subfields{l},'outlog')
     229                        %                                       if verbose > 4,
     230                        %                                               disp(sprintf("=[]=creating var for %s.%s[%i].%s", groups{i}, fields{j}, k, subfields{l}))
     231                        %                                       end
     232                        %                                       Var=md.(groups{i}).(fields{j}){k}.(subfields{l});
     233                        %                                       [DimSize,DimValue,varid]=CreateVar(ncid,Var,listgroupID,subfields{l},DimSize,DimValue);
     234                        %                                       if ~isempty(varid),
     235                        %                                               FillVar(Var,listgroupID,varid);
     236                        %                                       end
     237                        %                               end
     238                        %                       end
     239                        %               end
     240                                % else
     241                                %       subgroupID=netcdf.defGrp(groupID,fields{j});
     242                                %       klass=class(md.(groups{i}));
     243                                %       netcdf.putAtt(subgroupID,netcdf.getConstant('NC_GLOBAL'),'classtype',klass);
     244                                %       subfields=fieldnames(Var);
     245                                %       for k=1:length(subfields)
     246                                %               if ~strcmp(subfields{k},'outlog')
     247                                %                       if verbose > 4,
     248                                %                               disp(sprintf("====creating var for %s.%s.%s", groups{i}, fields{j}, subfields{k}))
     249                                %                       end
     250                                %                       Var=md.(groups{i}).(fields{j}).(subfields{k});
     251                                %                       [DimSize,DimValue,varid]=CreateVar(ncid,Var,subgroupID,subfields{k},DimSize,DimValue);
     252                                %                       if ~isempty(varid),
     253                                %                               FillVar(Var,subgroupID,varid);
     254                                %                       end
     255                                %               end
     256                                %       end
     257                                % end
     258                        elseif sum(strcmp(class(Var), issmclasses)) == 1,  % that is an issm class
     259                                if strcmp(class(Var), 'solution'),
     260                                        if verbose > 4,
     261                                                disp(sprintf("=$$=creating var for %s.%s",groups{i}, fields{j}))
     262                                                disp("NEED treatment")
     263                                        end
     264                                elseif strcmp(class(Var), 'dict'),  %we have potential for a dict in py not to sure what it translates to here.
     265                                        if verbose > 4,
     266                                                disp(sprintf("=WW=creating var for %s.%s",groups{i}, fields{j}))
     267                                                disp("NEED Treatment")
     268                                        end
     269
     270                                else
     271                                        klass=class(md.(groups{i}));
     272                                        klasstring = strcat(klass, '.',klass);
     273                                        netcdf.putAtt(groupID,netcdf.getConstant('NC_GLOBAL'),'classtype',klasstring);
     274                                        subgroupID=netcdf.defGrp(groupID,fields{j});
     275                                        klass=class(md.(groups{i}).(fields{j}));
     276                                        klasstring = strcat(klass, '.',klass);
     277                                        netcdf.putAtt(subgroupID,netcdf.getConstant('NC_GLOBAL'),'classtype',klasstring);
     278                                        subfields=fieldnames(Var);
     279                                        for k=1:length(subfields),
     280                                                if sum(strcmp(subfields{k},["outlog" "errlog"])) == 0,
     281                                                        if verbose > 4,
     282                                                                disp(sprintf("+==+creating var for %s.%s.%s",groups{i}, fields{j}, subfields{k}))
     283                                                        end
     284                                                        Var=md.(groups{i}).(fields{j}).(subfields{k});
     285                                                        [DimSize,DimValue,varid]=CreateVar(ncid,Var,subgroupID,subfields{k},DimSize,DimValue);
     286                                                        if ~isempty(varid),
     287                                                                FillVar(Var,subgroupID,varid);
    127288                                                        end
    128289                                                end
    129290                                        end
    130                                 else
    131                                         subgroupID=netcdf.defGrp(groupID,groupfields{j});
    132                                         classtype=class(Var);
    133                                         if strcmp(classtype,'struct')
    134                                                 classtype=groups{i};
    135                                         end
    136                                         netcdf.putAtt(subgroupID,netcdf.getConstant('NC_GLOBAL'),'classtype',classtype);
    137                                         subfields=fieldnames(Var);
    138                                         for m=1:length(subfields)
    139                                                 if ~strcmp(subfields{m},'outlog')
    140                                                         [DimSize,DimValue]=DefCreateVar(ncid,Var.(subfields{m}),subgroupID,subfields{m},DimSize,DimValue);
    141                                                 end
    142                                         end
     291
    143292                                end
    144293                        else
    145                                 netcdf.putAtt(groupID,netcdf.getConstant('NC_GLOBAL'),'classtype',class(md.(groups{i})));
    146                                 [DimSize,DimValue]=DefCreateVar(ncid,Var,groupID,groupfields{j},DimSize,DimValue);
     294                                disp(sprintf("WARNING, md.%s.%s is not treated as it does not fall in one of the existing cases with class '%s'.",groups{i}, fields{j}, class(Var)))
    147295                        end
    148296                end
    149  end
    150  netcdf.close(ncid);
     297        end
     298        netcdf.close(ncid);
    151299end
    152300
    153 function [DimSize,DimValue]=DefCreateVar(ncid,Var,groupID,field,DimSize,DimValue,last,md,midfield)
    154         varclass=class(Var);
     301function [DimSize,DimValue,varid]=CreateVar(ncid,Var,groupID,field,DimSize,DimValue)
     302% Grab dimensions
    155303        varsize=size(Var);
    156304        varlength=length(Var);
     305        % treating scalar string or bool as atribute
    157306        if isa(Var,'logical'),
    158307                if Var,
    … …  
    160309                else,
    161310                        LogicString='False';
    162         end
     311                end
    163312                netcdf.putAtt(groupID,netcdf.getConstant('NC_GLOBAL'),field,LogicString);
     313                varid=[];
     314
    164315        elseif isa(Var,'char'),
    165                 netcdf.putAtt(groupID,netcdf.getConstant('NC_GLOBAL'),field,Var);
     316                if strcmp(field,'name'),  % it looks like netCDF does not like attributes that are called "name"
     317                        field = 'varname';
     318                end
     319                if size(Var,1) <= 1  %that is a single string or empty
     320                        netcdf.putAtt(groupID,netcdf.getConstant('NC_GLOBAL'),field,Var);
     321                        varid=[];
     322                else  % that is a character array
     323                        [dims,DimSize,DimValue]=GetDims(ncid,Var,DimSize,DimValue);
     324                        varid = netcdf.defVar(groupID,field,'NC_CHAR',dims);
     325                end
     326
    166327        elseif isa(Var,'double'), %dealing with arrays
     328                if all(mod(Var, 1) == 0, 'all')  %those are actually integers,
     329                        [dims,DimSize,DimValue]=GetDims(ncid,Var,DimSize,DimValue);
     330                        varid = netcdf.defVar(groupID,field,'NC_INT64',dims);
     331                else
     332                        [dims,DimSize,DimValue]=GetDims(ncid,Var,DimSize,DimValue);
     333                        varid = netcdf.defVar(groupID,field,'NC_DOUBLE',dims);
     334                end
     335        elseif isa(Var,'cell'),
    167336                [dims,DimSize,DimValue]=GetDims(ncid,Var,DimSize,DimValue);
    168                 varid = netcdf.defVar(groupID,field,'NC_DOUBLE',dims);
     337                varid = netcdf.defVar(groupID,field,'NC_CHAR',dims);
     338
     339        elseif isa(Var,'struct'),
     340                if isempty(fieldnames(Var)),
     341                        netcdf.putAtt(groupID,netcdf.getConstant('NC_GLOBAL'),field,'emptystruct');
     342                        varid=[];
     343                else
     344                        %Start by getting the structure fields and size
     345                        [dims,DimSize,DimValue]=GetDims(ncid,Var,DimSize,DimValue);
     346                        varid = netcdf.defVar(groupID,field,'NC_CHAR',dims);
     347                end
     348        else
     349                disp(sprintf('no support for class %s of field %s',class(Var),field));
     350                varid=[];
     351        end
     352        return
     353end
     354
     355
     356function FillVar(Var,groupID,varid)
     357% Grab dimensions
     358        varsize=size(Var);
     359        varlength=length(Var);
     360        % treating scalar string or bool as atribute
     361        if isa(Var,'double'), %dealing with arrays
     362                if all(mod(Var, 1) == 0, 'all')  %those are actually integers,
     363                        Var = int64(Var);
     364                end
    169365                if length(Var)==0,
    170366                        netcdf.putVar(groupID,varid,NaN);
    … …  
    172368                        netcdf.putVar(groupID,varid,Var);
    173369                end
     370        elseif isa(Var,'char'),  % at this point this should be a character array
     371                netcdf.putVar(groupID,varid,Var);
    174372        elseif isa(Var,'cell'),
    175                 [dims,DimSize,DimValue]=GetDims(ncid,Var,DimSize,DimValue);
    176                 %dirty hack to be able to pass strings
    177                 varid = netcdf.defVar(groupID,field,'NC_CHAR',dims);
    178                 if length(Var)==0,
    179                         netcdf.putVar(groupID,varid,0,9,'emptycell')
    180                 else
    181                         for i=1:length(Var),
    182                                 if length(Var)>1,
    183                                         endpoint=[1,min(length(Var{i}),40)];
    184                                         startpoint=[1 0];
    185                                 else
    186                                         endpoint=min(length(Var{i}),40);
    187                                         startpoint=0;
    188                                 end
    189                                 if length(Var{i})>40,
    190                                         netcdf.putVar(groupID,varid,startpoint,extent,Var{i}(1:40))
    191                                         disp(sprintf('some variable have been truncated'));
    192                                 else
    193                                         netcdf.putVar(groupID,varid,startpoint,endpoint,Var{i})
     373                if ~isempty(Var),
     374                        if length(Var)==0,
     375                                netcdf.putVar(groupID,varid,0,9,'emptycell');
     376                        else
     377                                for i=1:length(Var),
     378                                        if length(Var)>1,
     379                                                count=[min(length(Var{i}),40), 1];
     380                                                startpoint=[0 i-1];
     381                                        else
     382                                                count=min(length(Var{i}),40);
     383                                                startpoint=0;
     384                                        end
     385
     386                                        if length(Var{i})>40,
     387                                                netcdf.putVar(groupID,varid,startpoint,count,Var{i}(1:40));
     388                                                disp(sprintf('some variable have been truncated'));
     389                                        else
     390                                                netcdf.putVar(groupID,varid,startpoint,count,Var{i});
     391                                        end
    194392                                end
    195393                        end
    … …  
    197395        elseif isa(Var,'struct'),
    198396                %Start by getting the structure fields and size
    199                 locfields=fieldnames(Var);
    200                 [dims,DimSize,DimValue]=GetDims(ncid,Var,DimSize,DimValue);
    201                 varid = netcdf.defVar(groupID,field,'NC_CHAR',dims);
    202                 if length(locfields)==0,
    203                         netcdf.putVar(groupID,varid,[0,0],[1,11],'emptystruct')
    204                 else
    205                         for i=1:length(locfields),
    206                                 for j=1:2,
    207                                         if j==1,
    208                                                 CharVar=locfields{i};
    209                                                 if length(CharVar)==0
    210                                                         CharVar='emptystruct';
    211                                                 end
    212                                                 startpoint=[i-1,0,0];
     397                locfields=fieldnames(Var)
     398                for i=1:length(locfields),
     399                        for j=1:2,
     400                                if j==1,
     401                                        CharVar=locfields{i};
     402                                        if length(CharVar)==0
     403                                                CharVar='emptystruct';
     404                                        end
     405                                        startpoint=[i-1,0,0];
     406                                else
     407                                        if isa(Var.(locfields{i}),'char'),
     408                                                CharVar=Var.(locfields{i});
    213409                                        else
    214                                                 if isa(Var.(locfields{i}),'char'),
    215                                                         CharVar=Var.(locfields{i});
    216                                                 else
    217                                                         CharVar=num2str(Var.(locfields{i}));
    218                                                 end
    219                                                 if length(CharVar)==0
    220                                                         CharVar='emptystruct';
    221                                                 end
    222                                                 startpoint=[i-1,1,0];
    223                                         end
    224 
    225                                         extent=[1,1,min(length(CharVar),40)];
    226                                         if length(CharVar)>40,
    227                                                 netcdf.putVar(groupID,varid,startpoint,extent,CharVar(1:40))
    228                                                 disp(sprintf('some variable have been truncated'));
    229                                         else
    230                                                 netcdf.putVar(groupID,varid,startpoint,extent,CharVar)
    231                                         end
     410                                                CharVar=num2str(Var.(locfields{i}));
     411                                        end
     412                                        if length(CharVar)==0
     413                                                CharVar='emptystruct';
     414                                        end
     415                                        startpoint=[i-1,1,0];
     416                                end
     417
     418                                extent=[1,1,min(length(CharVar),40)];
     419                                if length(CharVar)>40,
     420                                        netcdf.putVar(groupID,varid,startpoint,extent,CharVar(1:40));
     421                                        disp(sprintf('some variable have been truncated'));
     422                                else
     423                                        netcdf.putVar(groupID,varid,startpoint,extent,CharVar);
    232424                                end
    233425                        end
    234426                end
    235427        else
    236                 disp(sprintf('no support for class %s of field %s',varclass,field));
    237   end
     428                disp(sprintf('no support for class %s',class(Var)));
     429        end
    238430        return
    239431end
    … …  
    241433function [dims,DimSize,DimValue]=GetDims(ncid,Var,DimSize,DimValue)
    242434        dims=[];
    243         if isa(Var,'cell'),
    244                 varsize=size(Var');
    245         elseif isa(Var,'struct')
     435        % if isa(Var,'cell'),
     436        %       varsize=size(Var');
     437        % else
     438        if isa(Var,'struct')
    246439                varsize=length(fieldnames(Var));
    247440        else
    248441                varsize=size(Var);
    249442        end
    250         dim=sum(varsize>1);
    251         if dim>0
    252                 for i=1:dim
    253                         indsize=find(varsize(i)==DimValue);
    254                         if length(indsize)>0
    255                                 dims=[dims DimSize(indsize).index];
    256                         else
    257                                 indsize=length(DimSize)+1;
    258                                 DimSize(indsize).index=netcdf.defDim(ncid,['Dimension' num2str(indsize)],varsize(i));
    259                                 [DimSize(indsize).name,DimSize(indsize).value]=netcdf.inqDim(ncid,DimSize(indsize).index);
    260                                 DimValue(indsize)=DimSize(indsize).value;
    261                                 dims=[dims DimSize(indsize).index];
     443
     444        % dim=sum(varsize>1);
     445        dim=ndims(Var);
     446        if dim>0,
     447                for i=1:dim,
     448                        if size(Var, i) >1
     449                                indsize=find(varsize(i)==DimValue);
     450                                if length(indsize)>0
     451                                        dims=[dims DimSize(indsize).index];
     452                                else
     453                                        indsize=length(DimSize)+1;
     454                                        DimSize(indsize).index=netcdf.defDim(ncid,['DimNum' num2str(indsize)],varsize(i));
     455                                        [DimSize(indsize).name,DimSize(indsize).value]=netcdf.inqDim(ncid,DimSize(indsize).index);
     456                                        DimValue(indsize)=DimSize(indsize).value;
     457                                        dims=[dims DimSize(indsize).index];
     458                                end
    262459                        end
    263460                end
    264461        end
    265         %if we have a cell variable we need to add a stringlength dimension
     462        %if we have an empty  cell variable we need to add a stringlength for the no data
     463        % if isa(Var,'cell'),
     464        %       if isempty(dims)
     465        %               dims=[DimSize(4).index];
     466        %       else
     467        %               dims=[DimSize(4).index dims];
     468        %       end
     469        % end
     470        % struct also need an extra dimension 2, but only if non empty
    266471        if isa(Var,'struct'),
    267                 if DimValue(3)~=2
    268                         if DimValue(2)~=2
    269                                 dims=[dims DimSize(1).index];
    270                         else
    271                                 dims=[dims DimSize(2).index];
    272                         end
    273                 else
    274                         dims=[dims DimSize(3).index];
    275                 end
    276         end
    277         if isa(Var,'cell') || isa(Var,'struct'),
    278                 if DimValue(2)~=40
    279                         dims=[dims DimSize(1).index];
    280                 else
    281                         dims=[dims DimSize(2).index];
    282                 end
     472                dims=[dims DimSize(3).index DimSize(4).index];
    283473        end
    284474end

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

    r26566 r26761  
    6060def export_netCDF(md, filename):  # {{{
    6161    #verbosity of the code, 0 is no messages, 5 is chatty
    62     verbose = 0
     62    verbose = 5
    6363    if path.exists(filename):
    6464        print('File {} allready exist'.format(filename))
    … …  
    8888    dimindex = 2
    8989    #add mesh related dimension that we know are needed
    90     dimlist = [2, md.mesh.numberofelements, md.mesh.numberofvertices, np.shape(md.mesh.elements)[1]]
    91     dimnames = ['DictDummy', 'EltNum', 'VertNum', 'VertPerElt']
     90    dimlist = [2, 40, md.mesh.numberofelements, md.mesh.numberofvertices, np.shape(md.mesh.elements)[1]]
     91    dimnames = ['DictDummy', 'StringLength', 'EltNum', 'VertNum', 'VertPerElt']
    9292    if verbose > 0:
    9393        print('===Creating dimensions ===')
    … …  
    121121        # looping on fields in each group
    122122        for field in fields:
     123            Var = md.__dict__[group].__dict__[field]
    123124            # Special treatment for list fields
    124             if type(md.__dict__[group].__dict__[field]) == list:
     125            if type(Var) == list:
    125126                StdList = False
    126                 if len(md.__dict__[group].__dict__[field]) == 0:
     127                if len(Var) == 0:
    127128                    StdList = True  #this is an empty list
    128129                else:
    129130                    #returns False for exotic types (typicaly results)
    130                     StdList = type(md.__dict__[group].__dict__[field][0]) in typelist
     131                    StdList = type(Var[0]) in typelist
    131132                klass = type(md.__dict__[group]).__module__ + '.' + type(md.__dict__[group]).__name__
    132133                NCgroup.__setattr__('classtype', klass)
    … …  
    134135                    if verbose > 4:
    135136                        print("=££=creating var for {}.{} with classtype : {}".format(group, field, klass))
    136                     Var = md.__dict__[group].__dict__[field]
    137137                    Var = SqueezeVar(Var)
    138138                    DimDict, ncvar = CreateVar(NCData, Var, field, NCgroup, DimDict)
    … …  
    141141                else:  # this is a list of fields, specific treatment needed (usually results or outputdefinitions)
    142142                    if verbose > 4:
    143                         print("list of fields happens for {}.{} with classtype : {}".format(group, field, klass))
    144                     Listsize = len(md.__dict__[group].__dict__[field])
     143                        print("=??=we have a list of fields for {}.{} with classtype : {}".format(group, field, klass))
     144                    Listsize = len(Var)
    145145                    if group == 'results':  #for results we reshape the datas following time rather than subgrouping
    146146                        Subgroup = NCgroup.createGroup(str(field))
    … …  
    172172                                    StackedVar.update(Var)
    173173                                if verbose > 4:
    174                                     print("=$$=creating var for {}.{}.{}".format(group, field, subfield))
     174                                    print("=@@=creating var for {}.{}.{}".format(group, field, subfield))
    175175                                    print("last index of the list is {}".format(lastindex))
    176176                                StackedVar = SqueezeVar(StackedVar.finalize(int(lastindex)))
    … …  
    289289                                FillVar(ncvar, Var)
    290290            else:
    291                 print("WARNING, md.{}.{} is not treated as it does not fall in one of the existing cases.".format(str(group), str(field)))
     291                print("WARNING, md.{}.{} is not treated as it does not fall in one of the existing cases.".format(group, field))
    292292
    293293    NCData.close()
    … …  
    329329        if field == 'name':  # it looks like netCDF does not like attributes that are called "name"
    330330            field = 'varname'
    331         Group.__setattr__(str(field).swapcase(), str(var))
     331        #Group.__setattr__(str(field).swapcase(), str(var))
     332        Group.__setattr__(str(field), str(var))
    332333        ncvar = None
    333334    # numpy array of strings
    … …  
    337338            if field == 'name':
    338339                field = 'varname'
    339             Group.__setattr__(str(field).swapcase(), str(var[0]))
     340            #Group.__setattr__(str(field).swapcase(), str(var[0]))
     341            Group.__setattr__(str(field), str(var[0]))
    340342            ncvar = None
    341343        else:

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

    r26567 r26761  
    77from re import findall, split
    88import shelve
    9 from netCDF4 import Dataset
     9from netCDF4 import Dataset, chartostring
    1010import numpy as np
    1111from model import *
    … …  
    100100                    #that is the current treatment
    101101                    #here we have a more NC approach with time being a dimension
    102                     listtype = split(r'\.', classtype[mod][0])[0]
     102                    listtype = split(r'\.', classtype[mod][0])[1]
     103                    print(listtype)
    103104                    if len(NCFile.dimensions['Time']) == 1:
    104105                        nvdict['md'].__dict__[classtree[mod][0]].__dict__[classtree[mod][1]] = getattr(classtype[mod][1], listtype)()
    … …  
    137138                            Tree = nvdict['md'].__dict__[classtree[mod][0]].__dict__[classtree[mod][1]][:]
    138139                    #}}}
    139                 elif classtype[mod][0] == 'massfluxatgate':  #this is for output definitions {{{
     140                elif classtype[mod][0] == 'massfluxatgate.massfluxatgate':  #this is for output definitions {{{
    140141                    defname = split('Output|[0-9]+', classtree[mod][1])[1] + 's'
    141142                    defindex = int(findall('[0-9]+', classtree[mod][1])[0])
    142                     nvdict['md'].__dict__[classtree[mod][0]].__dict__[defname].append(getattr(classtype[mod][1], classtype[mod][0])())
     143                    nvdict['md'].__dict__[classtree[mod][0]].__dict__[defname].append(getattr(classtype[mod][1], 'massfluxatgate')())
    143144                    Tree = nvdict['md'].__dict__[classtree[mod][0]].__dict__[defname][defindex - 1]
    144145                #}}}
    … …  
    223224                        else:
    224225                            if vardim == 0:  #that is a scalar
    225                                 if str(varval[0]) == '':  #no value
     226                                if str(varval[0]) == '' or str(varval[0]) == '--':  #no value
    226227                                    Tree.__dict__[str(var)] = []
    227228                                elif varval[0] == 'True':  #treatin bool
    … …  
    233234
    234235                            elif vardim == 1:  #that is a vector
     236                                if debug:
     237                                    print("   for variable {} type is {}".format(str(var), varval.dtype))
    235238                                if varval.dtype == str:
    236239                                    if varval.shape[0] == 1:
    … …  
    240243                                    else:
    241244                                        Tree.__dict__[str(var)] = [str(vallue) for vallue in varval[:]]
     245                                elif varval.dtype == "|S1":  #that is for matlab chararcter arrays
     246                                    stringlist = chartostring(varval[:])
     247                                    Tree.__dict__[str(var)] = [stringlist.tolist(), ]
    242248                                else:
    243249                                    try:
    … …  
    253259                            elif vardim == 2:
    254260                                #dealling with dict
     261                                if debug:
     262                                    print("   for variable {} type is {}".format(str(var), varval.dtype))
    255263                                if varval.dtype == str:  #that is for dictionaries
    256264                                    if any(varval[:, 0] == 'toolkit'):  #toolkit definition have to be first
    … …  
    263271                                        strings2 = [str(arg[1]) for arg in varval]
    264272                                        Tree.__dict__[str(var)] = OrderedDict(list(zip(strings1, strings2)))
     273                                elif varval.dtype == "|S1":  #that is for matlab chararcter arrays
     274                                    stringlist = chartostring(varval[:, :])
     275                                    stringlist = [string.strip() for string in stringlist]
     276                                    Tree.__dict__[str(var)] = stringlist
    265277                                else:
    266278                                    if type(Tree) == list:
    … …  
    282294                        print("      ==> treating attribute {}".format(attr))
    283295                    if attr != 'classtype':  #classtype is for treatment, don't get it back
    284                         attribute = str(attr).swapcase()  #there is a reason for swapcase, no sure what it isanymore
    285                         if attr == 'VARNAME':
     296                        # attribute = str(attr).swapcase()  #there is a reason for swapcase, no sure what it isanymore
     297                        # if attr == 'VARNAME':
     298                        #     attribute = 'name'
     299                        if attr == 'varname':
    286300                            attribute = 'name'
     301                        else:
     302                            attribute = attr
    287303                        if type(Tree) == list:
    288304                            if debug:
Note: See TracChangeset for help on using the changeset viewer.