Friday, January 7, 2011

Write netcdf file using matlab native functions

Check the help page from Mathworks for new version Matlab
Please leave your comments to help improve the code in the future. Your help is greatly appreciated!

function SaveNcVarMatlab(ncfile,ncvar,MagicTag,myValue,myDType)
% to save a variable to a netcdf file using native matlab netcdf-functions
% Usage: SaveNcVarMatlab(nc-filename,nc-varname,MagicTag,myValue,[myDType])
% e.g.,
% 1: to save a 4D data
%     SaveNcVarMatlab('mask.nc','tmask',{'t','z','y','x'},4D-data,'float');
%     if the nc-file exist, the code will try to append data to old
%     file as long as the dimensions of input data is consistent
%     with existing one. If necessary, new dimensions will be added
%     automatically.
% 2: to add an attribute to an variable in existing nc file
%     SaveNcVarMatlab('mask.nc','tmask','att','long_name','mask at t-point')
%     Will overwrite existing one automatically. If no value set to
%     the attribute, the code will try to delete that attribute.
% 3: add/delete a global attribute, similar to the second example
%     but won't use ncvar here.
% Note:
%    Overwriting an existing variable is possible but disabled in
%    this version. Also unlimited dimension is not considered
%    here.
% copyright @ http://scriptdemo.blogspot.com
% ref: matlab netcdf help

if nargin<4 help SaveNcVarMatlab; return; end

if iscell(MagicTag)
    if nargin~=5
        help SaveNcVarMatlab
        return
    end
elseif ischar(MagicTag)
    if ischar(myValue)
       tmpstr=myValue; clear myValue
       myValue.name=tmpstr; clear tmpstr
       if nargin==5
          myValue.value=myDType;
       else
           myValue.value='';
       end
    elseif isstruct(myValue)
       if ~isfield(myValue,'name')
           help SaveNcVarMatlab;
           return;
       else
           if ~isfield(myValue,'value')
              if nargin==5
                 myValue.value=myDType;
              else
                 myValue.value='';
              end
           end
       end
    elseif ~isstruct(myValue)
        disp('Not valid value for myValue.')
        help SaveNcVarMatlab
        return
    end
end

% to support my work, please keep the "creator" global attribute
mycreatorstr='SaveNcVarMatlab.m from http://scriptdemo.blogspot.com';

if iscell(MagicTag) % add values
  [ncfid,IsAppend,OldNC,IsOldVar,OldVar]=CheckNCExist(ncfile,ncvar);
  if length(MagicTag)==1
      myValue=reshape(myValue,[],1);
      vardiminput=length(myValue);
  else
     vardiminput=size(myValue);
  end

  if IsAppend
     for nd=1:length(MagicTag)
        if sum(ismember(OldNC.dimname,MagicTag{nd}))==0
           myDim.name{nd}=MagicTag{nd};
           myDim.len(nd)=vardiminput(nd);
           myDim.IsNew(nd)=1;
        else
           dim_loc=find(ismember(OldNC.dimname,MagicTag{nd})==1);
           if OldNC.dimlen(dim_loc)~=vardiminput(nd)
              disp(['The length of input dimension ',MagicTag{nd},' is not the same as the exist one'])
              return
              % should consider unlimite dimension in the future
           else
              myDim.name{nd}=MagicTag{nd};
              myDim.ID(nd)=OldNC.dimid(dim_loc);
              myDim.IsNew(nd)=0;
           end
        end
     end
  else
     for nd=1:length(MagicTag)
         myDim.name{nd}=MagicTag{nd};
         myDim.len(nd)=vardiminput(nd);
         myDim.IsNew(nd)=1;
     end
  end
  if IsOldVar
     % support overwrite???
     netcdf.abort(ncfid); return
  else
     if IsAppend netcdf.reDef(ncfid); end
     % define some new dimension
     for nd=1:length(myDim.name)
         if myDim.IsNew(nd)
            myDim.ID(nd)=netcdf.defDim(ncfid,myDim.name{nd},myDim.len(nd));
         end
     end
     myvarid=netcdf.defVar(ncfid,ncvar,myDType,fliplr(myDim.ID));
     netcdf.putAtt(ncfid,netcdf.getConstant('NC_GLOBAL'),'creator',mycreatorstr);
     netcdf.endDef(ncfid);
     %netcdf.putVar(ncid,varid,start,count,stride,data)
     myValue=permute(myValue,length(size(myValue)):-1:1);
     netcdf.putVar(ncfid,myvarid,myValue);
  end

elseif ischar(MagicTag)
  [ncfid,IsAppend,OldNC,IsOldVar,OldVar]=CheckNCExist(ncfile,ncvar);
  if IsAppend
     switch lower(MagicTag)
       case {'a','attribute','var_attribute','vatt','att'}
         if IsOldVar
            OldVar=GetVAtts(ncfid,OldVar);
            netcdf.reDef(ncfid);
            if isfield(OldVar,'attname')
               if ~isstruct(myValue) %delete attribute
                   if sum(ismember(OldVar.attname,myValue))==1
                      netcdf.delAtt(ncfid,OldVar.varid,myValue);
                   else
                       disp(['attribute: ',myValue,' not exist. Nothing is done'])
                       netcdf.abort(ncfid); return;
                   end
               else
                  if sum(ismember(OldVar.attname,myValue.name))==1
                     %netcdf.delAtt(ncfid,OldVar.varid,myValue.name);
                     if isempty(myValue.value)
                        netcdf.delAtt(ncfid,OldVar.varid,myValue.name);
                     else
                        % overwrite old one automatically
                        netcdf.putAtt(ncfid,OldVar.varid,myValue.name,myValue.value);
                     end
                  else
                     %new attribute
                     if ~isempty(myValue.value)
                       netcdf.putAtt(ncfid,OldVar.varid,myValue.name,myValue.value);
                     end
                  end
               end
            else
              %new attribute
              if ~isempty(myValue.value)
                 netcdf.putAtt(ncfid,OldVar.varid,myValue.name,myValue.value);
              end
            end
         else
            disp(['var: ',ncvar,' isn''t exist in ',ncfile]);
            netcdf.abort(ncfid);
            return
         end
       case {'ga','global_attribute','gatt'}
         OldNC=GetGAtts(ncfid,OldNC);
         netcdf.reDef(ncfid);
         if ~isstruct(myValue)
             if isfield(OldNC,'attname')
                if sum(ismember(OldNC.attname,myValue))==1
                   netcdf.delAtt(ncfid,netcdf.getConstant('GLOBAL'),myValue);
                else
                   disp(['attribute: ',myValue,' not exist. Nothing is done'])
                   netcdf.abort(ncfid); return;
                end
             else
                disp('No global attributes available.')
                netcdf.abort(ncfid); return;
             end
         else
            if isfield(OldNC,'attname')
               if sum(ismember(OldNC.attname,myValue.name))==1
                  netcdf.delAtt(ncfid,netcdf.getConstant('GLOBAL'),myValue.name);
                  if ~isempty(myValue.value)
                     netcdf.putAtt(ncfid,netcdf.getConstant('NC_GLOBAL'),myValue.name,myValue.value);
                  end
               else
                  %new global attribute
                  if ~isempty(myValue.value)
                     netcdf.putAtt(ncfid,netcdf.getConstant('NC_GLOBAL'),myValue.name,myValue.value);
                  end
               end
            else
               % also new global attribute
               if ~isempty(myValue.value)
                  netcdf.putAtt(ncfid,netcdf.getConstant('NC_GLOBAL'),myValue.name,myValue.value);
               else
                  netcdf.abort(ncfid); return;
               end
            end
         end
       otherwise
           disp(['Not defined attribute type: MagicTag: ',MagicTag])
           netcdf.abort(ncfid); return;
     end
  else
     netcdf.abort(ncfid);
     disp([ncfile,' dos not exist, global attribute can''t be added.'])
     return
  end
  netcdf.putAtt(ncfid,netcdf.getConstant('NC_GLOBAL'),'creator',mycreatorstr);
  netcdf.endDef(ncfid);
else
  disp('Not valid value for MagicTag')
  return;
end

%close the file
netcdf.close(ncfid);
end

function [ncfid,IsAppend,OldNC,IsOldVar,OldVar]=CheckNCExist(ncfile,ncvar)
  IsAppend=0; IsOldVar=0;
  if exist(ncfile,'file')
     disp([ncfile,' already exist, try the append mode']);
     IsAppend=1;
  end

  % Open the nc file
  if IsAppend
     ncfid=netcdf.open(ncfile,'NC_WRITE');
     [OldNC.numdims, OldNC.numvars, OldNC.numGAtt, unlimdimID] = netcdf.inq(ncfid);
     for NDim=1:OldNC.numdims
        [OldNC.dimname{NDim}, OldNC.dimlen(NDim)] = netcdf.inqDim(ncfid,NDim-1);
        OldNC.dimid(NDim)=NDim-1;
     end
     for NVar=1:OldNC.numvars
        [tmpvarname, tmpvarxtype,tmpvardims, tmpvarnumatts] = netcdf.inqVar(ncfid,NVar-1);
        if strcmp(tmpvarname,ncvar)
            IsOldVar=1;
            OldVar.varid=NVar-1;
            OldVar.varname=tmpvarname;
            OldVar.varxtype=tmpvarxtype;
            OldVar.vardimids=tmpvardims;
            OldVar.numatts=tmpvarnumatts;
            break;
        end
     end
     if IsOldVar==0 OldVar=nan; end
  else
     ncfid=netcdf.create(ncfile,'NC_NOCLOBBER');
     OldNC=nan; OldVar=nan;
  end
end

function OldVar=GetVAtts(ncfid,OldVar)
  for numA=1:OldVar.numatts
      OldVar.attname{numA}=netcdf.inqAttName(ncfid,OldVar.varid,numA-1);
      % OldVar.attid(numA)=netcdf.inqAttID(ncfid,varid,OldVar.attname{numA});
  end
end
function OldNC=GetGAtts(ncfid,OldNC)
  for numA=1:OldNC.numGAtt
      OldNC.attname{numA}= netcdf.inqAttName(ncfid,netcdf.getConstant('NC_GLOBAL'),numA-1);
  end
end

No comments:

ShowCalendar