matlab_class/gcmfaces_IO/write2nctiles.m

function [dimOut]=write2nctiles(fileOut,fldIn,doCreate,varargin);
% WRITE2NCTILES writes a gcmfaces object to nctiles, tiled netcdf, files
%
% inputs : 
%    - fileOut is the name of the file to be created
%    - fldIn is the array to write to disk
%    - doCreate: if 1 (default) then create file ; if 0 then append to file. 
%
% additional paramaters can be provided afterwards in the {'name',value} form:
%         'fldName' is the nc variable name for fld (see "notes" regarding default)
%         'longName' is the corresponding long name ('' by default).
%         'units' is the unit of fld ('(unknown)' by default).
%         'missval' is the missing value (NaN by default).
%         'fillval' is the fill value (NaN by default).
%         'tileNo' is a map of tile indices (face # by default)
%         'coord' is auxilliary coordinates attribute (e.g. 'lon lat dep')
%         'dimsize' is the array size associated with 'coord'
%         'xtype' is the variable ('double' by default)
%         'rdm' is the extended estimate description ('' by default).
%         'descr' is the file description ('' by default).
%
% notes: 
%    - if fldName is not explicitly specified then the input variable 
%      name (if ok) or file name (otherwise) is used as fldName.
%    - netcdf dimensions are simply set to 'i1,i2,...'
%
% examples:
%
%    write2nctiles('test1',mygrid.RAC,1,{'fldName','RAC'});
%    RAC=read_nctiles('test1','RAC');
%
%    write2nctiles('test2',mygrid.RAC);
%    RAC=read_nctiles('test2','test2');
%    RAC=read_nctiles('test2');

gcmfaces_global;
if ~(myenv.useNativeMatlabNetcdf);
  error('only native matlab nectdf is supported in write2nctiles');
end;

doCheck=0;%set to one to print stuff to screen

if isempty(whos('fileOut')); error('fileOut must be specified'); end;
if isempty(whos('fldIn')); error('fldIn must be specified'); end;
if isempty(whos('doCreate')); doCreate=1; end;

fldInIsaGcmfaces=isa(fldIn,'gcmfaces');
if fldInIsaGcmfaces;
  fldInIsEmpty=isempty(fldIn{1});
else;
  fldInIsEmpty=isempty(fldIn);
end;
    
%set more optional paramaters to default values
fldName=inputname(2); 
if isempty(fldName); 
    [tmp1,fldName,tmp2] = fileparts(fileOut); 
end;
longName=''; units='(unknown)'; missval=NaN; fillval=NaN; dimIn=[];
tileNo=mygrid.XC; for ff=1:mygrid.nFaces; tileNo{ff}(:)=ff; end;
coord=''; dimsize=[]; xtype='double'; descr=''; rdm='';

%set more optional paramaters to user defined values
for ii=1:nargin-3;
    if ~iscell(varargin{ii});
        warning('inputCheck:write2nctiles_1',...
            ['write2nctiles expects \n'...
            '         its optional parameters as cell arrays. \n'...
            '         Argument no. ' num2str(ii+1) ' was ignored \n'...
            '         Type ''help write2nctiles'' for details.']);
    elseif ~ischar(varargin{ii}{1});
        warning('inputCheck:write2nctiles_2',...
            ['write2nctiles expects \n'...
            '         its optional parameters cell arrays \n'...
            '         to start with character string. \n'...
            '         Argument no. ' num2str(ii+1) ' was ignored \n'...
            '         Type ''help write2nctiles'' for details.']);
    else;
        if strcmp(varargin{ii}{1},'descr')|...
                strcmp(varargin{ii}{1},'rdm')|...
                strcmp(varargin{ii}{1},'fldName')|...
                strcmp(varargin{ii}{1},'longName')|...
                strcmp(varargin{ii}{1},'units')|...
                strcmp(varargin{ii}{1},'missval')|...
                strcmp(varargin{ii}{1},'fillval')|...
                strcmp(varargin{ii}{1},'tileNo')|...
                strcmp(varargin{ii}{1},'coord')|...
                strcmp(varargin{ii}{1},'dimsize')|...
                strcmp(varargin{ii}{1},'xtype')|...
                strcmp(varargin{ii}{1},'dimIn');
            eval([varargin{ii}{1} '=varargin{ii}{2};']);
        else;
            warning('inputCheck:write2nctiles_3',...
                ['unknown option ''' varargin{ii}{1} ''' was ignored']);
        end;
    end;
end;

%split fldIn (if isa gcmfaces) into tiles
tileList=unique(convert2vector(tileNo)); 
tileList=tileList(~isnan(tileList));
ntile=length(tileList);
%
if fldInIsaGcmfaces&~fldInIsEmpty;
 for ff=1:ntile; 
   tmp1=[];
   for gg=1:mygrid.nFaces;
     [tmpi,tmpj]=find(tileNo{gg}==ff);
     if ~isempty(tmpi);
       tmpi=[min(tmpi(:)):max(tmpi(:))];
       tmpj=[min(tmpj(:)):max(tmpj(:))];
       tmp1=fldIn{gg}(tmpi,tmpj,:,:);
     end;
   end;
   fldTiles{ff}=tmp1;
 end;
 clear fldIn;
elseif fldInIsaGcmfaces;
 for ff=1:ntile;
   fldTiles{ff}=[];
 end;
end;

%start processing loop
for ff=1:ntile;

if fldInIsaGcmfaces;
  fldTile=fldTiles{ff};
  %reverse order of dimensions
  nn=length(size(fldTile));
  fldTile=permute(fldTile,[nn:-1:1]);
  if ntile==1; clear fldTiles; end;
else;
  fldTile=fldIn;
end;
fileTile=[fileOut sprintf('.%04d.nc',ff)];

%select dimensions of relevance:
if ~fldInIsEmpty;
  nDim=length(size(fldTile));
  dimsize=size(fldTile);
elseif ~isempty(dimsize);
  nDim=length(dimsize);
else;
  error('undertermined array size');
end;

for iDim=1:nDim;
if dimsize(iDim)~=1;
    dimlist{iDim}=['i' num2str(iDim)];
    dimName{iDim}=['array index ' num2str(iDim)];
    eval(['dimvec.i' num2str(iDim) '=[1:dimsize(iDim)];']);
end;
end;

%omit singleton dimensions:
ii=find(dimsize~=1);
dimsize=dimsize(ii);
dimlist={dimlist{ii}};
dimName={dimName{ii}};

%check : 
if doCheck;
whos fldTile
descr
fldName
longName
units
missval
fillval
dimlist
dimName
dimsize
dimvec
keyboard;
end;

%output dimension information
dimOut{ff}=dimlist;

if doCreate;
  %create netcdf file:
  %-------------------
  if prod(dimsize)*4/1e9<1.5;%use (always available) basic netcdf:
    mode='clobber';
  else;%to allow for large file:
    mode='NETCDF4';
  end;
  ncid=nccreateFile(fileTile,mode);
  nc_global=netcdf.getConstant('NC_GLOBAL');

  if ~isempty(rdm); 
    descr2=[descr ' -- ' rdm{1}];
  else;
    descr2=descr;
  end;
  ncputAtt(ncid,'','description',descr2);
  for pp=2:length(rdm);
  tmp1=char(pp+63);
  netcdf.putAtt(ncid,nc_global,tmp1,rdm{pp});
  end;
  %append readme
  if length(rdm)>0; pp=length(rdm)+1; tmp1=char(pp+63); else; tmp1='A'; end;
  netcdf.putAtt(ncid,nc_global,tmp1,'file created using gcmfaces_IO/write2nctiles.m');
  ncputAtt(ncid,'','date',date);
  netcdf.putAtt(ncid,nc_global,'Conventions','CF-1.6')   

  ncputAtt(ncid,'','_FillValue',fillval);
  ncputAtt(ncid,'','missing_value',missval);

  %ncdefDim(ncid,'itxt',30);
  for dd=1:length(dimlist); ncdefDim(ncid,dimlist{dd},dimsize(dd)); end;

  for dd=1:length(dimlist);
    ncdefVar(ncid,dimlist{dd},'double',{dimlist{dd}});
    ncputAtt(ncid,dimlist{dd},'long_name',dimName{dd});
    ncputAtt(ncid,dimlist{dd},'units','1');
  end;
  ncclose(ncid);

  %fill in the dimensions dimensions values vectors:
  %-------------------------------------------------
  ncid=ncopen(fileTile,'write');
  for dd=1:length(dimlist);
    ncputvar(ncid,dimlist{dd},getfield(dimvec,dimlist{dd}));
  end;
  ncclose(ncid);
end;

%use dimentsion specified by user
if ~isempty(dimIn); dimlist=dimIn{ff}; end;

%define and fill field:
%----------------------
ncid=ncopen(fileTile,'write');
%
netcdf.reDef(ncid);
ncdefVar(ncid,fldName,xtype,flipdim(dimlist,2));%note the direction flip
if ~isempty(longName); ncputAtt(ncid,fldName,'long_name',longName); end;
if ~isempty(units); ncputAtt(ncid,fldName,'units',units); end;
if ~isempty(coord); ncputAtt(ncid,fldName,'coordinates',coord); end; 
if strcmp(fldName,'lon'); ncputAtt(ncid,fldName,'standard_name','longitude'); end;
if strcmp(fldName,'lat'); ncputAtt(ncid,fldName,'standard_name','latitude'); end;
if strcmp(fldName,'dep'); ncputAtt(ncid,fldName,'standard_name','depth'); end;
if strcmp(fldName,'tim'); ncputAtt(ncid,fldName,'standard_name','time'); end;
if strcmp(fldName,'land'); ncputAtt(ncid,fldName,'standard_name','land_binary_mask'); end;
if strcmp(fldName,'area'); ncputAtt(ncid,fldName,'standard_name','cell_area'); end;
if strcmp(fldName,'thic'); ncputAtt(ncid,fldName,'standard_name','cell_thickness'); end;
netcdf.endDef(ncid);
%
if ~fldInIsEmpty; ncputvar(ncid,fldName,fldTile); end;
%
ncclose(ncid);

end;%for ff=1:ntile;

1	gforget	1.1	function [dimOut]=write2nctiles(fileOut,fldIn,doCreate,varargin);
2	gforget	1.10	% WRITE2NCTILES writes a gcmfaces object to nctiles, tiled netcdf, files
3			%
4			% inputs :
5			% - fileOut is the name of the file to be created
6			% - fldIn is the array to write to disk
7			% - doCreate: if 1 (default) then create file ; if 0 then append to file.
8			%
9			% additional paramaters can be provided afterwards in the {'name',value} form:
10			% 'fldName' is the nc variable name for fld (see "notes" regarding default)
11			% 'longName' is the corresponding long name ('' by default).
12			% 'units' is the unit of fld ('(unknown)' by default).
13			% 'missval' is the missing value (NaN by default).
14			% 'fillval' is the fill value (NaN by default).
15			% 'tileNo' is a map of tile indices (face # by default)
16			% 'coord' is auxilliary coordinates attribute (e.g. 'lon lat dep')
17			% 'dimsize' is the array size associated with 'coord'
18			% 'xtype' is the variable ('double' by default)
19			% 'rdm' is the extended estimate description ('' by default).
20			% 'descr' is the file description ('' by default).
21			%
22			% notes:
23			% - if fldName is not explicitly specified then the input variable
24			% name (if ok) or file name (otherwise) is used as fldName.
25			% - netcdf dimensions are simply set to 'i1,i2,...'
26			%
27			% examples:
28			%
29			% write2nctiles('test1',mygrid.RAC,1,{'fldName','RAC'});
30			% RAC=read_nctiles('test1','RAC');
31			%
32			% write2nctiles('test2',mygrid.RAC);
33			% RAC=read_nctiles('test2','test2');
34			% RAC=read_nctiles('test2');
35	gforget	1.1
36			gcmfaces_global;
37	gforget	1.5	if ~(myenv.useNativeMatlabNetcdf);
38			error('only native matlab nectdf is supported in write2nctiles');
39			end;
40	gforget	1.1
41			doCheck=0;%set to one to print stuff to screen
42
43	gforget	1.10	if isempty(whos('fileOut')); error('fileOut must be specified'); end;
44			if isempty(whos('fldIn')); error('fldIn must be specified'); end;
45			if isempty(whos('doCreate')); doCreate=1; end;
46
47	gforget	1.7	fldInIsaGcmfaces=isa(fldIn,'gcmfaces');
48	gforget	1.8	if fldInIsaGcmfaces;
49			fldInIsEmpty=isempty(fldIn{1});
50			else;
51			fldInIsEmpty=isempty(fldIn);
52			end;
53	gforget	1.10
54	gforget	1.1	%set more optional paramaters to default values
55	gforget	1.10	fldName=inputname(2);
56			if isempty(fldName);
57			[tmp1,fldName,tmp2] = fileparts(fileOut);
58			end;
59			longName=''; units='(unknown)'; missval=NaN; fillval=NaN; dimIn=[];
60	gforget	1.4	tileNo=mygrid.XC; for ff=1:mygrid.nFaces; tileNo{ff}(:)=ff; end;
61	gforget	1.10	coord=''; dimsize=[]; xtype='double'; descr=''; rdm='';
62	gforget	1.1
63			%set more optional paramaters to user defined values
64			for ii=1:nargin-3;
65			if ~iscell(varargin{ii});
66			warning('inputCheck:write2nctiles_1',...
67			['write2nctiles expects \n'...
68			' its optional parameters as cell arrays. \n'...
69			' Argument no. ' num2str(ii+1) ' was ignored \n'...
70			' Type ''help write2nctiles'' for details.']);
71			elseif ~ischar(varargin{ii}{1});
72			warning('inputCheck:write2nctiles_2',...
73			['write2nctiles expects \n'...
74			' its optional parameters cell arrays \n'...
75			' to start with character string. \n'...
76			' Argument no. ' num2str(ii+1) ' was ignored \n'...
77			' Type ''help write2nctiles'' for details.']);
78			else;
79			if strcmp(varargin{ii}{1},'descr')\|...
80	gforget	1.2	strcmp(varargin{ii}{1},'rdm')\|...
81	gforget	1.1	strcmp(varargin{ii}{1},'fldName')\|...
82			strcmp(varargin{ii}{1},'longName')\|...
83			strcmp(varargin{ii}{1},'units')\|...
84			strcmp(varargin{ii}{1},'missval')\|...
85			strcmp(varargin{ii}{1},'fillval')\|...
86	gforget	1.4	strcmp(varargin{ii}{1},'tileNo')\|...
87	gforget	1.5	strcmp(varargin{ii}{1},'coord')\|...
88	gforget	1.8	strcmp(varargin{ii}{1},'dimsize')\|...
89			strcmp(varargin{ii}{1},'xtype')\|...
90	gforget	1.1	strcmp(varargin{ii}{1},'dimIn');
91			eval([varargin{ii}{1} '=varargin{ii}{2};']);
92			else;
93			warning('inputCheck:write2nctiles_3',...
94			['unknown option ''' varargin{ii}{1} ''' was ignored']);
95			end;
96			end;
97			end;
98
99	gforget	1.4	%split fldIn (if isa gcmfaces) into tiles
100			tileList=unique(convert2vector(tileNo));
101			tileList=tileList(~isnan(tileList));
102			ntile=length(tileList);
103			%
104	gforget	1.8	if fldInIsaGcmfaces&~fldInIsEmpty;
105	gforget	1.4	for ff=1:ntile;
106			tmp1=[];
107			for gg=1:mygrid.nFaces;
108			[tmpi,tmpj]=find(tileNo{gg}==ff);
109			if ~isempty(tmpi);
110			tmpi=[min(tmpi(:)):max(tmpi(:))];
111			tmpj=[min(tmpj(:)):max(tmpj(:))];
112			tmp1=fldIn{gg}(tmpi,tmpj,:,:);
113			end;
114			end;
115			fldTiles{ff}=tmp1;
116			end;
117	gforget	1.7	clear fldIn;
118	gforget	1.8	elseif fldInIsaGcmfaces;
119			for ff=1:ntile;
120			fldTiles{ff}=[];
121			end;
122	gforget	1.4	end;
123
124			%start processing loop
125			for ff=1:ntile;
126	gforget	1.1
127	gforget	1.7	if fldInIsaGcmfaces;
128	gforget	1.4	fldTile=fldTiles{ff};
129	gforget	1.1	%reverse order of dimensions
130			nn=length(size(fldTile));
131			fldTile=permute(fldTile,[nn:-1:1]);
132	gforget	1.7	if ntile==1; clear fldTiles; end;
133	gforget	1.1	else;
134			fldTile=fldIn;
135			end;
136	gforget	1.4	fileTile=[fileOut sprintf('.%04d.nc',ff)];
137	gforget	1.1
138			%select dimensions of relevance:
139	gforget	1.8	if ~fldInIsEmpty;
140			nDim=length(size(fldTile));
141			dimsize=size(fldTile);
142			elseif ~isempty(dimsize);
143			nDim=length(dimsize);
144			else;
145			error('undertermined array size');
146			end;
147
148	gforget	1.1	for iDim=1:nDim;
149	gforget	1.8	if dimsize(iDim)~=1;
150	gforget	1.1	dimlist{iDim}=['i' num2str(iDim)];
151			dimName{iDim}=['array index ' num2str(iDim)];
152	gforget	1.8	eval(['dimvec.i' num2str(iDim) '=[1:dimsize(iDim)];']);
153	gforget	1.1	end;
154			end;
155
156			%omit singleton dimensions:
157			ii=find(dimsize~=1);
158			dimsize=dimsize(ii);
159			dimlist={dimlist{ii}};
160			dimName={dimName{ii}};
161
162			%check :
163			if doCheck;
164			whos fldTile
165			descr
166			fldName
167			longName
168			units
169			missval
170			fillval
171			dimlist
172			dimName
173			dimsize
174			dimvec
175			keyboard;
176			end;
177
178			%output dimension information
179			dimOut{ff}=dimlist;
180
181			if doCreate;
182			%create netcdf file:
183			%-------------------
184	gforget	1.8	if prod(dimsize)*4/1e9<1.5;%use (always available) basic netcdf:
185	gforget	1.7	mode='clobber';
186			else;%to allow for large file:
187			mode='NETCDF4';
188			end;
189	gforget	1.9	ncid=nccreateFile(fileTile,mode);
190	gforget	1.6	nc_global=netcdf.getConstant('NC_GLOBAL');
191	gforget	1.1
192	gforget	1.3	if ~isempty(rdm);
193			descr2=[descr ' -- ' rdm{1}];
194			else;
195			descr2=descr;
196			end;
197	gforget	1.4	ncputAtt(ncid,'','description',descr2);
198	gforget	1.2	for pp=2:length(rdm);
199			tmp1=char(pp+63);
200			netcdf.putAtt(ncid,nc_global,tmp1,rdm{pp});
201			end;
202			%append readme
203	gforget	1.10	if length(rdm)>0; pp=length(rdm)+1; tmp1=char(pp+63); else; tmp1='A'; end;
204	gforget	1.2	netcdf.putAtt(ncid,nc_global,tmp1,'file created using gcmfaces_IO/write2nctiles.m');
205	gforget	1.1	ncputAtt(ncid,'','date',date);
206	gforget	1.5	netcdf.putAtt(ncid,nc_global,'Conventions','CF-1.6')
207
208			ncputAtt(ncid,'','_FillValue',fillval);
209			ncputAtt(ncid,'','missing_value',missval);
210	gforget	1.1
211	gforget	1.11	%ncdefDim(ncid,'itxt',30);
212	gforget	1.1	for dd=1:length(dimlist); ncdefDim(ncid,dimlist{dd},dimsize(dd)); end;
213
214			for dd=1:length(dimlist);
215			ncdefVar(ncid,dimlist{dd},'double',{dimlist{dd}});
216			ncputAtt(ncid,dimlist{dd},'long_name',dimName{dd});
217	gforget	1.5	ncputAtt(ncid,dimlist{dd},'units','1');
218	gforget	1.1	end;
219			ncclose(ncid);
220
221			%fill in the dimensions dimensions values vectors:
222			%-------------------------------------------------
223			ncid=ncopen(fileTile,'write');
224			for dd=1:length(dimlist);
225			ncputvar(ncid,dimlist{dd},getfield(dimvec,dimlist{dd}));
226			end;
227			ncclose(ncid);
228			end;
229
230			%use dimentsion specified by user
231			if ~isempty(dimIn); dimlist=dimIn{ff}; end;
232
233			%define and fill field:
234			%----------------------
235			ncid=ncopen(fileTile,'write');
236	gforget	1.5	%
237			netcdf.reDef(ncid);
238	gforget	1.8	ncdefVar(ncid,fldName,xtype,flipdim(dimlist,2));%note the direction flip
239	gforget	1.1	if ~isempty(longName); ncputAtt(ncid,fldName,'long_name',longName); end;
240			if ~isempty(units); ncputAtt(ncid,fldName,'units',units); end;
241	gforget	1.5	if ~isempty(coord); ncputAtt(ncid,fldName,'coordinates',coord); end;
242			if strcmp(fldName,'lon'); ncputAtt(ncid,fldName,'standard_name','longitude'); end;
243			if strcmp(fldName,'lat'); ncputAtt(ncid,fldName,'standard_name','latitude'); end;
244			if strcmp(fldName,'dep'); ncputAtt(ncid,fldName,'standard_name','depth'); end;
245			if strcmp(fldName,'tim'); ncputAtt(ncid,fldName,'standard_name','time'); end;
246			if strcmp(fldName,'land'); ncputAtt(ncid,fldName,'standard_name','land_binary_mask'); end;
247			if strcmp(fldName,'area'); ncputAtt(ncid,fldName,'standard_name','cell_area'); end;
248			if strcmp(fldName,'thic'); ncputAtt(ncid,fldName,'standard_name','cell_thickness'); end;
249			netcdf.endDef(ncid);
250			%
251	gforget	1.8	if ~fldInIsEmpty; ncputvar(ncid,fldName,fldTile); end;
252			%
253	gforget	1.1	ncclose(ncid);
254
255	gforget	1.4	end;%for ff=1:ntile;
256	gforget	1.1