[Dart-dev] [3963] DART/trunk/matlab: Removing dependency on CSIRO (getnc) and netcdf_toolbox.
nancy at ucar.edu
nancy at ucar.edu
Wed Jul 15 16:24:44 MDT 2009
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/dart-dev/attachments/20090715/151740ef/attachment-0001.html
-------------- next part --------------
Modified: DART/trunk/matlab/CheckModel.m
===================================================================
--- DART/trunk/matlab/CheckModel.m 2009-07-10 21:24:09 UTC (rev 3962)
+++ DART/trunk/matlab/CheckModel.m 2009-07-15 22:24:43 UTC (rev 3963)
@@ -8,7 +8,7 @@
% vars = CheckModel(fname)
% Data Assimilation Research Testbed -- DART
-% Copyright 2004-2007, Data Assimilation Research Section
+% Copyright 2004-2009, Data Assimilation Research Section
% University Corporation for Atmospheric Research
% Licensed under the GPL -- www.gpl.org/licenses/gpl.html
%
@@ -18,38 +18,26 @@
% $Revision$
% $Date$
-if ( exist(fname) ~= 2 )
- error(sprintf('%s does not exist.',fname))
-end
+if ( exist(fname,'file') ~= 2 ), error('%s does not exist.',fname); end
% Get some information from the file
-f = netcdf(fname);
-model = f.model(:);
+model = nc_attget(fname,nc_global,'model');
-atts = dim(f{'copy'}); num_copies = length(atts{1}); % determine # of ensemble members
-atts = dim(f{'time'}); num_times = length(atts{1}); % determine # of output times
+num_copies = dim_length(fname,'copy'); % determine # of ensemble members
+num_times = dim_length(fname,'time'); % determine # of output times
if (isempty(model))
- error(sprintf('%s has no ''model'' global attribute.',fname))
+ error('%s has no ''model'' global attribute.',fname)
end
-if (prod(size(num_copies)) > 1 )
- error(sprintf('%s has no ''copy'' dimension.',fname))
-end
-if (prod(size(num_times)) > 1 )
- error(sprintf('%s has no ''time'' dimension.',fname))
-end
-copy = getnc(fname,'copy');
+copy = nc_varget(fname,'copy');
switch lower(model)
case {'9var','lorenz_63','lorenz_84','ikeda'}
- atts = dim(f{'StateVariable'}); num_vars = length(atts{1}); % # of state varbls
- if (prod(size(num_vars)) > 1 )
- error(sprintf('%s has no ''StateVariable'' dimension.',fname))
- end
- StateVariable = getnc(fname,'StateVariable');
+ num_vars = dim_length(fname,'StateVariable'); % determine # of state varbls
+ StateVariable = nc_varget(fname,'StateVariable');
def_state_vars = zeros(1,num_vars); % for use as a subscript array,
def_state_vars(:) = StateVariable(:); % def_state_vars must be a row vector.
@@ -67,11 +55,8 @@
case {'lorenz_96', 'lorenz_04'}
- atts = dim(f{'StateVariable'}); num_vars = length(atts{1}); % # of state varbls
- if (prod(size(num_vars)) > 1 )
- error(sprintf('%s has no ''StateVariable'' dimension.',fname))
- end
- StateVariable = getnc(fname,'StateVariable');
+ num_vars = dim_length(fname,'StateVariable'); % determine # of state varbls
+ StateVariable = nc_varget(fname,'StateVariable');
% The only trick is to pick an equally-spaced subset of state
% variables for the default.
@@ -93,17 +78,14 @@
% This model has the state variables replicated, so there is a difference
% between num_state_vars and the length of the state variable.
- forcing = f.model_forcing(:);
- delta_t = f.model_delta_t(:);
- time_step_days = f.model_time_step_days(:);
- time_step_seconds = f.model_time_step_seconds(:);
- num_model_vars = cast(f.model_num_state_vars(:),'double'); % ACTUAL state vars, not
+ forcing = nc_attget(fname, nc_global, 'model_forcing');
+ delta_t = nc_attget(fname, nc_global, 'model_delta_t');
+ time_step_days = nc_attget(fname, nc_global, 'model_time_step_days');
+ time_step_seconds = nc_attget(fname, nc_global, 'model_time_step_seconds');
+ num_model_vars = nc_attget(fname, nc_global, 'model_num_state_vars');
- atts = dim(f{'StateVariable'}); num_vars = length(atts{1}); % # of state varbls
- if (prod(size(num_vars)) > 1 )
- error(sprintf('%s has no ''StateVariable'' dimension.',fname))
- end
- StateVariable = getnc(fname,'StateVariable');
+ num_vars = dim_length(fname,'StateVariable'); % determine # of state varbls
+ StateVariable = nc_varget(fname,'StateVariable');
% The only trick is to pick an equally-spaced subset of state
% variables for the default.
@@ -127,17 +109,11 @@
case 'lorenz_96_2scale'
- atts = dim(f{'Xdim'}); num_X = length(atts{1}); % # of X variables
- if (prod(size(num_X)) > 1 )
- error(sprintf('%s has no ''Xdim'' dimension.',fname))
- end
- Xdim = getnc(fname,'Xdim');
+ num_X = dim_length(fname,'Xdim'); % # of X variables
+ Xdim = nc_varget(fname,'Xdim');
- atts = dim(f{'Ydim'}); num_Y = length(atts{1}); % # of Y variables
- if (prod(size(num_Y)) > 1 )
- error(sprintf('%s has no ''Ydim'' dimension.',fname))
- end
- Ydim = getnc(fname,'Ydim');
+ num_Y = dim_length(fname,'Ydim'); % # of Y variables
+ Ydim = nc_varget(fname,'Ydim');
% The only trick is to pick an equally-spaced subset of state
% variables for the default.
@@ -158,16 +134,16 @@
case 'simple_advection'
- loc1d = getnc(fname,'loc1d');
- num_locs = length(loc1d);
+ num_locs = dim_length(fname,'loc1d'); % # of X variables
+ loc1d = nc_varget(fname,'loc1d');
- if ( isempty(f{'state'}) )
+ if ( nc_isvar(fname,'state') )
+ varnames = {'state'};
+ def_inds = [1 13 27];
+ else
varnames = {'concentration','source','wind', ...
'mean_source','source_phase'};
def_inds = round([1 , num_locs/3 , 2*num_locs/3]);
- else
- varnames = {'state'};
- def_inds = [1 13 27];
end
vars = struct('model' ,model, ...
@@ -191,19 +167,17 @@
% name(bob{1}) % is the variable name string
% bob{1}(:) % is the value of the netcdf variable (no offset/scale)
- varnames = {'ps','t','u','v'};
+ varnames = {'ps','t','u','v'};
num_vars = length(varnames);
- nlevels = ncsize(f('lev')); % determine # of state variables
- if (prod(size(nlevels)) > 1 )
- error(sprintf('%s has no ''lev'' dimension.',fname))
- end
-% times = getnc(fname,'time');
-% TmpI = getnc(fname,'TmpI'); % longitude
-% TmpJ = getnc(fname,'TmpJ'); % latitude
-% levels = getnc(fname,'level');
-% VelI = getnc(fname,'VelI'); % longitude
-% VelJ = getnc(fname,'VelJ'); % latitude
+ nlevels = dim_length(fname,'lev'); % determine # of state variables
+% times = nc_varget(fname,'time');
+% TmpI = nc_varget(fname,'TmpI'); % longitude
+% TmpJ = nc_varget(fname,'TmpJ'); % latitude
+% levels = nc_varget(fname,'level');
+% VelI = nc_varget(fname,'VelI'); % longitude
+% VelJ = nc_varget(fname,'VelJ'); % latitude
+
vars = struct('model',model, ...
'num_state_vars',num_vars, ...
'num_ens_members',num_copies, ...
@@ -223,10 +197,7 @@
varnames = {'PS','T','U','V','Q','CLDLIQ','CLDICE'};
num_vars = length(varnames);
- nlevels = ncsize(f('lev')); % determine # of state variables
- if (prod(size(nlevels)) > 1 )
- error(sprintf('%s has no ''lev'' dimension.',fname))
- end
+ nlevels = dim_length(fname,'lev'); % determine # of state variables
vars = struct('model',model, ...
'num_state_vars',num_vars, ...
@@ -244,15 +215,12 @@
% name(bob{1}) % is the variable name string
% bob{1}(:) % is the value of the netcdf variable (no offset/scale)
- num_vars = 22; % ps, t, u, v
- z_level = ncsize(f('z_level')); % determine # of state variables
- sl_level = ncsize(f('sl_level')); % determine # of state variables
- if (prod(size(z_level)) > 1 )
- error(sprintf('%s has no ''z_level'' dimension.',fname))
- end
- times = getnc(fname,'time');
- z_level = getnc(fname,'z_level');
- sl_level = getnc(fname,'sl_level');
+ num_vars = 22; % ps, t, u, v
+ z_level = dim_length(fname, 'z_level'); % determine # of state variables
+ sl_level = dim_length(fname,'sl_level'); % determine # of state variables
+ times = nc_varget(fname,'time');
+ z_level = nc_varget(fname,'z_level');
+ sl_level = nc_varget(fname,'sl_level');
vars = struct('model',model, ...
'num_state_vars',num_vars, ...
@@ -268,8 +236,8 @@
% so there is a separate function (GetPe2lyrInfo.m) that gets
% the information for each specific plot type.
- varnames = {'u','v','z'};
- num_vars = length(varnames);
+ varnames = {'u','v','z'};
+ num_vars = length(varnames);
vars = struct('model',model, ...
'num_state_vars',num_vars, ...
@@ -290,10 +258,7 @@
varnames = {'S','T','U','V','SSH'};
num_vars = length(varnames);
- nlevels = ncsize(f('ZG')); % determine # of state variables
- if (prod(size(nlevels)) > 1 )
- error(sprintf('%s has no ''ZG'' dimension.',fname))
- end
+ nlevels = dim_length(fname,'ZG'); % determine # of state variables
vars = struct('model',model, ...
'num_state_vars',num_vars, ...
@@ -306,8 +271,19 @@
otherwise
- error(sprintf('model %s unknown',model))
+ error('model %s unknown',model)
end
-close(f);
+
+
+
+function x = dim_length(fname,dimname)
+
+y = nc_isvar(fname,dimname);
+if (y < 1)
+ error('%s has no %s dimension/coordinate variable',fname,varname)
+end
+bob = nc_getdiminfo(fname,dimname);
+x = bob.Length;
+
Modified: DART/trunk/matlab/CheckModelCompatibility.m
===================================================================
--- DART/trunk/matlab/CheckModelCompatibility.m 2009-07-10 21:24:09 UTC (rev 3962)
+++ DART/trunk/matlab/CheckModelCompatibility.m 2009-07-15 22:24:43 UTC (rev 3963)
@@ -1,15 +1,15 @@
-function pinfo_out = CheckModelCompatibility(arg1, arg2);
+function pinfo_out = CheckModelCompatibility(arg1, arg2)
% CheckModelCompatibility tries to ensure that two netcdf files can be compared.
% There are 2 ways to call this: with 2 filenames, or with an already existing
% pinfo struct (with 2 filenames and 2 2-vector arrays for start/stop times).
% This routine fills in the 2-vectors with the time overlap region in a
% pinfo struct.
-% If the time indices are common between the 2 files it returns [1,nlength]
-% for both; if no overlap [-1,-1] for both; otherwise the [start,end] indices
-% for each array.
+% If the time indices are common between the 2 files it returns the
+% [start,count] indices for each array (indexing starts at 1,N).
+% It is an error situation if there is no overlap ([-1,-1] for both).
% Data Assimilation Research Testbed -- DART
-% Copyright 2004-2007, Data Assimilation Research Section
+% Copyright 2004-2009, Data Assimilation Research Section
% University Corporation for Atmospheric Research
% Licensed under the GPL -- www.gpl.org/licenses/gpl.html
%
@@ -35,77 +35,57 @@
error('Wrong number of arguments: must be 1 (pinfo) or 2 (file1,file2)')
end
+if ( exist(file1,'file') ~= 2 ), error('(file1) %s does not exist.',file1); end
+if ( exist(file2,'file') ~= 2 ), error('(file2) %s does not exist.',file2); end
+
% set this up for later
pinfo_out = setfield(pinfo_out, 'truth_time', [-1,-1]);
pinfo_out = setfield(pinfo_out, 'diagn_time', [-1,-1]);
-
-if ( exist(file1) ~= 2 )
- error(sprintf('(file1) %s does not exist.',file1))
-end
-if ( exist(file2) ~= 2 )
- error(sprintf('(file2) %s does not exist.',file2))
-end
-
% Get some information from the file1
-% dimensions are f1(xxxx), variables are f1{xxxx}
-f1 = netcdf(file1);
-tmodel = f1.model(:);
+tmodel = nc_attget(file1,nc_global,'model');
+
if (isempty(tmodel))
- error(sprintf('%s has no ''model'' global attribute.',file1))
+ error('%s has no ''model'' global attribute.',file1)
end
-if ( VarExist(f1,'copy') )
- tnum_copies = length(f1('copy')); % determine # of ensemble members
-else
- error(sprintf('%s has no ''copy'' dimension.',file1))
-end
-if ( VarExist(f1,'time') )
- tnum_times = length(f1('time')); % determine # of output times
-else
- error(sprintf('%s has no ''time'' dimension.',file1))
-end
-ttimes = f1{'time'}(:);
+tnum_copies = dim_length(file1,'copy');
+tnum_times = dim_length(file1,'time');
+ttimes = nc_varget(file1,'time');
-[tnum_vars,tdims] = ModelDimension(f1,tmodel);
+[tnum_vars,tdims] = ModelDimension(file1,tmodel);
if ( tnum_vars <= 0 )
- error(sprintf('Unable to determine resolution of %s.',file1))
+ error('Unable to determine resolution of %s.',file1)
end
-close(f1);
% Get some information from the file2
-f2 = netcdf(file2);
-dmodel = f2.model(:);
+dmodel = nc_attget(file1,nc_global,'model');
+
if (isempty(dmodel))
- error(sprintf('%s has no ''model'' global attribute.',file2))
+ error('%s has no ''model'' global attribute.',file2)
end
-if (VarExist(f2,'copy'))
- dnum_copies = length(f2('copy')); % determine # of ensemble members
-else
- error(sprintf('%s has no ''copy'' dimension.',file2))
-end
-if (VarExist(f2,'time'))
- dnum_times = length(f2('time')); % determine # of output times
-else
- error(sprintf('%s has no ''time'' dimension.',file2))
-end
-dtimes = f2{'time'}(:);
-[dnum_vars,ddims] = ModelDimension(f2,dmodel);
+
+dnum_copies = dim_length(file2,'copy');
+dnum_times = dim_length(file2,'time');
+dtimes = nc_varget(file2,'time');
+
+[dnum_vars,ddims] = ModelDimension(file2,dmodel);
if ( dnum_vars <= 0 )
- error(sprintf('Unable to determine resolution of %s.',file2))
+ error('Unable to determine resolution of %s.',file2)
end
-close(f2);
% rudimentary bulletproofing
if (strcmp(tmodel,dmodel) ~= 1)
- disp(sprintf('%s has model %s ',file1,tmodel))
- disp(sprintf('%s has model %s ',file2,dmodel))
+ fprintf('%s has model %s\n',file1,tmodel)
+ fprintf('%s has model %s\n',file2,dmodel)
error('no No NO ... models must be the same')
end
+pinfo_out = setfield(pinfo_out, 'model', tmodel);
+
if (prod(tnum_vars) ~= prod(dnum_vars))
- disp(sprintf('%s has %d state variables',file1,prod(tnum_vars)))
- disp(sprintf('%s has %d state variables',file2,prod(dnum_vars)))
+ fprintf('%s has %d state variables\n',file1,prod(tnum_vars))
+ fprintf('%s has %d state variables\n',file2,prod(dnum_vars))
error('no No NO ... both files must have same shape of state variables.')
end
@@ -123,10 +103,10 @@
( pinfo_out.truth_time(2) == -1 ) || ...
( pinfo_out.diagn_time(1) == -1 ) || ...
( pinfo_out.diagn_time(2) == -1 ))
- disp(sprintf('%s has %d timesteps, from %f to %f', ...
- file1,tnum_times,ttimes(1), ttimes(tnum_times)))
- disp(sprintf('%s has %d timesteps, from %f to %f', ...
- file2,dnum_times,dtimes(1), dtimes(dnum_times)))
+ fprintf('%s has %d timesteps, from %f to %f\n', ...
+ file1,tnum_times,ttimes(1), ttimes(tnum_times))
+ fprintf('%s has %d timesteps, from %f to %f\n', ...
+ file2,dnum_times,dtimes(1), dtimes(dnum_times))
error('These files have no timesteps in common')
end
@@ -148,9 +128,6 @@
pret.truth_time = [-1,-1];
pret.diagn_time = [-1,-1];
-%disp('at start')
-%pret
-
% ensure times are increasing and monotonic, and do they need to be
% a constant delta or not? compute delta array and validate those match?
% (to within an epsilon with floating pt roundoff)
@@ -159,16 +136,14 @@
% watch out for the floating point compares, and the min/max are probably
% redundant with the (1) and (l) comparisons, but until we put in checks
% for monotonicity, it's a cheap safety check.
-l = length(times1);
+len = length(times1);
if ( (length(times1) == length(times2)) ...
&& (abs(min(times1) - min(times2)) < epsilon) ...
&& (abs(max(times1) - max(times2)) < epsilon) ...
&& (times1(1) == times2(1)) ...
- && (times1(l) == times2(l)))
- pret.truth_time = [1,l];
- pret.diagn_time = [1,l];
-%disp('equal')
-%pret
+ && (times1(len) == times2(len)))
+ pret.truth_time = [1,len]; % start/count
+ pret.diagn_time = [1,len]; % start/count
return
end
@@ -191,9 +166,9 @@
minB = min(B);
maxA = max(A);
maxB = max(B);
-if (abs(minA - minB) < epsilon) , minB = minA; , end
-if (abs(maxA - minB) < epsilon) , minB = maxA; , end
-if (abs(maxA - maxB) < epsilon) , maxB = maxA; , end
+if (abs(minA - minB) < epsilon) , minB = minA; end
+if (abs(maxA - minB) < epsilon) , minB = maxA; end
+if (abs(maxA - maxB) < epsilon) , maxB = maxA; end
% case 1: disjoint regions; simply return here because
% return struct was initialized to the 'no intersection' case.
@@ -231,112 +206,71 @@
end
% now put the indices in the return struct and we are done.
-pret.truth_time = [min1,max1];
-pret.diagn_time = [min2,max2];
+pret.truth_time = [min1, max1-min1+1]; % start,count
+pret.diagn_time = [min2, min2-max2+1]; % start,count
% return here
-function x = VarExist(ncid,varname)
+function x = dim_length(fname,dimname)
-x = 0; % false ... assumed not to exist.
-variables = var(ncid);
-for i=1:length(variables)
- if ( strmatch(name(variables{i}), varname) == 1 )
- x = 1; % true ... variables exists in the netcdf file.
- end
+y = nc_isvar(fname,dimname);
+if (y < 1)
+ error('%s has no %s dimension/coordinate variable',fname,dimname)
end
+bob = nc_getdiminfo(fname,dimname);
+x = bob.Length;
-function x = DimExist(ncid,dimname)
-x = 0; % false ... assumed not to exist.
-dimensions = dim(ncid);
-for i=1:length(dimensions)
- if ( strmatch(name(dimensions{i}), dimname) == 1 )
- x = 1; % true ... variables exists in the netcdf file.
- end
-end
+function [x,y] = ModelDimension(fname,modelname)
-
-function [x,y] = ModelDimension(ncid,modelname)
-
x = 0;
y = NaN;
-%disp(sprintf('working with %s',modelname))
-
switch lower(modelname)
case 'cam'
- lonexist = VarExist(ncid,'lon');
- latexist = VarExist(ncid,'lat');
- lvlexist = VarExist(ncid,'lev');
- if ( latexist && lonexist && lvlexist )
- dnum_lons = prod(size(ncid('lon')));
- dnum_lats = prod(size(ncid('lat')));
- dnum_lvls = prod(size(ncid('lev')));
+ dnum_lons = dim_length(fname,'lon');
+ dnum_lats = dim_length(fname,'lat');
+ dnum_lvls = dim_length(fname,'lev');
x = 3;
y = [dnum_lons dnum_lats dnum_lvls];
- end
case 'pe2lyr'
- lonexist = VarExist(ncid,'lon');
- latexist = VarExist(ncid,'lat');
- lvlexist = VarExist(ncid,'level');
- if ( latexist && lonexist && lvlexist )
- dnum_lons = prod(size(ncid('lon')));
- dnum_lats = prod(size(ncid('lat')));
- dnum_lvls = prod(size(ncid('lev')));
+ dnum_lons = dim_length(fname,'lon');
+ dnum_lats = dim_length(fname,'lat');
+ dnum_lvls = dim_length(fname,'lev');
x = 3;
y = [dnum_lons dnum_lats dnum_lvls];
- end
case 'fms_bgrid'
- lonexist = VarExist(ncid,'TmpI' );
- latexist = VarExist(ncid,'TmpJ' );
- lvlexist = VarExist(ncid,'level');
- if ( latexist && lonexist && lvlexist )
- dnum_lons = prod(size(ncid('TmpI')));
- dnum_lats = prod(size(ncid('TmpJ')));
- dnum_lvls = prod(size(ncid('level')));
+ dnum_lons = dim_length(fname,'TmpI');
+ dnum_lats = dim_length(fname,'TmpJ');
+ dnum_lvls = dim_length(fname,'lev');
x = 3;
y = [dnum_lons dnum_lats dnum_lvls];
- end
case 'mitgcm_ocean'
- lonexist = VarExist(ncid,'XG' );
- latexist = VarExist(ncid,'YG' );
- lvlexist = VarExist(ncid,'ZG');
- if ( latexist && lonexist && lvlexist )
- dnum_lons = prod(size(ncid('XG')));
- dnum_lats = prod(size(ncid('YG')));
- dnum_lvls = prod(size(ncid('ZG')));
+ dnum_lons = dim_length(fname,'XG');
+ dnum_lats = dim_length(fname,'YG');
+ dnum_lvls = dim_length(fname,'ZG');
x = 3;
y = [dnum_lons dnum_lats dnum_lvls];
- end
case 'lorenz_96_2scale'
- Xexist = VarExist(ncid,'Xdim');
- Yexist = VarExist(ncid,'Ydim');
- if ( Xexist && Yexist )
- dnum_X = prod(size(ncid('Xdim')));
- dnum_Y = prod(size(ncid('Ydim')));
+ dnum_X = dim_length(fname,'Xdim');
+ dnum_Y = dim_length(fname,'Ydim');
x = 2;
y = [dnum_X dnum_Y];
- end
case 'simple_advection'
- if ( VarExist(ncid,'loc1d'))
- y = prod(size(ncid('loc1d')));
+ y = dim_length(fname,'loc1d');
x = 1;
- end
otherwise
- if ( VarExist(ncid,'StateVariable'))
- y = prod(size(ncid('StateVariable')));
+ y = dim_length(fname,'StateVariable');
x = 1;
- end
end
Modified: DART/trunk/matlab/ChecknetCDFuse.m
===================================================================
--- DART/trunk/matlab/ChecknetCDFuse.m 2009-07-10 21:24:09 UTC (rev 3962)
+++ DART/trunk/matlab/ChecknetCDFuse.m 2009-07-15 22:24:43 UTC (rev 3963)
@@ -12,7 +12,7 @@
%
% Data Assimilation Research Testbed -- DART
-% Copyright 2004-2007, Data Assimilation Research Section
+% Copyright 2004-2009, Data Assimilation Research Section
% University Corporation for Atmospheric Research
% Licensed under the GPL -- www.gpl.org/licenses/gpl.html
%
@@ -22,15 +22,17 @@
% $Revision$
% $Date$
+if ( exist(fname,'file') ~= 2 ), error('%s does not exist.',fname); end
+
fid = fopen(fname,'wt');
-fprintf(fid,'%s \n',datestr(now));
-fprintf(fid,'%s \n',pwd);
-if (exist('getnc') ~= 2)
- fprintf(fid,'%s \n','No matlab netcdf operators, not using matlab.');
+fprintf(fid,'%s\n',datestr(now));
+fprintf(fid,'%s\n',pwd);
+if (exist('nc_varget') ~= 2)
+ fprintf(fid,'%s \n','No matlab snctools, not using matlab.');
fprintf(fid,'%d \n',-1);
x = -1;
else
- fprintf(fid,'%s \n','Found matlab netcdf operators, using matlab.');
+ fprintf(fid,'%s \n','Found matlab snctools, using matlab.');
fprintf(fid,'%d \n',0);
x = 0;
end
Modified: DART/trunk/matlab/CombineStructs.m
===================================================================
--- DART/trunk/matlab/CombineStructs.m 2009-07-10 21:24:09 UTC (rev 3962)
+++ DART/trunk/matlab/CombineStructs.m 2009-07-15 22:24:43 UTC (rev 3963)
@@ -6,7 +6,7 @@
%
% Data Assimilation Research Testbed -- DART
-% Copyright 2004-2007, Data Assimilation Research Section
+% Copyright 2004-2009, Data Assimilation Research Section
% University Corporation for Atmospheric Research
% Licensed under the GPL -- www.gpl.org/licenses/gpl.html
%
Modified: DART/trunk/matlab/DiffnetCDFstate.m
===================================================================
--- DART/trunk/matlab/DiffnetCDFstate.m 2009-07-10 21:24:09 UTC (rev 3962)
+++ DART/trunk/matlab/DiffnetCDFstate.m 2009-07-15 22:24:43 UTC (rev 3963)
@@ -17,7 +17,7 @@
% DiffnetCDFstate(file1,file2,outfile)
% Data Assimilation Research Testbed -- DART
-% Copyright 2004-2007, Data Assimilation Research Section
+% Copyright 2004-2009, Data Assimilation Research Section
% University Corporation for Atmospheric Research
% Licensed under the GPL -- www.gpl.org/licenses/gpl.html
%
@@ -31,30 +31,23 @@
% Check for the existence of the two input files.
%----------------------------------------------------------------------
-if (exist(file1,'file') ~= 2)
- error(sprintf('%s does not exist.',file1))
-end
-if (exist(file2,'file') ~= 2)
- error(sprintf('%s does not exist.',file2))
-end
+if (exist(file1,'file') ~= 2), error('%s does not exist.',file1); end
+if (exist(file2,'file') ~= 2), error('%s does not exist.',file2); end
%----------------------------------------------------------------------
% if the 'state' variable exists, read in the WHOLE THING.
%----------------------------------------------------------------------
-f1 = netcdf(file1);
-f2 = netcdf(file2);
-
-if ( VarExist(f1,'state'))
- a = getnc(file1,'state');
+if ( nc_isvar(file1,'state'))
+ a = nc_varget(file1,'state');
else
- error(sprintf('%s has no ''state'' variable.',file1))
+ error('%s has no ''state'' variable.',file1)
end
-if ( VarExist(f2,'state'))
- b = getnc(file2,'state');
+if ( nc_isvar(file2,'state'))
+ b = nc_varget(file2,'state');
else
- error(sprintf('%s has no ''state'' variable.',file2))
+ error('%s has no ''state'' variable.',file2)
end
% string them into 1D arrays and take the difference
@@ -65,15 +58,15 @@
a = min(c);
b = max(c);
-clear c f1 f2
+clear c
% Write the min and max of the differences to a file
fid = fopen(outfile,'wt');
-fprintf(fid,'%s \n',datestr(now));
-fprintf(fid,'%s \n',pwd);
-fprintf(fid,'%s \n',file1);
-fprintf(fid,'%s \n',file2);
+fprintf(fid,'%s\n',datestr(now));
+fprintf(fid,'%s\n',pwd);
+fprintf(fid,'%s\n',file1);
+fprintf(fid,'%s\n',file2);
if ( (a == 0.0) && (b == 0.0) )
fprintf(fid,'0 %e %e\n',a,b);
@@ -83,14 +76,3 @@
fclose(fid);
-
-function x = VarExist(ncid,varname)
-
-x = 0; % false ... assumed not to exist.
-variables = var(ncid);
-for i=1:length(variables)
- if ( strmatch(name(variables{i}), varname) == 1 )
- x = 1; % true ... variables exists in the netcdf file.
- end
-end
-
Modified: DART/trunk/matlab/GetBgridInfo.m
===================================================================
--- DART/trunk/matlab/GetBgridInfo.m 2009-07-10 21:24:09 UTC (rev 3962)
+++ DART/trunk/matlab/GetBgridInfo.m 2009-07-15 22:24:43 UTC (rev 3963)
@@ -9,7 +9,7 @@
% routine name of subsequent plot routine.
% Data Assimilation Research Testbed -- DART
-% Copyright 2004-2007, Data Assimilation Research Section
+% Copyright 2004-2009, Data Assimilation Research Section
% University Corporation for Atmospheric Research
% Licensed under the GPL -- www.gpl.org/licenses/gpl.html
%
@@ -19,22 +19,22 @@
% $Revision$
% $Date$
-pinfo = pinfo_in;
-ft = netcdf(fname);
-model = ft.model(:);
-close(ft)
+if ( exist(fname,'file') ~= 2 ), error('%s does not exist.',fname); end
+pinfo = pinfo_in;
+model = nc_attget(fname, nc_global, 'model');
+
if strcmp(lower(model),'fms_bgrid') ~= 1
error('Not so fast, this is not a bgrid model.')
end
-copy = getnc(fname,'copy');
-times = getnc(fname,'time');
-levels = getnc(fname,'lev');
-TmpI = getnc(fname,'TmpI'); % temperature/pressure grid longitude
-TmpJ = getnc(fname,'TmpJ'); % temperature/pressure grid latitude
-VelI = getnc(fname,'VelI'); % velocity grid longitude
-VelJ = getnc(fname,'VelJ'); % velocity grid latitude
+copy = nc_varget(fname,'copy');
+times = nc_varget(fname,'time');
+levels = nc_varget(fname,'lev');
+TmpI = nc_varget(fname,'TmpI'); % temperature/pressure grid longitude
+TmpJ = nc_varget(fname,'TmpJ'); % temperature/pressure grid latitude
+VelI = nc_varget(fname,'VelI'); % velocity grid longitude
+VelJ = nc_varget(fname,'VelJ'); % velocity grid latitude
prognostic_vars = {'ps','t','u','v'};
@@ -131,7 +131,7 @@
pinfo = setfield(pinfo, 'model' , model);
pinfo = setfield(pinfo, 'var_names' , pgvar);
- pinfo = setfield(pinfo, 'truth_file' , []);
+ %pinfo = setfield(pinfo, 'truth_file' , []);
%pinfo = setfield(pinfo, 'prior_file' , pinfo.prior_file);
%pinfo = setfield(pinfo, 'posterior_file', pinfo.posterior_file);
pinfo = setfield(pinfo, 'level' , level);
@@ -209,8 +209,8 @@
for i = 2:length(prognostic_vars),
str = sprintf(' %s %s ',str,prognostic_vars{i});
end
-disp(sprintf('Default variable is ''%s'', if this is OK, <cr>;',pgvar))
-disp(sprintf('If not, please enter one of: %s',str))
+fprintf('Default variable is ''%s'', if this is OK, <cr>;\n',pgvar)
+fprintf('If not, please enter one of: %s\n',str)
varstring = input('(no syntax required)\n','s');
if ~isempty(varstring), pgvar = deblank(varstring); end
@@ -221,9 +221,9 @@
%----------------------------------------------------------------------
if (nargin == 3), time = deftime; else time = mean(times); end
-disp(sprintf('Default time is %f, if this is OK, <cr>;',time))
-disp(sprintf('If not, enter a time between %.4f and %.4f, we use the closest.', ...
- min(times),max(times)))
+fprintf('Default time is %f, if this is OK, <cr>;\n',time)
+fprintf('If not, enter a time between %.4f and %.4f, we use the closest.\n', ...
+ min(times),max(times))
varstring = input('(no syntax required)\n','s');
if ~isempty(varstring), time = str2num(varstring); end
@@ -245,9 +245,9 @@
level = 1;
lvlind = 1;
else
- disp(sprintf('Default level is %d, if this is OK, <cr>;',level))
- disp(sprintf('If not, enter a level between %d and %d, inclusive ...', ...
- min(levels),max(levels)))
+ fprintf('Default level is %d, if this is OK, <cr>;\n',level)
+ fprintf('If not, enter a level between %d and %d, inclusive ...\n', ...
+ min(levels),max(levels))
varstring = input('we''ll use the closest (no syntax required)\n','s');
if ~isempty(varstring), level = str2num(varstring); end
@@ -271,9 +271,9 @@
lons = VelI;
end
-disp(sprintf('Default longitude is %f, if this is OK, <cr>;',lon))
-disp(sprintf('If not, enter a longitude between %.2f and %.2f, we use the closest.', ...
- min(lons),max(lons)))
+fprintf('Default longitude is %f, if this is OK, <cr>;\n',lon)
+fprintf('If not, enter a longitude between %.2f and %.2f, we use the closest.\n', ...
+ min(lons),max(lons))
varstring = input('(no syntax required)\n','s');
if ~isempty(varstring), lon = str2num(varstring); end
@@ -296,9 +296,9 @@
lats = VelJ;
end
-disp(sprintf('Default latitude is %f, if this is OK, <cr>;',lat))
-disp(sprintf('If not, enter a latitude between %.2f and %.2f, we use the closest.', ...
- min(lats),max(lats)))
+fprintf('Default latitude is %f, if this is OK, <cr>;\n',lat)
+fprintf('If not, enter a latitude between %.2f and %.2f, we use the closest.\n', ...
+ min(lats),max(lats))
varstring = input('(no syntax required)\n','s');
if ~isempty(varstring), lat = str2num(varstring); end
Modified: DART/trunk/matlab/GetCalendarDate.m
===================================================================
--- DART/trunk/matlab/GetCalendarDate.m 2009-07-10 21:24:09 UTC (rev 3962)
+++ DART/trunk/matlab/GetCalendarDate.m 2009-07-15 22:24:43 UTC (rev 3963)
@@ -9,7 +9,7 @@
% mydate = GetCalendarDate(82761,148520);
% Data Assimilation Research Testbed -- DART
-% Copyright 2004-2007, Data Assimilation Research Section
+% Copyright 2004-2009, Data Assimilation Research Section
% University Corporation for Atmospheric Research
% Licensed under the GPL -- www.gpl.org/licenses/gpl.html
%
@@ -31,8 +31,8 @@
case 'gregorian'
mytime = datenum(1601,1,1) + days + seconds/86400;
h = datestr(mytime);
- disp(sprintf('DART time (%d s, %d d) is %s %s', ...
- seconds,days,h,calendartype))
+ fprintf('DART time (%d s, %d d) is %s %s\n', ...
+ seconds,days,h,calendartype)
case 'noleap'
error('noleap not supported yet')
case 'thirty_day_months'
Modified: DART/trunk/matlab/GetCamInfo.m
===================================================================
--- DART/trunk/matlab/GetCamInfo.m 2009-07-10 21:24:09 UTC (rev 3962)
+++ DART/trunk/matlab/GetCamInfo.m 2009-07-15 22:24:43 UTC (rev 3963)
@@ -8,7 +8,7 @@
% routine name of subsequent plot routine.
% Data Assimilation Research Testbed -- DART
-% Copyright 2004-2007, Data Assimilation Research Section
+% Copyright 2004-2009, Data Assimilation Research Section
% University Corporation for Atmospheric Research
% Licensed under the GPL -- www.gpl.org/licenses/gpl.html
%
@@ -18,30 +18,28 @@
% $Revision$
% $Date$
-if ( exist(pstruct.truth_file) )
+if ( exist(pstruct.truth_file,'file') )
fname = pstruct.truth_file;
-elseif ( exist(pstruct.diagn_file) )
+elseif ( exist(pstruct.diagn_file,'file') )
fname = pstruct.diagn_file;
-elseif ( exist(pstruct.prior_file) )
+elseif ( exist(pstruct.prior_file,'file') )
fname = pstruct.prior_file;
-elseif ( exist(pstruct.posterior_file) )
+elseif ( exist(pstruct.posterior_file,'file') )
fname = pstruct.posterior_file;
end
-ft = netcdf(fname);
-model = ft.model(:);
-close(ft)
+model = nc_attget(fname,nc_global,'model');
if strcmp(lower(model),'cam') ~= 1
error('Not so fast, this is not a cam model.')
end
-copy = getnc(fname,'copy');
-times = getnc(fname,'time');
-ilevel = getnc(fname,'ilev'); % interfaces
-levels = getnc(fname, 'lev'); % midpoints
-lon = getnc(fname, 'lon');
-lat = getnc(fname, 'lat');
+copy = nc_varget(fname,'copy');
+times = nc_varget(fname,'time');
+ilevel = nc_varget(fname,'ilev'); % interfaces
+levels = nc_varget(fname, 'lev'); % midpoints
+lon = nc_varget(fname, 'lon');
+lat = nc_varget(fname, 'lat');
% A more robust way would be to use the netcdf low-level ops:
% bob = var(f); % bob is a cell array of ncvars
@@ -192,8 +190,8 @@
for i = 2:length(prognostic_vars),
str = sprintf(' %s %s ',str,prognostic_vars{i});
end
-disp(sprintf('Default variable is ''%s'', if this is OK, <cr>;',pgvar))
-disp(sprintf('If not, please enter one of: %s',str))
+fprintf('Default variable is ''%s'', if this is OK, <cr>;\n',pgvar)
+fprintf('If not, please enter one of: %s\n',str)
varstring = input('(no syntax required)\n','s');
if ~isempty(varstring), pgvar = strtrim(varstring); end
@@ -206,9 +204,9 @@
%----------------------------------------------------------------------
if (nargin == 3), time = deftime; else time = mean(times); end
-disp(sprintf('Default time is %f, if this is OK, <cr>;',time))
-disp(sprintf('If not, enter a time between %.4f and %.4f, we use the closest.', ...
- min(times),max(times)))
+fprintf('Default time is %f, if this is OK, <cr>;\n',time)
+fprintf('If not, enter a time between %.4f and %.4f, we use the closest.\n', ...
+ min(times),max(times))
varstring = input('(no syntax required)\n','s');
if ~isempty(varstring), time = str2num(varstring); end
@@ -230,9 +228,9 @@
level = 1;
lvlind = 1;
else
- disp(sprintf('Default level (index) is %d, if this is OK, <cr>;',lvlind))
- disp(sprintf('If not, enter a level between %d and %d, inclusive ...', ...
- 1,length(levels)))
+ fprintf('Default level (index) is %d, if this is OK, <cr>;\n',lvlind)
+ fprintf('If not, enter a level between %d and %d, inclusive ...\n', ...
+ 1,length(levels))
varstring = input('we''ll use the closest (no syntax required)\n','s');
if ~isempty(varstring), lvlind = str2num(varstring); end
@@ -249,9 +247,9 @@
%----------------------------------------------------------------------
if (nargin == 3), lon = deflon; else lon = 255.0; end
-disp(sprintf('Default longitude is %f, if this is OK, <cr>;',lon))
-disp(sprintf('If not, enter a longitude between %.2f and %.2f, we use the closest.', ...
- min(lons),max(lons)))
+fprintf('Default longitude is %f, if this is OK, <cr>;\n',lon)
+fprintf('If not, enter a longitude between %.2f and %.2f, we use the closest.\n', ...
+ min(lons),max(lons))
varstring = input('(no syntax required)\n','s');
if ~isempty(varstring), lon = str2num(varstring); end
@@ -267,9 +265,9 @@
%----------------------------------------------------------------------
if (nargin == 3), lat = deflat; else lat = 40.0; end
-disp(sprintf('Default latitude is %f, if this is OK, <cr>;',lat))
-disp(sprintf('If not, enter a latitude between %.2f and %.2f, we use the closest.', ...
- min(lats),max(lats)))
+fprintf('Default latitude is %f, if this is OK, <cr>;\n',lat)
+fprintf('If not, enter a latitude between %.2f and %.2f, we use the closest.\n', ...
+ min(lats),max(lats))
varstring = input('(no syntax required)\n','s');
if ~isempty(varstring), lat = str2num(varstring); end
Modified: DART/trunk/matlab/GetNCindices.m
===================================================================
--- DART/trunk/matlab/GetNCindices.m 2009-07-10 21:24:09 UTC (rev 3962)
+++ DART/trunk/matlab/GetNCindices.m 2009-07-15 22:24:43 UTC (rev 3963)
@@ -1,23 +1,29 @@
-function [corner, endpnt] = GetNCindices(pinfo, whichfile, varname);
-% GETNCindices returns a corner,endpnt array for use with getnc.
+function [start, count] = GetNCindices(pinfo, whichfile, varname);
+% GETNCindices returns a start,count array for use with nc_getvar.
% At present, all times, all copies for a specific level,lat,lon.
% Does not assume anything about the dimension of the variable.
%
% USAGE:
-% [corner, endpnt] = GetNCindices(pinfo, whichfile, varname);
+% [start, count] = GetNCindices(pinfo, whichfile, varname);
+%
+% The structure 'pinfo' must have ONE of the following components:
+% pinfo.[prior,posterior,truth,diagn,fname]
+% and may have
+% pinfo.timeindex
+% pinfo.copyindex
+% pinfo.levelindex
+% pinfo.latindex
+% pinfo.lonindex
+% pinfo.stateindex
+%
% whichfile is a character string specifying which
-% component of 'pinfo' will be used.
-% ['prior','posterior','truth']
+% filename component of 'pinfo' will be used.
+% ['prior','posterior','truth','diagn','fname']
% varname is the netcdf variable being extracted.
+%
-% The structure 'pinfo' must have the following components:
-%pinfo.latindex
-%pinfo.lonindex
-%pinfo.levelindex
-% and the value
-
% Data Assimilation Research Testbed -- DART
-% Copyright 2004-2007, Data Assimilation Research Section
+% Copyright 2004-2009, Data Assimilation Research Section
% University Corporation for Atmospheric Research
% Licensed under the GPL -- www.gpl.org/licenses/gpl.html
%
@@ -30,11 +36,11 @@
% GetNCindices replaces the following hardwired piece of code.
%
% if ( strcmp(lower(vname),'ps') ==1 ) % PS(time, copy, lat, lon)
-% corner = [ 1 NaN pinfo.latindex pinfo.lonindex];
-% endpnt = [ -1 NaN pinfo.latindex pinfo.lonindex];
+% start = [ 1 NaN pinfo.latindex pinfo.lonindex];
+% count = [ -1 NaN pinfo.latindex pinfo.lonindex];
% else % U(time, copy, lev, lat, lon)
-% corner = [ 1 NaN pinfo.levelindex pinfo.latindex pinfo.lonindex ];
-% endpnt = [ -1 NaN pinfo.levelindex pinfo.latindex pinfo.lonindex ];
+% start = [ 1 NaN pinfo.levelindex pinfo.latindex pinfo.lonindex ];
+% count = [ -1 NaN pinfo.levelindex pinfo.latindex pinfo.lonindex ];
% end
switch lower(whichfile)
@@ -50,72 +56,122 @@
fname = pinfo.fname;
end
-corner = [];
-endpnt = [];
+if ( exist(fname,'file') ~= 2 ), error('%s does not exist.',fname); end
-% This block uses the lowest level netcdf operator I know of.
+% If the structure has subsetting information, use it.
+% Otherwise, use the whole extent.
-[cdfid, rcode] = ncmex('OPEN',fname,'NC_NOWRITE');
-[varid, rcode] = ncmex('VARID', cdfid, varname);
-[var, datatype, ndims, dim, natts, status] = ncmex('VARINQ', cdfid, varid);
+lat1 = 0; latN = -1;
+lon1 = 0; lonN = -1;
+time1 = 0; timeN = -1;
+copy1 = 0; copyN = -1;
+level1 = 0; levelN = -1;
+state1 = 0; stateN = -1;
-% var is the name of the variable as it is in the netCDF file
-% datatype is
-% ndims is the 'rank' of the variable 2D, 3D, 4D, etc.
-% dims is an array of length 'ndims' specifying the
-% dimension ID's used - IN ORDER.
-% natts is the number of attributes for the variable
-% status is obvious
+if (isfield(pinfo,'timeindex'))
+ time1 = pinfo.timeindex - 1;
+ timeN = 1;
+end
+if (isfield(pinfo,'levelindex'))
+ level1 = pinfo.levelindex - 1;
+ levelN = 1;
+end
+if (isfield(pinfo,'latindex'))
+ lat1 = pinfo.latindex - 1;
+ latN = 1;
+end
+if (isfield(pinfo,'lonindex'))
+ lon1 = pinfo.lonindex - 1;
+ lonN = 1;
+end
+if (isfield(pinfo,'stateindex'))
+ state1 = pinfo.stateindex - 1;
+ stateN = 1;
+end
+if (isfield(pinfo,'copyindex'))
+ copy1 = pinfo.copyindex - 1;
+ copyN = 1;
+end
-% make sure it has a time and copy index
+% Determine shape of variable in question.
+varinfo = nc_getvarinfo(fname,varname);
+ndims = length(varinfo.Dimension);
+start = zeros(1,ndims);
+count = zeros(1,ndims);
+
+% varinfo.Dimension is a cell array of the Dimension strings
+% varinfo.Size is an N-D array describing the variable shape
@@ Diff output truncated at 40000 characters. @@
More information about the Dart-dev
mailing list