Index: /issm/trunk-jpl/src/m/contrib/defleurian/netCDF/export_netCDF.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/defleurian/netCDF/export_netCDF.m (revision 27890) +++ /issm/trunk-jpl/src/m/contrib/defleurian/netCDF/export_netCDF.m (revision 27890) @@ -0,0 +1,515 @@ +function export_netCDF(md,filename) +%verbosity of the code, 0 is no messages, 5 is chatty + verbose = 5; + if exist(filename), + delete(filename) + % disp(sprintf('File %s allready exist', filename)); + % prompt = 'Give a new name or "delete" to replace: '; + % newname = input(prompt,'s'); + % if strcmp(newname,'delete') + % delete(filename) + % else + % disp(sprintf('New file name is %s ', newname)); + % filename=newname + % end + end + %open file and write description + mode = netcdf.getConstant('NC_NETCDF4'); + mode = bitor(mode,netcdf.getConstant('NC_NOCLOBBER')); %NOCLOBBER to avoid overwrite + ncid = netcdf.create(filename,mode); + netcdf.putAtt(ncid,netcdf.getConstant('NC_GLOBAL'),'description',['Results for run ' md.miscellaneous.name]); + netcdf.putAtt(ncid,netcdf.getConstant('NC_GLOBAL'),'history',['Created ' datestr(now)]); + + %gather geometry and timestepping as dimensions + if isempty(fieldnames(md.results)), + %results as no field so no time is present + Duration = 0; + else + resfields = fieldnames(md.results); + Duration = size(eval(['md.results. ' resfields{1} ]),2); + end + if Duration>0, + StepNum = Duration; + else + StepNum=1; + end + DimSize(1).index=netcdf.defDim(ncid,'Time',StepNum); %time is the first dimension + [DimSize(1).name,DimSize(1).value]=netcdf.inqDim(ncid,DimSize(1).index); + DimValue(1)=DimSize(1).value; + DimSize(2).index=netcdf.defDim(ncid,'UnLim',netcdf.getConstant('NC_UNLIMITED')); % we add an unlimited dimension if needed + [DimSize(2).name,DimSize(2).value]=netcdf.inqDim(ncid,DimSize(2).index); + DimValue(2)=DimSize(2).value; + % adding mesh related dimensions + dimlist=[2 40 md.mesh.numberofelements md.mesh.numberofvertices size(md.mesh.elements,2)]; + dimnames=["DictDummy" "StringLength" "EltNum" "VertNum" "VertPerElt"]; + if isprop(md.mesh, 'edges'), + dimlist(end+1)=md.mesh.numberofedges; + dimnames(end+1)="EdgeNum"; + else + dimlist(end+1)=0; + dimnames(end+1)="EdgeNum"; + end + if verbose > 0, + disp('===Creating dimensions ==='); + end + %define netcdf dimensions + for i=1:length(dimlist) + % do not add the dimension if it exists already + if sum(dimlist(i) == DimValue) == 0 + DimSize(i+2).index=netcdf.defDim(ncid,dimnames(i),dimlist(i)); + [DimSize(i+2).name,DimSize(i+2).value]=netcdf.inqDim(ncid,DimSize(i+2).index); + DimValue(i+2)=DimSize(i+2).value; + end + end + issmclasses = fieldnames(md)'; + typelist={'half', 'single','double','int8','int16'... + ,'int32','int64','uint8','uint16','uint32'... + ,'uint64','logical','char','string'}; %all malab types that are 0D + + for cl=1:length(issmclasses), + if isempty(md.(issmclasses{cl})), + disp(sprintf("md.%s is empty and will be left as default",issmclasses{cl})); + else + subclasses=fieldnames(md.(issmclasses{cl}))'; + for sc=1:length(subclasses), + if sum(strcmp(class(md.(issmclasses{cl}).(subclasses{sc})), typelist)) == 0, + issmclasses = [issmclasses class(md.(issmclasses{cl}).(subclasses{sc}))]; + end + end + end + end + %get all model classes and create respective groups + groups=fieldnames(md); + if verbose > 0, + disp('===Creating and populating groups==='); + end + for i=1:length(groups), + if verbose >1, + disp(sprintf('===Now treating %s===',groups{i})); + end + if any(strcmp(groups{i}, {'qmu'})), + disp('qmu is skipped until it is more stable'); + continue + end + if any(strcmp(groups{i},{'radaroverlay'})), + disp(sprintf('%s is skipped.',groups{i})); + continue + end + groupID=netcdf.defGrp(ncid,groups{i}); + %In each group gather the fields of the class + if isempty(md.(groups{i})), + disp(sprintf("WARNING: md.%s is empty, we skip it.",groups{i})) + continue + end + fields=fieldnames(md.(groups{i})); + if isempty(fields), + disp(sprintf("WARNING: md.%s as no fields, we skip it.",groups{i})) + continue + end + %looping on fields in each group + for j=1:length(fields), + Var=md.(groups{i}).(fields{j}); + %treatment for lists + if isa(Var,'cell') + Stdlist=false; %first assume it is not a standard list + if length(Var) == 0 + Stdlist=true; %It is empty and so standard (for us) + else + for k=1:length(typelist) + if isa(Var{1},typelist{k}) + Stdlist=true; %if the list is of a known type (to matlab) if not it is probably some exotic ISSM stuff + end + end + end + %print the issm class as a classtype attribute + klass = class(md.(groups{i})); + klasstring = strcat(klass, '.',klass); + netcdf.putAtt(groupID,netcdf.getConstant('NC_GLOBAL'),'classtype',klasstring); + if(Stdlist) % this is a standard or empty list just proceed + if verbose > 4, + disp(sprintf("=££=creating var for %s.%s with classtype : %s",groups{i}, fields{j}, klasstring)) + end + if ~isempty(Var) && isa(Var{1}, 'char'), % we have a char array, pad it to a given length + Var=char(Var)'; + end + [DimSize,DimValue,varid]=CreateVar(ncid,Var,groupID,fields{j},DimSize,DimValue); + if ~isempty(varid), + FillVar(Var,groupID,varid); + end + + else % this is a list of fields, specific treatment needed (perhaps) + if verbose > 4, + disp(sprintf("=??=we have a list of fields for %s.%s with classtype : %s",groups{i}, fields{j}, klasstring)); + end + if strcmp(groups{i}, 'outputdefinition'), + listsize=length(Var); + for k=1:listsize, + subgroupname=md.(groups{i}).(fields{j}){k}.definitionstring; + subgroupID=netcdf.defGrp(groupID,subgroupname); + klass=class(md.(groups{i}).(fields{j}){k}); + klasstring = strcat(klass, '.',klass); + netcdf.putAtt(subgroupID,netcdf.getConstant('NC_GLOBAL'),'classtype',klasstring); + subfields=fieldnames(md.(groups{i}).(fields{j}){k}); + for l=1:length(subfields) + if verbose > 4, + disp(sprintf("=--=creating var for %s.%s[%i].%s",groups{i}, fields{j}, k, subfields{l})); + end + Var = md.(groups{i}).(fields{j}){k}.(subfields{l}); + if sum(numel(Var) == size(Var)) == 0, %this is a 2D array or more (and not a vector with dimension 2 = 1) + Var = Var'; + end + [DimSize,DimValue,varid]=CreateVar(ncid,Var,subgroupID,subfields{l},DimSize,DimValue); + if ~isempty(varid), + FillVar(Var,subgroupID,varid); + end + end + end + else + disp(sprintf("WARNING: unknown treatment for md.%s",groups{i})); + end + end + elseif sum(strcmp(class(Var), typelist))==1, %this is a standard matlab class with no subgrouping + if verbose > 4, + disp(sprintf("====creating var for %s.%s", groups{i}, fields{j})) + end + klass=class(md.(groups{i})); + klasstring = strcat(klass, '.',klass); + netcdf.putAtt(groupID,netcdf.getConstant('NC_GLOBAL'),'classgroup',groups{i}); + netcdf.putAtt(groupID,netcdf.getConstant('NC_GLOBAL'),'classtype',klasstring); + if sum(numel(Var) == size(Var)) == 0, %this is a 2D array or more (and not a vector with dimension 2 = 1) + Var = Var'; + end + + [DimSize,DimValue,varid]=CreateVar(ncid,Var,groupID,fields{j},DimSize,DimValue); + if ~isempty(varid), + FillVar(Var,groupID,varid); + end + + + elseif isa(Var,'struct') % structures need special treatment + if strcmp(groups{i}, 'results'), + klasstring=strcat(groups{i} ,'.', groups{i}); + netcdf.putAtt(groupID,netcdf.getConstant('NC_GLOBAL'),'classtype',klasstring); + Listsize= length(md.(groups{i}).(fields{j})); + subgroupname=fields{j}; + subgroupID=netcdf.defGrp(groupID,subgroupname); + klasstring='results.solutionstep'; + netcdf.putAtt(subgroupID,netcdf.getConstant('NC_GLOBAL'),'classtype',klasstring); + subfields=fieldnames(md.(groups{i}).(fields{j})); + if isempty(subfields), + disp(sprintf("WARNING: md.%s.%s as no subfields, we skip it.",groups{i}, fields{j})); + continue + end + for k=1:length(subfields), + if ~ismember(subfields{k}, {'errlog', 'outlog', 'SolutionType'}) + StackedVar=restable(); + for l=1:Listsize, + Var = md.(groups{i}).(fields{j})(l).(subfields{k}); + if length(Var) == 0, + %Some variables only have data on the first step + break + end + lastindex=l; + StackedVar=StackedVar.update(Var); + end + if verbose > 4, + disp(sprintf("=$$=creating var for %s.%s.%s",groups{i}, fields{j}, subfields{k})); + disp(sprintf("last index on the list is %i",lastindex)); + end + StackedVar=StackedVar.finalize(lastindex); + [DimSize,DimValue,varid]=CreateVar(ncid,StackedVar,subgroupID,subfields{k},DimSize,DimValue); + if ~isempty(varid), + FillVar(StackedVar,subgroupID,varid); + end + elseif ismember(subfields{k}, {'SolutionType'}) + %We just add solution type once as an attribute + Var = md.(groups{i}).(fields{j})(1).(subfields{k}); + [DimSize,DimValue,varid]=CreateVar(ncid,Var,subgroupID,subfields{k},DimSize,DimValue); + if ~isempty(varid), + FillVar(Var,subgroupID,varid); + end + + end + end + elseif strcmp(groups{i}, 'toolkits'), + klasstring=strcat(groups{i} ,'.', groups{i}); + netcdf.putAtt(groupID,netcdf.getConstant('NC_GLOBAL'),'classtype',klasstring); + if verbose > 4, + disp(sprintf("=}{=creating var for %s.%s",groups{i}, fields{j})); + end + + [DimSize,DimValue,varid]=CreateVar(ncid,Var,groupID,fields{j},DimSize,DimValue); + if ~isempty(varid), + FillVar(Var,groupID,varid); + end + + elseif isempty(fieldnames(md.(groups{i}).(fields{j}))) % this is an empty struct, jus treat it as normal + klass=class(md.(groups{i})); + klasstring = strcat(klass, '.',klass); + netcdf.putAtt(groupID,netcdf.getConstant('NC_GLOBAL'),'classtype',klasstring); + if verbose > 4, + disp(sprintf("=[]=creating var for %s.%s",groups{i}, fields{j})); + end + + [DimSize,DimValue,varid]=CreateVar(ncid,Var,groupID,fields{j},DimSize,DimValue); + if ~isempty(varid), + FillVar(Var,groupID,varid); + end + + else + 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})))) + end + elseif sum(strcmp(class(Var), issmclasses)) == 1, % that is an issm class + if strcmp(class(Var), 'solution'), + if verbose > 4, + disp(sprintf("=$$=creating var for %s.%s",groups{i}, fields{j})) + disp("NEED treatment") + end + elseif strcmp(class(Var), 'dict'), %we have potential for a dict in py not to sure what it translates to here. + if verbose > 4, + disp(sprintf("=WW=creating var for %s.%s",groups{i}, fields{j})) + disp("NEED Treatment") + end + + else + klass=class(md.(groups{i})); + klasstring = strcat(klass, '.',klass); + netcdf.putAtt(groupID,netcdf.getConstant('NC_GLOBAL'),'classtype',klasstring); + subgroupID=netcdf.defGrp(groupID,fields{j}); + klass=class(md.(groups{i}).(fields{j})); + klasstring = strcat(klass, '.',klass); + netcdf.putAtt(subgroupID,netcdf.getConstant('NC_GLOBAL'),'classtype',klasstring); + subfields=fieldnames(Var); + for k=1:length(subfields), + if sum(strcmp(subfields{k},["outlog" "errlog"])) == 0, + if verbose > 4, + disp(sprintf("+==+creating var for %s.%s.%s",groups{i}, fields{j}, subfields{k})) + end + Var=md.(groups{i}).(fields{j}).(subfields{k}); + [DimSize,DimValue,varid]=CreateVar(ncid,Var,subgroupID,subfields{k},DimSize,DimValue); + if ~isempty(varid), + FillVar(Var,subgroupID,varid); + end + end + end + + end + else + 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))) + end + end + end + netcdf.close(ncid); +end + +function [DimSize,DimValue,varid]=CreateVar(ncid,Var,groupID,field,DimSize,DimValue) +% Grab dimensions + varsize=size(Var); + varlength=length(Var); + % treating scalar string or bool as atribute + if isa(Var,'logical'), + if Var, + LogicString='True'; + else, + LogicString='False'; + end + netcdf.putAtt(groupID,netcdf.getConstant('NC_GLOBAL'),field,LogicString); + varid=[]; + + elseif isa(Var,'char'), + if strcmp(field,'name'), % it looks like netCDF does not like attributes that are called "name" + field = 'varname'; + end + if size(Var,1) <= 1 %that is a single string or empty + netcdf.putAtt(groupID,netcdf.getConstant('NC_GLOBAL'),field,Var); + varid=[]; + else % that is a character array + [dims,DimSize,DimValue]=GetDims(ncid,Var,DimSize,DimValue); + varid = netcdf.defVar(groupID,field,'NC_CHAR',dims); + if numel(Var)>1 + netcdf.defVarDeflate(groupID,varid,true,true,4); + end + end + + elseif isa(Var,'double'), %dealing with arrays + if all(mod(Var, 1) == 0, 'all') %those are actually integers, + [dims,DimSize,DimValue]=GetDims(ncid,Var,DimSize,DimValue); + varid = netcdf.defVar(groupID,field,'NC_INT64',dims); + if numel(Var)>1 + netcdf.defVarDeflate(groupID,varid,true,true,4); + end + else + [dims,DimSize,DimValue]=GetDims(ncid,Var,DimSize,DimValue); + varid = netcdf.defVar(groupID,field,'NC_DOUBLE',dims); + if numel(Var)>1 + netcdf.defVarDeflate(groupID,varid,true,true,4); + end + end + elseif isa(Var,'cell'), + % cells can be a range of things, what are we dealing with here + if isempty(Var), + netcdf.putAtt(groupID,netcdf.getConstant('NC_GLOBAL'),field,'emptycell'); + varid=[]; + else + [dims,DimSize,DimValue]=GetDims(ncid,Var,DimSize,DimValue); + if isa(Var{1}, 'double'), + varid = netcdf.defVar(groupID,field,'NC_DOUBLE',dims); + if numel(Var)>1 + netcdf.defVarDeflate(groupID,varid,true,true,4); + end + else + varid = netcdf.defVar(groupID,field,'NC_CHAR',dims); + if numel(Var)>1 + netcdf.defVarDeflate(groupID,varid,true,true,4); + end + end + end + elseif isa(Var,'struct'), + if isempty(fieldnames(Var)), + netcdf.putAtt(groupID,netcdf.getConstant('NC_GLOBAL'),field,'emptystruct'); + varid=[]; + else + %Start by getting the structure fields and size + [dims,DimSize,DimValue]=GetDims(ncid,Var,DimSize,DimValue); + varid = netcdf.defVar(groupID,field,'NC_CHAR',dims); + if numel(Var)>1 + netcdf.defVarDeflate(groupID,varid,true,true,4); + end + end + else + disp(sprintf('no support for class %s of field %s',class(Var),field)); + varid=[]; + end + return +end + + +function FillVar(Var,groupID,varid) +% Grab dimensions + varsize=size(Var); + varlength=length(Var); + % treating scalar string or bool as atribute + if isa(Var,'double'), %dealing with arrays + if all(mod(Var, 1) == 0, 'all') %those are actually integers, + Var = int64(Var); + end + if length(Var)==0, + netcdf.putVar(groupID,varid,NaN); + else + netcdf.putVar(groupID,varid,Var); + end + elseif isa(Var,'char'), % at this point this should be a character array + netcdf.putVar(groupID,varid,Var); + elseif isa(Var,'cell'), % there can be a number of things in a cell array + for i=1:length(Var), + if isa(Var{i},'char') %for characters we limit the size to 40 for now + if length(Var)>1, + count=[min(length(Var{i}),40), 1]; + startpoint=[0 i-1]; + else + count=min(length(Var{i}),40); + startpoint=0; + end + + if length(Var{i})>40, + netcdf.putVar(groupID,varid,startpoint,count,Var{i}(1:40)); + disp(sprintf('some variable have been truncated')); + else + netcdf.putVar(groupID,varid,startpoint,count,Var{i}); + end + elseif isa(Var{i},'double') + startpoint=[i-1]; + count=[1 length(Var{i}) ndims(Var{i})]; + for j=1:ndims(Var{i}), + startpoint=[startpoint 0]; + end + netcdf.putVar(groupID,varid,startpoint,count,Var{i}); + else + disp(sprintf("WARNING: cell of class %s is not supported.",class(Var{i}))) + end + end + elseif isa(Var,'struct'), + %Start by getting the structure fields and size + locfields=fieldnames(Var); + for i=1:length(locfields), + for j=1:2, + if j==1, + CharVar=locfields{i}'; + disp(size(CharVar)) + if length(CharVar)==0 + CharVar='emptystruct'; + end + startpoint=[0,0,i-1] + else + if isa(Var.(locfields{i}),'char'), + CharVar=Var.(locfields{i})'; + else + CharVar=num2str(Var.(locfields{i}))'; + end + if length(CharVar)==0 + CharVar='emptystruct'; + end + startpoint=[0,1,i-1] + end + + extent=[min(length(CharVar),40), 1, 1] + if length(CharVar)>40, + netcdf.putVar(groupID,varid,startpoint,extent,CharVar(1:40)); + disp(sprintf('some variable have been truncated')); + else + netcdf.putVar(groupID,varid,startpoint,extent,CharVar); + end + end + end + else + disp(sprintf('no support for class %s',class(Var))); + end + return +end + +function [dims,DimSize,DimValue]=GetDims(ncid,Var,DimSize,DimValue) + dims=[]; + celldims=[]; + dim=ndims(Var); + if isa(Var,'struct'), + varsize=length(fieldnames(Var)); + else + varsize=size(Var); + if isa(Var, 'cell') + %we add the dimension of the cells themselves, + %that will most probably fail if cells have different sizes + for i=1:dim, + newdim=size(Var{i}); + if ~ismember(newdim, celldims), + celldims=[celldims newdim]; + end + end + end + end + varsize=[varsize celldims]; + alldim=length(varsize); + if dim>0, + for i=1:alldim, + if size(Var, i)>1 || i>dim || isa(Var, 'struct'), %we skip dimensions with zero lenght but want to add dimensions from cells + indsize=find(varsize(i)==DimValue); + if length(indsize)>0 + dims=[dims DimSize(indsize).index]; + else + indsize=length(DimSize)+1; + DimSize(indsize).index=netcdf.defDim(ncid,['DimNum' num2str(indsize)],varsize(i)); + [DimSize(indsize).name,DimSize(indsize).value]=netcdf.inqDim(ncid,DimSize(indsize).index); + DimValue(indsize)=DimSize(indsize).value; + dims=[dims DimSize(indsize).index]; + end + end + end + end + if isa(Var, 'cell') && isa(Var{1}, 'char'), + %if we have an cell variable with strings we need to add a stringlength + dims=[dims DimSize(4).index]; + end + % struct also need an extra dimension 2, but only if non empty + if isa(Var,'struct'), + dims=[DimSize(4).index DimSize(3).index, dims]; + end +end Index: /issm/trunk-jpl/src/m/contrib/defleurian/netCDF/export_netCDF.py =================================================================== --- /issm/trunk-jpl/src/m/contrib/defleurian/netCDF/export_netCDF.py (revision 27890) +++ /issm/trunk-jpl/src/m/contrib/defleurian/netCDF/export_netCDF.py (revision 27890) @@ -0,0 +1,536 @@ +from netCDF4 import Dataset +import numpy as np +import numpy.ma as ma +import time +import collections +from inspect import isclass +from os import path, remove + + +class ResTable: + def __init__(self): + self.data = [] + self.sizes = [] + + def update(self, stepvar): + #if we have a scalar we just add it to the end + #we save the size of the current step for further treatment + if len(np.shape(stepvar)) == 0: + self.sizes.append(1) + self.data.append(stepvar) + # if it is an array we add the values one by one + #we save the size of the current step for further treatment + else: + self.sizes.append([np.shape(stepvar)]) + stackdat = np.squeeze(stepvar.flatten()) + self.data = np.hstack((self.data, stackdat)) + + def finalize(self, rows): + #we have more scalars than steps, so we have an array + if len(self.data) > rows: + datasize = np.squeeze(self.sizes) + maxsize = [] + try: + dims = np.arange(np.shape(datasize)[1]) + for dim in dims: + maxsize.append(np.nanmax(datasize[:, dim])) + except IndexError: + if datasize.ndim == 0: + maxsize.append(datasize) + else: + maxsize.append(np.nanmax(datasize[:])) + findim = np.insert(maxsize, 0, rows) + #first check if all steps are the same size + if datasize.ndim == 0: + SameSize = True + else: + SameSize = np.sum(np.abs(datasize - datasize[0])) == 0 + if SameSize: + #same size for all steps, just reshape + return np.reshape(self.data, newshape=(findim)) + else: + #different sizes at each steps, first create a table big enough for the biggest step + startpoint = 0 + outdat = np.nan * np.ones(findim) + for step in range(rows): + #slicer is the data slice in the final array + slicer = [slice(0, d) for d in datasize[step, :]] + slicer = np.insert(slicer, 0, step) + curlen = int(np.prod(datasize[step, :])) + outdat[tuple(slicer)] = np.reshape(self.data[startpoint:startpoint + curlen], newshape=(datasize[step, :])) + startpoint += curlen + #outmasked = ma.masked_array(outdat, mask=np.where(np.isnan(outdat), 1, 0)) + return outdat + #as much scalars as steps (or less) so just one value per step + else: + return np.squeeze(np.asarray(self.data)) + + +def export_netCDF(md, filename): # {{{ + #verbosity of the code, 0 is no messages, 5 is chatty + verbose = 0 + if path.exists(filename): + print('File {} allready exist'.format(filename)) + newname = input('Give a new name or "delete" to replace: ') + if newname == 'delete': + remove(filename) + else: + print(('New file name is {}'.format(newname))) + filename = newname + #create file and define it + NCData = Dataset(filename, 'w', format='NETCDF4') + NCData.description = 'Results for run' + md.miscellaneous.name + NCData.history = 'Created ' + time.ctime(time.time()) + # define netCDF dimensions + #grab time from Transient if it exists + try: + StepNum = len(md.results.TransientSolution) + except AttributeError: #no transient so just one timestep + StepNum = 1 + except TypeError: #this isnot a result so no results in there + StepNum = 0 + TimeDim = NCData.createDimension('Time', StepNum) # time is first + DimDict = {len(TimeDim): 'Time'} + dimindex = 1 + UnlimDim = NCData.createDimension('Unlim', None) # unlimited dimension if needed + DimDict[len(UnlimDim)] = {'Inf'} + dimindex = 2 + #add mesh related dimension that we know are needed + dimlist = [2, 40, md.mesh.numberofelements, md.mesh.numberofvertices, np.shape(md.mesh.elements)[1]] + dimnames = ['DictDummy', 'StringLength', 'EltNum', 'VertNum', 'VertPerElt'] + try: + dimlist = dimlist + [md.mesh.numberofedges] + dimnames = dimnames + ['EdgeNum'] + except AttributeError: + #no edges on this mesh, we fix it at 0 + dimlist += [0] + dimnames += ['EdgeNum'] + if verbose > 0: + print('===Creating dimensions ===') + for i, newdim in enumerate(dimlist): + if newdim not in list(DimDict.keys()): + dimindex += 1 + NewDim = NCData.createDimension(dimnames[i], newdim) + DimDict[len(NewDim)] = dimnames[i] + + typelist = [bool, str, int, float, complex, + collections.OrderedDict, + np.int64, np.ndarray, np.float64] + + # get all model classes and create respective groups + groups = dict.keys(md.__dict__) + if verbose > 0: + print('===Creating and populating groups===') + for group in groups: + if verbose > 1: + print('===Now treating {}==='.format(group)) + if group == 'qmu': + print("qmu is skipped until it is more stable") + continue + NCgroup = NCData.createGroup(str(group)) + # In each group gather the fields of the class + try: + fields = dict.keys(md.__dict__[group].__dict__) + except AttributeError: + print("WARNING: md.{} as no fields, we skip it.".format(group)) + continue + # looping on fields in each group + for field in fields: + Var = md.__dict__[group].__dict__[field] + # Special treatment for list fields + if type(Var) == list: + StdList = False + if len(Var) == 0: + StdList = True #this is an empty list + else: + #returns False for exotic types (typicaly results) + StdList = type(Var[0]) in typelist + klass = type(md.__dict__[group]).__module__ + '.' + type(md.__dict__[group]).__name__ + NCgroup.__setattr__('classtype', klass) + if StdList: # this is a standard or empty list just proceed + if verbose > 4: + print("=££=creating var for {}.{} with classtype : {}".format(group, field, klass)) + Var = SqueezeVar(Var) + DimDict, ncvar = CreateVar(NCData, Var, field, NCgroup, DimDict) + if ncvar is not None: + FillVar(ncvar, Var) + else: # this is a list of fields, specific treatment needed (usually results or outputdefinitions) + if verbose > 4: + print("=??=we have a list of fields for {}.{} with classtype : {}".format(group, field, klass)) + Listsize = len(Var) + if group == 'results': #for results we reshape the datas following time rather than subgrouping + Subgroup = NCgroup.createGroup(str(field)) + try: + #take the class of the first element to define nc class and get the list of variables + klass = type(md.__dict__[group].__dict__[field][0]).__module__ + '.' + type(md.__dict__[group].__dict__[field][0]).__name__ + Subgroup.__setattr__('classtype', klass) + subfields = dict.keys(md.__dict__[group].__dict__[field][0].__dict__) + except (IndexError, AttributeError): + klass = type(md.__dict__[group].__dict__[field]).__module__ + '.' + type(md.__dict__[group].__dict__[field]).__name__ + Subgroup.__setattr__('classtype', klass) + subfields = dict.keys(md.__dict__[group].__dict__[field].__getitem__(0)) + for subfield in subfields: + if subfield not in ['errlog', 'outlog']: + StackedVar = ResTable() + #first loop over the field (result type) to find the index of the last subfield (variable) + for listindex in range(0, Listsize): + try: + Var = md.__dict__[group].__dict__[field].__getitem__(listindex).__dict__[subfield] + lastindex = listindex + 1 + except AttributeError: + Var = md.__dict__[group].__dict__[field].__getitem__(listindex)[subfield] + except KeyError: + #Some fields only exist for the first step + lastindex = listindex + continue + #Add the subfield at the current step + Var = SqueezeVar(Var) + StackedVar.update(Var) + if verbose > 4: + print("=@@=creating var for {}.{}.{}".format(group, field, subfield)) + print("last index of the list is {}".format(lastindex)) + StackedVar = SqueezeVar(StackedVar.finalize(int(lastindex))) + DimDict, ncvar = CreateVar(NCData, StackedVar, subfield, Subgroup, DimDict) + #and fill it up + if ncvar is not None: + FillVar(ncvar, StackedVar) + elif group == 'outputdefinition': #for outputdefinition we keep a subgroup format + for listindex in range(0, Listsize): + Subgroupname = str(md.__dict__[group].__dict__[field][listindex].definitionstring) + Subgroup = NCgroup.createGroup(Subgroupname) + klass = type(md.__dict__[group].__dict__[field][listindex]).__module__ + '.' + type(md.__dict__[group].__dict__[field][listindex]).__name__ + Subgroup.__setattr__('classtype', klass) + subfields = dict.keys(md.__dict__[group].__dict__[field][listindex].__dict__) + for subfield in subfields: + Var = md.__dict__[group].__dict__[field].__getitem__(listindex).__dict__[subfield] + Var = SqueezeVar(Var) + if verbose > 4: + print("=--=creating var for {}.{}[{}].{}".format(group, field, listindex, subfield)) + DimDict, ncvar = CreateVar(NCData, Var, subfield, Subgroup, DimDict) + #and fill it up + if ncvar is not None: + FillVar(ncvar, Var) + else: + print("WARNING: unknown treatment for md.{}".format(group)) + # No subgroup, we directly treat the variable + elif type(md.__dict__[group].__dict__[field]) in typelist or field == 'bamg': + klass = type(md.__dict__[group]).__module__ + '.' + type(md.__dict__[group]).__name__ + NCgroup.__setattr__('classtype', klass) + Var = md.__dict__[group].__dict__[field] + Var = SqueezeVar(Var) + if verbose > 4: + print("====creating var for {}.{}".format(group, field)) + DimDict, ncvar = CreateVar(NCData, Var, field, NCgroup, DimDict) + if ncvar is not None: + FillVar(ncvar, Var) + # empty field, do nothing + elif md.__dict__[group].__dict__[field] is None: + print('field md.{}.{} is None'.format(group, field)) + # if it is a masked array + elif type(md.__dict__[group].__dict__[field]) is np.ma.core.MaskedArray: + klass = type(md.__dict__[group]).__module__ + '.' + type(md.__dict__[group]).__name__ + NCgroup.__setattr__('classtype', klass) + Var = md.__dict__[group].__dict__[field].data + Var = SqueezeVar(Var) + if verbose > 4: + print("=++=creating var for {}.{}".format(group, field)) + DimDict, ncvar = CreateVar(NCData, Var, field, NCgroup, DimDict) + if ncvar is not None: + FillVar(ncvar, Var) + # this is an issm class + elif isclass(type(md.__dict__[group].__dict__[field])): + if type(md.__dict__[group].__dict__[field]).__name__ == 'solution': + #for results we reshape the datas following time rather than subgrouping + Listsize = len(md.__dict__[group].__dict__[field]) + Subgroup = NCgroup.createGroup(str(field)) + try: + #take the class of the first element to define nc class and get the list of variables + klass = type(md.__dict__[group].__dict__[field][0]).__module__ + '.' + type(md.__dict__[group].__dict__[field][0]).__name__ + Subgroup.__setattr__('classtype', klass) + subfields = dict.keys(md.__dict__[group].__dict__[field][0].__dict__) + except (IndexError, AttributeError): + klass = type(md.__dict__[group].__dict__[field]).__module__ + '.' + type(md.__dict__[group].__dict__[field]).__name__ + Subgroup.__setattr__('classtype', klass) + subfields = dict.keys(md.__dict__[group].__dict__[field].__getitem__(0)) + for subfield in subfields: + if subfield not in ['errlog', 'outlog']: + StackedVar = ResTable() + for listindex in range(0, Listsize): + try: + Var = md.__dict__[group].__dict__[field].__getitem__(listindex).__dict__[subfield] + lastindex = listindex + 1 + except AttributeError: + Var = md.__dict__[group].__dict__[field].__dict__[subfield] + lastindex = listindex + except KeyError: + #Some fields only exist for the first step + lastindex = listindex + break + Var = SqueezeVar(Var) + StackedVar.update(Var) + if verbose > 4: + print("=$$=creating var for {}.{}.{}".format(group, field, subfield)) + print("last index of the list is {}".format(lastindex)) + StackedVar = SqueezeVar(StackedVar.finalize(int(lastindex))) + DimDict, ncvar = CreateVar(NCData, StackedVar, subfield, Subgroup, DimDict) + #and fill it up + if ncvar is not None: + FillVar(ncvar, StackedVar) + elif type(md.__dict__[group].__dict__[field]).__name__ == 'dict': + # designed for a dict in dummy but might be used elsewhere + # there is no subgroup + klass = type(md.__dict__[group]).__module__ + '.' + type(md.__dict__[group]).__name__ + NCgroup.__setattr__('classtype', klass) + Var = md.__dict__[group].__dict__[field] + Var = SqueezeVar(Var) + if verbose > 4: + print("=WW=creating var for {}.{}".format(group, field)) + DimDict, ncvar = CreateVar(NCData, Var, field, NCgroup, DimDict) + if ncvar is not None: + FillVar(ncvar, Var) + else: + klass = type(md.__dict__[group]).__module__ + '.' + type(md.__dict__[group]).__name__ + NCgroup.__setattr__('classtype', klass) + Subgroup = NCgroup.createGroup(str(field)) + klass = type(md.__dict__[group].__dict__[field]).__module__ + '.' + type(md.__dict__[group].__dict__[field]).__name__ + Subgroup.__setattr__('classtype', klass) + subfields = dict.keys(md.__dict__[group].__dict__[field].__dict__) + for subfield in subfields: + if str(subfield) not in ['errlog', 'outlog']: + Var = md.__dict__[group].__dict__[field].__dict__[subfield] + Var = SqueezeVar(Var) + if verbose > 4: + print("+==+creating var for {}.{}.{}".format(group, field, subfield)) + DimDict, ncvar = CreateVar(NCData, Var, subfield, Subgroup, DimDict) + if ncvar is not None: + FillVar(ncvar, Var) + else: + print("WARNING, md.{}.{} is not treated as it does not fall in one of the existing cases.".format(group, field)) + + NCData.close() + +# }}} + + +def CreateVar(NCData, var, field, Group, DimDict, *SupDim): # {{{ + #================================================================= + # Define the variables + #================================================================= + # grab type + try: + val_type = str(var.dtype) + if val_type.startswith(' 0 + except TypeError: + #check if we are integer and under C long overflow + IsInt = np.mod(var, 1) == 0 and abs(var) < 2147483647 + if IsInt: + val_type = 'int64' + if val_shape in [(), (0,), 0] and not SupDim: + ncvar = Group.createVariable(str(field), TypeDict[val_type], zlib=True) + else: + dimensions, DimDict = GetDim(NCData, val_shape, val_type, DimDict, val_dim) + if SupDim: + dimensions = SupDim + dimensions + ncvar = Group.createVariable(str(field), TypeDict[val_type], dimensions=dimensions, zlib=True) + elif val_type in [int, 'int64']: + if val_shape in [(), (0,), 0] and not SupDim: + ncvar = Group.createVariable(str(field), TypeDict[val_type], zlib=True) + else: + dimensions, DimDict = GetDim(NCData, val_shape, val_type, DimDict, val_dim) + if SupDim: + dimensions = SupDim + dimensions + ncvar = Group.createVariable(str(field), TypeDict[val_type], dimensions=dimensions, zlib=True) + else: + print(('WARNING type "{}" is unknown for "{}.{}"'.format(val_type, Group.name, field))) + ncvar = None + return DimDict, ncvar +# }}} + + +def FillVar(ncvar, invar, *UnlimIndex): # {{{ + #================================================================= + # Define the variables + #================================================================= + # grab type + try: + val_type = str(invar.dtype) + if val_type.startswith(' 0: + index = len(DimDict) + 1 + NewDim = NCData.createDimension('DimNum' + str(index), (val_shape[dim])) + DimDict[len(NewDim)] = 'DimNum' + str(index) + output = output + [str(DimDict[val_shape[dim]])] + return tuple(output), DimDict +# }}} + + +def SqueezeVar(Var): # {{{ + vardim = len(np.shape(Var)) + if vardim > 1: + Var = np.squeeze(Var) + + return Var +# }}} + + +def grow(self, row): # {{{ + np.append(self.data, row) +# }}} Index: /issm/trunk-jpl/src/m/contrib/defleurian/netCDF/read_netCDF.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/defleurian/netCDF/read_netCDF.m (revision 27890) +++ /issm/trunk-jpl/src/m/contrib/defleurian/netCDF/read_netCDF.m (revision 27890) @@ -0,0 +1,110 @@ +function self=read_netCDF(filename) + + % Different types in the netcdf standard are: + % 2 for char + % 4 for integer + % 6 for doubles + + ncid=netcdf.open(filename,'NC_NOWRITE'); + groupIDs=netcdf.inqGrps(ncid);%retrieve group IDs + self=model; + %loop on groups + for i=1:length(groupIDs) + whichclass = netcdf.getAtt(groupIDs(i),netcdf.getConstant('NC_GLOBAL'),'classtype'); + groupName = netcdf.inqGrpName(groupIDs(i)); + %results needs a special treatment as it is a structure + if strcmp(whichclass,'results'), + subgroupIDs=netcdf.inqGrps(groupIDs(i));%retrieve group IDs + %define the model structure + self=setfield(self,groupName,struct); + for j=1:length(subgroupIDs) + subclass = netcdf.getAtt(subgroupIDs(j),netcdf.getConstant('NC_GLOBAL'),'classtype'); + self.results=setfield(self.results,subclass,struct); + [ndims nvar natts]=netcdf.inq(subgroupIDs(j)); + varIDs=netcdf.inqVarIDs(subgroupIDs(j)); + %first loop on group atributes + for k=1:natts, + attname = netcdf.inqAttName(subgroupIDs(j),netcdf.getConstant('NC_GLOBAL'),k-1); + [xtype,attlen] = netcdf.inqAtt(subgroupIDs(j),netcdf.getConstant('NC_GLOBAL'),attname); + disp(sprintf('In %s, Treating attribute %s of type %i',subclass,attname,xtype)); + %classtype have done is job, no need to keep it any more + if ~strcmp(attname,'classtype'), + attval=netcdf.getAtt(subgroupIDs(j),netcdf.getConstant('NC_GLOBAL'),attname); + if strcmp(attval,'False'), + self.(groupName).(subclass).(attname)=false; + elseif strcmp(attval,'True') + self.(groupName).(subclass).(attname)=true; + else + self.(groupName).(subclass).(attname)=attval; + end + end + end + %now loop on variable in group + count=0; + for k=1:length(varIDs), + [varname, xtype, varDimIDs, varAtts] =netcdf.inqVar(subgroupIDs(j),varIDs(k)); + disp(sprintf('In %s, Treating variable %s of type %i',whichclass,varname,xtype)); + %time dimension seems to be last in our construction + for l=1:length(varDimIDs), + [dimname, dimlen] = netcdf.inqDim(ncid,varDimIDs(l)); + count(l)=[dimlen]; + end + startpoint=zeros(size(varDimIDs)); + timestep=count(end); + count(end)=1; + for l=1:timestep, + data=netcdf.getVar(subgroupIDs(j),varIDs(k),startpoint,count); + self.(groupName).(subclass)(l).(varname)=data; + startpoint(end)=startpoint(end)+1; + self.(groupName).(subclass)(l).('errlog')=''; + self.(groupName).(subclass)(l).('outlog')=''; + self.(groupName).(subclass)(l).('SolutionType')=subclass; + end + count=0; + end + end + else, + %define the model structure + self.(groupName)=eval(whichclass); + varIDs=netcdf.inqVarIDs(groupIDs(i)); + [ndims nvar natts]=netcdf.inq(groupIDs(i)); + %first loop on group atributes + for j=1:natts, + attname = netcdf.inqAttName(groupIDs(i),netcdf.getConstant('NC_GLOBAL'),j-1); + [xtype,attlen] = netcdf.inqAtt(groupIDs(i),netcdf.getConstant('NC_GLOBAL'),attname); + disp(sprintf('In %s, Treating attribute %s of type %i',whichclass,attname,xtype)); + %classtype have done is job, no need to keep it any more + if ~strcmp(attname,'classtype'), + attval=netcdf.getAtt(groupIDs(i),netcdf.getConstant('NC_GLOBAL'),attname); + if strcmp(attval,'False'), + self.(groupName).(attname)=false; + elseif strcmp(attval,'True') + self.(groupName).(attname)=true; + else + self.(groupName).(attname)=attval; + end + end + end + %now loop on variable in group + for j=1:length(varIDs), + [varname, xtype, varDimIDs, varAtts] =netcdf.inqVar(groupIDs(i),varIDs(j)); + disp(sprintf('In %s, Treating variable %s of type %i',whichclass,varname,xtype)); + %if the value is a single string, we need to transpose it (cross check with python file is necessary) + if xtype==2 + varval=netcdf.getVar(groupIDs(i),varIDs(j))'; + varval=cellstr(varval)'; + if strcmp(varval{1},'emptystruct'), + self.(groupName).(varname)=struct; + elseif strcmp(varval{1},'emptycell'), + self.(groupName).(varname)=cell(0,0); + else + self.(groupName).(varname)=varval; + end + else + self.(groupName).(varname)=netcdf.getVar(groupIDs(i),varIDs(j)); + end + end + end + end + netcdf.close(ncid) +end Index: /issm/trunk-jpl/src/m/contrib/defleurian/netCDF/read_netCDF.py =================================================================== --- /issm/trunk-jpl/src/m/contrib/defleurian/netCDF/read_netCDF.py (revision 27890) +++ /issm/trunk-jpl/src/m/contrib/defleurian/netCDF/read_netCDF.py (revision 27890) @@ -0,0 +1,22 @@ +from netCDF4 import Dataset +from os import path + + +def netCDFRead(filename): + def walktree(data): + keys = list(data.groups.keys()) + yield keys + for key in keys: + for children in walktree(data.groups[str(key)]): + yield children + + if path.exists(filename): + print(('Opening {} for reading '.format(filename))) + NCData = Dataset(filename, 'r') + class_dict = {} + + for children in walktree(NCData): + for child in children: + class_dict[str(child)] = str(getattr(NCData.groups[str(child)], 'classtype') + '()') + + print(class_dict) Index: /issm/trunk-jpl/src/m/contrib/defleurian/netCDF/restable.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/defleurian/netCDF/restable.m (revision 27890) +++ /issm/trunk-jpl/src/m/contrib/defleurian/netCDF/restable.m (revision 27890) @@ -0,0 +1,60 @@ +classdef restable + properties(SetAccess=public) + data=[]; + sizes=[]; + end + methods + function self = update(self, stepvar) + if length(stepvar) == 1 + %if we have a scalar we just add it to the end + %we save the size of the current step for further treatment + self.sizes=[self.sizes;1]; + self.data=[self.data;stepvar]; + else + % if it is an array we add the values one by one + %we save the size of the current step for further treatment + self.sizes=[self.sizes;fliplr(size(stepvar))]; + %we need to transpose to follow the indexing + flatdata=reshape(stepvar', 1, []); + self.data=[self.data,flatdata]; + end + end + function outdat = finalize(self, rows) + if length(self.data)>rows, + if size(self.sizes, 1)==1, + %just one step, data don't need treatment + outdat=self.data; + else, + %we have more scalars than steps, so we have an array + maxsize=[]; + for d=1:size(self.sizes,2) + maxsize=[maxsize,max(self.sizes(:,d))]; + end + findim=[maxsize, rows]; + %first check if all steps are the same size + SameSize = sum(self.sizes - self.sizes(1, :))==0; + if SameSize, + %same size for all steps, just reshape + outdat=reshape(self.data, findim); + else, + %different sizes at each steps, first create a table big enough for the biggest step + startpoint=1; + datadim=ndims(self.data); + outdat=nan(findim); + for r=1:rows + curlen = prod(self.sizes(r, :)); + outdat(1:self.sizes(r,1), 1:self.sizes(r,2), r) = reshape(self.data(startpoint:startpoint+ ... + curlen-1),self.sizes(r,:)); + startpoint = startpoint+curlen; + end + + end + end, + else, + %as much scalars as steps (or less) so just one value per step + outdat=self.data; + + end + end + end +end