read_netCDF_FVCOM.m 29.8 KB
Newer Older
1
function [data,selection] =read_netCDF_FVCOM(varargin)
2
% Function to extract data from a Netcdf file output from FVCOM.
3
%
4
% data = read_netCDF_FVCOM(varargin)
5 6
%
% DESCRIPTION:
7
%    Function to extract data from a netCDF file output from FVCOM. Outputs
8
%    data in cell array.
9
%
10
% INPUT [keyword pairs]:
11 12
%   Options are passed in pairs.
%
13 14 15 16 17 18 19 20 21 22 23 24
%   The list of keywords is:
%       - 'time'
%       - 'data_dir'
%       - 'file_netcdf'
%       - 'varnames'
%       - 'nele_idx'
%       - 'node_idx'
%       - 'siglay_idx'
%       - 'siglev_idx'
%
%   'time' - {'30/01/06 00:00:00', '01/02/06 23:00:00'} or -1 to extract
%   all times.
25
%
26 27
%   'data_dir' - '/home/fvcom/results/...' directory where netCDF file is.
%   Default value is ../fvcom_postproc/netcdf
28
%
29 30
%   'file_netcdf' - 'filename.nc'. Default value is '*.nc', but it only
%   access the first file in alphabetical order in the directory.
31
%
32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49
%   'varnames' - Cell array of variable names to read from the netCDF file.
%   The variables need to exist in the file but they are case insensitive.
%   Choose FVCOM output variable names. For example:
%       - 'Itime'
%       - 'Itime2'
%       - 'xc'
%       - 'yc'
%       - 'art1'
%       - 'art2'
%       - 'h'
%       - 'siglay'
%       - 'siglev'
%       - 'nv'
%       - 'zeta'
%       - 'ua'
%       - 'va'
%   The complete list for a given file is given by running this script with
%   varnames set to [].
50
%
51 52 53 54 55 56 57
%   The variables can be restricted in five possible dimensions:
%       - 'node_idx'
%       - 'nele_idx'
%       - 'siglev_idx'
%       - 'siglay_idx'
%       - 'time_idx'
%   Default values cause the script to extract all available data for all
58 59
%   possible dimensions. No checks are done on the bounds of each dimension
%   so make sure you choose them right!
60 61
%
% OUTPUT:
62 63
%    data = struct with fields whose names match those from the list of
%    input variables extracted ('varnames').
64 65
%
% EXAMPLE USAGE
66 67 68 69 70 71 72 73 74 75
%   vars = {'Times', 'xc', 'yc', 'h', 'siglay', 'nv', 'zeta', 'ua', 'va'};
%   date_range = {'30/01/06 00:00:00', '15/02/06 23:00:00'};
%   node_idx = [10:30, 40:50]; % zero referenced!
%   data_dir = '/home/fvcom/results/output/';
%   FVCOM = read_netCDF_FVCOM('data_dir', data_dir, ...
%       'file_netcdf', 'casename_0001.nc', ...
%       'time', date_range, ...
%       'siglev_idx', 1, ...
%       'node_idx', node_idx, ...
%       'varnames', vars);
76
%
77
% Author(s):
78
%    Ricardo Torres - Plymouth Marine Laboratory 2012
79
%    Hakeem Johnson - CH2M
80
%    Pierre Cazenave - Plymouth Marine Laboratory
81
%
82
% Revision history:
83
%   v0 March 2012
84 85
%   2014-03-06 - Add the global verbose flag. Also tidy up the help a bit.
%   Also change some verbose statements to use fprintf instead of disp for
86 87 88
%   better control over formatting. Also fixed a bug where if a 2D array
%   was requested after a 3D array, the 2D array would cause the function
%   to crash (because it was using a 3D index for getVar).
89 90 91 92
%   2014-08-20 - Complete the functionality to be able to slice the data
%   along any dimension (siglay, time, node etc.).
%   2014-10-17 - Fix ability to slice with any combination of space
%   (horizontal and vertical) and time.
93 94
%
%==========================================================================
95 96 97 98 99 100 101 102 103

global ftbverbose
subname = 'read_netCDF_FVCOM';

if ftbverbose
    fprintf('\nbegin : %s \n', subname)
end

%--------------------------------------------------------------------------
104
%  Parse input arguments
105
%--------------------------------------------------------------------------
106

107
params_opts = {'time', 'data_dir', 'file_netcdf', 'varnames', 'nele_idx', ...
108
    'node_idx', 'siglay_idx', 'siglev_idx', 'timestride'};
109 110

if ftbverbose
111
    fprintf('Input parameters being used are:\n')
112 113
end
var_in_list = {'all_data', 'netfile_dir', 'file_netcdf', 'varnames', ...
114
    'nele_idx', 'node_idx', 'siglay_idx', 'siglev_idx', 'timestrd'};
115 116 117 118 119
all_data = 1;
netfile_dir = '../fvcom_postproc/netcdf';
file_netcdf='*.nc';
siglay_idx=-1;
siglev_idx=-1;
120 121
nele_idx=-1;
node_idx=-1;
122 123
time_idx=-1;
varnames={};
124
timestrd=1;
125 126
for aa=1:2:nargin
    res=strcmp(varargin(aa),params_opts);
127
    if sum(res)
128
        eval([var_in_list{res},' = varargin{aa+1};'])
129 130 131
        if ftbverbose
            fprintf(' %s\n', params_opts{res})
        end
132 133
    end
end
134 135 136 137

%--------------------------------------------------------------------------
% Sort (and remove repeats) for all indices elements, nodes or layers
%--------------------------------------------------------------------------
138 139 140 141
nele_idx=unique(nele_idx);
node_idx=unique(node_idx);
siglay_idx=unique(siglay_idx);
siglev_idx=unique(siglev_idx);
142

143 144
RestrictDims.Name={'node' 'nele' 'siglay' 'siglev' 'time'};
RestrictDims.idx={node_idx, nele_idx, siglay_idx, siglev_idx, time_idx};
145

146 147 148
if ~isempty(varnames)
    nvarnames = length(varnames);
    for nn=1:nvarnames
149
        data.(varnames{nn}) = [];
150 151
    end
end
152 153

%--------------------------------------------------------------------------
154
% Open netcdf file
155 156
%--------------------------------------------------------------------------
file_netcdf=fullfile(netfile_dir, file_netcdf);
157 158 159
filesINdir=dir(file_netcdf);
file_netcdf= fullfile(netfile_dir,filesINdir(1).name);
nc = netcdf.open(file_netcdf, 'NC_NOWRITE');
160 161 162 163 164 165 166 167
if ftbverbose
    if length(file_netcdf) > 50
        % Truncate output file name to display.
        fprintf('NetCDF file ...%s opened successfully.\n', file_netcdf(end-70:end))
    else
        fprintf('NetCDF file %s opened successfully.\n', file_netcdf)
    end
end
168 169 170 171 172 173 174
% Get information from netcdf file
info=ncinfo(file_netcdf);
% Extract all possible dimensions in file
DimsAll=info.Dimensions;
% Extract variable names in  nc file
Vars=struct2cell(info.Variables);
vars = squeeze(Vars(1,:,:));
175 176 177 178 179 180 181 182

%--------------------------------------------------------------------------
% Find variable Itime
%--------------------------------------------------------------------------
if ftbverbose
    fprintf('Using date conversion of +678942 days to go from FVCOM time (Modified Julian Day) to MATLAB time.\n')
end
time_offset = 678942;
183 184 185
idx=find(strcmpi(cat(1,{DimsAll.Name}),'time'));
last_entry=DimsAll(idx).Length;
Itime=[];Itime2=[];
186
% tic
187 188 189
try
    Itime.idx=find(strcmpi(vars,'Itime'));
    Itime.ID=netcdf.inqVarID(nc,'Itime');
190 191 192 193
    Itime.Data(1)  = netcdf.getVar(nc,Itime.ID,0,1,'int32');
    Itime.Data(2)  = netcdf.getVar(nc,Itime.ID,last_entry-1,1,'int32');
    Itime2.Data(1)  = netcdf.getVar(nc,Itime.ID+1,0,1,'int32');
    Itime2.Data(2)  = netcdf.getVar(nc,Itime.ID+1,last_entry-1,1,'int32');
Pierre Cazenave's avatar
Pierre Cazenave committed
194

195 196
    [start_d(1),end_d(1)] = deal(double(Itime.Data(1))+time_offset,double(Itime.Data(end))+time_offset);
    [start_d(2),end_d(2)] = deal(double(Itime2.Data(1)),double(Itime2.Data(end)));
Pierre Cazenave's avatar
Pierre Cazenave committed
197

198 199
    start_date=sum(start_d.*[1 1/(24*60*60*1000)]);     %hkj missing 1000 inserted
    end_date = sum(end_d.*[1 1/(24*60*60*1000)]);       %hkj missing 1000 inserted
200 201
    var_time =  netcdf.getVar(nc,Itime.ID,[0],[10],'double')+time_offset+...
        netcdf.getVar(nc,Itime.ID+1,0,10,'double')./(24*600*6000) ;
Pierre Cazenave's avatar
Pierre Cazenave committed
202

203 204
    DeltaT=median(diff(var_time));
    var_time = start_date:DeltaT:(end_date-DeltaT);
Pierre Cazenave's avatar
Pierre Cazenave committed
205

206
catch me
207
    if ftbverbose
208
        warning('No ''Itime'' and/or ''Itime2'' variables, using less precise ''time'' instead.\n(%s)\n', me.message)
209
    end
210 211
    Itime.idx=find(strcmpi(vars,'time'));
    Itime.ID=netcdf.inqVarID(nc,'time');
212 213 214
    Itime.Data(1)  = netcdf.getVar(nc,Itime.ID,0,1,'double');
    Itime.Data(2)  = netcdf.getVar(nc,Itime.ID,last_entry-1,1,'double');
    [start_date,end_date] = deal(Itime.Data(1)+time_offset,Itime.Data(end)+time_offset);
215
    DeltaT=(end_date-start_date)./last_entry;
216
    var_time = start_date:DeltaT:(end_date-DeltaT);
217
end
218
% toc
219
if length(all_data) == 2
220 221 222
    req_st = datenum(all_data{1},'dd/mm/yy HH:MM:SS');
    req_end = datenum(all_data{2},'dd/mm/yy HH:MM:SS');
else
223
    req_st = start_date;
224 225 226
    req_end =end_date;
end
time_idx = find(req_st <= var_time &   var_time <= req_end );
227
time_idx = time_idx(1:timestrd:end);
228
% Add correct time_idx to RestrictDims
229
RestrictDims.idx{end}=time_idx;
230 231 232 233 234
if ftbverbose
    fprintf('Start and end of file: %s - %s\n', datestr(start_date), datestr(end_date))
end

%--------------------------------------------------------------------------
235
% Return information about file to the screen
236
%--------------------------------------------------------------------------
237

238 239 240 241
if ftbverbose
    fprintf('Possible variables to extract are:\n')
end
for ii = 1:length(vars)
242 243 244
    if ftbverbose
        fprintf(' %s\n', vars{ii})
    end
245 246
end
if isempty(varnames)
247
    data = 0;
248
    netcdf.close(nc)
249
    error('Stopping. Choose a variable from the list above.')
250
end
251 252

%--------------------------------------------------------------------------
253
% Re-organise RestrictDims to follow order of dimensions in nc file from
254 255
% FVCOM
%--------------------------------------------------------------------------
256 257 258 259 260 261 262 263 264
cc=1;
for dd=1:length(DimsAll)
    idx=find(strcmpi(RestrictDims.Name,DimsAll(dd).Name));
    if ~isempty(idx)
        TEMP{cc}=RestrictDims.Name{idx};
        TEMPidx{cc}=RestrictDims.idx{idx};
        cc=cc+1;
    end
end
265 266 267
RestrictDims.Name = TEMP;
RestrictDims.idx = TEMPidx;
clear TEMP TEMPidx
268 269

%--------------------------------------------------------------------------
270
% Start Processing extraction of data from NC file
271
%--------------------------------------------------------------------------
272
selection=[];
273
for aa=1:length(varnames)
274 275
    selection.(varnames{aa}).start=-1;
    selection.(varnames{aa}).count=-1;
276 277 278
    %----------------------------------------------------------------------
    % Extract number of dimensions, lengths and names of all variables
    %----------------------------------------------------------------------
279
% tic
280 281 282
    if ftbverbose
        fprintf('Processing variable %s: ', varnames{aa})
    end
283 284
    % Tidy up the previous iteration's variables so we don't get confused.
    clear dimName dimLength
285

286 287
    TF = strcmpi(varnames{aa},vars);
    if ~isempty(find(TF));
288 289 290
        varidx(aa) = find(TF);
        TF = sum(TF);
        dimens=ndims(aa);
291 292 293
    else
        netcdf.close(nc)
        varargout{1} = 0;
294
        if ftbverbose; fprintf('\n'); end
295
        error('Variable %s NOT found in file. Stopping. Check input variable names.', varnames{aa})
296 297
    end
    varID=netcdf.inqVarID(nc,vars{varidx(aa)});
298

299 300
    [name,xtype,dimids,natts] = netcdf.inqVar(nc,varID);
    dimens=length(dimids);
301

302
    for dd=1:length(dimids)
303 304 305 306
        [dimName{dd}, dimLength(dd)] = netcdf.inqDim(nc,dimids(dd));
        if ftbverbose
            if dd == 1
                if length(dimids) == 1
307 308 309
                    if ftbverbose
                        fprintf('%i dimension: %s ', dimens, dimName{dd})
                    end
310
                else
311 312 313
                    if ftbverbose
                        fprintf('%i dimensions: %s ', dimens, dimName{dd})
                    end
314 315
                end
            else
316 317 318
                if ftbverbose
                    fprintf('%s ', dimName{dd})
                end
319 320
            end
        end
321
    end
322
    if ftbverbose; fprintf('\n'); end
323

324 325 326
    %----------------------------------------------------------------------
    % Get the data!
    %----------------------------------------------------------------------
327

328 329 330 331 332 333
    switch dimens
        case 1
            % only one dimension present in variable
            switch dimName{1}
                case 'time'
                    if time_idx>=0
334
                        % Only restrict data on access if dimension is TIME
335

336 337 338
                        % hkj it appears the first value in matlab netcdf
                        % interface is 0.
                        % hkj time_idx(1) CORRECTED TO time_idx(1)-1.
339
                        eval([varnames{aa},'=netcdf.getVar(nc,varID,time_idx(1)-1,length(time_idx),timestrd,''double'');'])
340 341
                    end
                case 'nele'
342
                    eval([varnames{aa},'=netcdf.getVar(nc,varID,''double'');'])
343 344 345 346
                    if nele_idx>=0
                        eval([varnames{aa},' = ',varnames{aa},'(nele_idx);'])
                    end
                case 'node'
347
                    eval([varnames{aa},'=netcdf.getVar(nc,varID,''double'');'])
348 349 350 351
                    if node_idx>=0
                        eval([varnames{aa},' = ',varnames{aa},'(node_idx);'])
                    end
                otherwise
352 353 354
                    if ftbverbose
                        fprintf('Unkown dimension for variable %s. Skipping to next one in function call.\n', name);
                    end
355 356 357 358
            end
        otherwise
            % identified dimensions to restrict
            do_restrict=zeros(size(dimName));
359
            dimidx=nan(size(dimName));
360
            clear start count stride
361
            for dd=1:length(dimName)
362 363
                start.(dimName{dd})=[];
                count.(dimName{dd})=[];
364
                stride.(dimName{dd})=[];
365
                test=find(strcmpi(RestrictDims.Name,dimName{dd}));
366
                if ~isempty(test); dimidx(dd)=test; end
367 368 369
            end
            % create start index for dimensions of the variable to
            % access
370
            if any(isfinite(dimidx))
371
                % we have at least two valid dimension indices, proceed
372
                for dd=1:length(dimidx)
373 374
                    % restrict time as range but node and nele dims are
                    % considered as stations rather than ranges.
375 376
                    % if restriction is not -1 then select specified
                    % indices otherwise read all
377
                    if ~isnan(dimidx(dd)) && RestrictDims.idx{dimidx(dd)}(1)>=0
378
                        if (strcmpi(dimName(dd),'time'))
379
                            start.(dimName{dd})=RestrictDims.idx{dimidx(dd)}(1)-1;
380 381
                            count.(dimName{dd})=length(RestrictDims.idx{dimidx(dd)});
                            stride.(dimName{dd})=timestrd;
Pierre Cazenave's avatar
Pierre Cazenave committed
382

383
                        else
384 385 386
                            for ss=1:length(RestrictDims.idx{dimidx(dd)})
                                start.(dimName{dd})(ss)=RestrictDims.idx{dimidx(dd)}(ss)-1;
                                count.(dimName{dd})(ss)=1;
387
                                stride.(dimName{dd})=1;
388
                            end
389
                        end
390
                        do_restrict(dd)=1;
391 392 393
                    else
                        start.(dimName{dd})=0;
                        count.(dimName{dd})=dimLength(dd);
394
                        stride.(dimName{dd})=1;
395 396 397
                    end
                end
            else
398 399 400
                if ftbverbose
                    fprintf('Wrong selection of dimensions to extract.\nExtracting all values in current variable.\n');
                end
401
            end
402 403
            %
            %             eval([varnames{aa},'=netcdf.getVar(nc,varID,start,count,''double'');'])
404
            cc_names=fieldnames(count);
405
            clear read_start read_count read_stride
406 407 408 409 410 411 412 413
            switch sum(do_restrict) % there are dimensions to restrict
                case 1 % only one dimension to restrict
                    switch find(do_restrict) % find position of restrictive variable
                        case 1 % restrict the first variable
                            % but the variable can have more than 2 dimensions
                            switch dimens
                                % initialise variable
                                case 2
414
                                    rr=[min(sum(count.(cc_names{1})),dimLength(1)) min(sum(count.(cc_names{2})),dimLength(2))];
415 416 417
                                case 3
                                    rr=[min(sum(count.(cc_names{1})),dimLength(1)),...
                                        min(sum(count.(cc_names{2})),dimLength(2)),...
418
                                        min(sum(count.(cc_names{3})),dimLength(3))];
419
                            end
420

421
                            eval([varnames{aa},'=nan(rr);'])
422 423 424
                            % reorganize start and count arrays
                            read_start(find(~do_restrict))=start.(cc_names{find(~do_restrict)});
                            read_count(find(~do_restrict))=count.(cc_names{find(~do_restrict)});
425
                            read_stride(find(~do_restrict))=stride.(cc_names{find(~do_restrict)});
Pierre Cazenave's avatar
Pierre Cazenave committed
426

427 428 429
                            for cc=1:length(start.(cc_names{find(do_restrict)}))
                                read_start(find(do_restrict))=start.(cc_names{find(do_restrict)})(cc);
                                read_count(find(do_restrict))=count.(cc_names{find(do_restrict)})(cc);
430
                                read_stride(find(do_restrict))=stride.(cc_names{find(do_restrict)});
Pierre Cazenave's avatar
Pierre Cazenave committed
431

432
                                var_dump=netcdf.getVar(nc,varID,read_start,read_count,read_stride,'double');
Pierre Cazenave's avatar
Pierre Cazenave committed
433

434 435 436 437 438 439 440 441
                                eval([varnames{aa},'(cc,:)=var_dump;'])
                                clear var_dump
                            end
                        case 2 % restrict the second variable (ie depth)
                            % but the variable can have more than 2 dimensions
                            switch dimens
                                % initialise variable
                                case 2
442
                                    rr=[min(sum(count.(cc_names{1})),dimLength(1)) min(sum(count.(cc_names{2})),dimLength(2))];
443 444 445
                                case 3
                                    rr=[min(sum(count.(cc_names{1})),dimLength(1)),...
                                        min(sum(count.(cc_names{2})),dimLength(2)),...
446
                                        min(sum(count.(cc_names{3})),dimLength(3))];
447
                            end
448

449
                            eval([varnames{aa},'=nan(rr);'])
450 451 452
                            % reorganize start and count arrays
                            read_start(find(~do_restrict))=start.(cc_names{find(~do_restrict)});
                            read_count(find(~do_restrict))=count.(cc_names{find(~do_restrict)});
453
                            read_stride(find(~do_restrict))=stride.(cc_names{find(~do_restrict)});
454

455
                            for cc=1:length(start.(cc_names{logical(do_restrict)}))
456 457
                                read_start(find(do_restrict))=start.(cc_names{find(do_restrict)})(cc);
                                read_count(find(do_restrict))=count.(cc_names{find(do_restrict)})(cc);
458 459
                                read_stride(find(do_restrict))=stride.(cc_names{find(do_restrict)});
                                var_dump=netcdf.getVar(nc,varID,read_start,read_count,read_stride,'double');
460 461 462 463
                                try
                                    eval([varnames{aa},'(:,cc)=var_dump;'])
                                catch
                                    eval([varnames{aa},'(:,:,cc)=var_dump;'])
464

465
                                end
466 467 468 469 470 471
                                clear var_dump
                            end
                        case 3 % restrict the second variable (ie depth)
                            % but the variable needs to have at least 3 dimensions
                            rr=[min(sum(count.(cc_names{1})),dimLength(1)),...
                                min(sum(count.(cc_names{2})),dimLength(2)),...
472
                                min(sum(count.(cc_names{3})),dimLength(3))];
473

474
                            eval([varnames{aa},'=nan(rr);'])
475 476 477 478 479
                            % reorganize start and count arrays
                            % There are now 2 unrestricted dimensions
                            for tt=find(~do_restrict)
                                read_start(tt)=start.(cc_names{tt});
                                read_count(tt)=count.(cc_names{tt});
480
                                read_stride(tt)=stride.(cc_names{tt});
Pierre Cazenave's avatar
Pierre Cazenave committed
481

482
                            end
483

484 485 486 487 488 489
                            % check if time is one of them
                            if ~isempty(find(dimidx==5))
                                do_time = find(dimidx==5); % 5 is the index for time
                                % reorganize start and count arrays
                                read_start(do_time)=start.(cc_names{do_time});
                                read_count(do_time)=count.(cc_names{do_time});
490 491
                                read_stride(do_time)=stride.(cc_names{do_time});
                                eval([varnames{aa},'=netcdf.getVar(nc,varID,read_start,read_count,read_stride,''double'');'])
492 493 494 495
                            else % we are looking at stations or depth layers
                                for cc=1:length(start.(cc_names{(do_restrict)}))
                                    read_start(find(do_restrict))=start.(cc_names{find(do_restrict)})(cc);
                                    read_count(find(do_restrict))=count.(cc_names{find(do_restrict)})(cc);
496 497
                                    read_stride(find(do_restrict))=stride.(cc_names{find(do_restrict)});
                                    var_dump=netcdf.getVar(nc,varID,read_start,read_count,read_stride,'double');
Pierre Cazenave's avatar
Pierre Cazenave committed
498

499 500 501 502 503 504 505 506
                                    switch dimName(find(do_restrict))
                                        case 'node' | 'nele'
                                            eval([varnames{aa},'(cc,:,:)=var_dump;'])
                                        case 'siglay' | 'siglev'
                                            eval([varnames{aa},'(:,cc,:)=var_dump;'])
                                    end
                                    clear var_dump
                                end
507

508 509 510 511
                            end
                    end
                    eval(['selection.',varnames{aa},'.start=start;'])
                    eval(['selection.',varnames{aa},'.count=count;'])
512

513 514 515 516 517
                case 2 % Two dimension to restrict!
                    % but the variable can have more than 2 dimensions
                    switch dimens
                        % initialise variable
                        case 2
518
                            rr=[min(sum(count.(cc_names{1})),dimLength(1)) min(sum(count.(cc_names{2})),dimLength(2))];
519 520 521
                        case 3
                            rr=[min(sum(count.(cc_names{1})),dimLength(1)),...
                                min(sum(count.(cc_names{2})),dimLength(2)),...
522
                                min(sum(count.(cc_names{3})),dimLength(3))];
523
                    end
524

525
                    eval([varnames{aa},'=nan(rr);'])
526 527 528 529 530 531
                    % check if time is one of them
                    if ~isempty(find(dimidx==5))
                        do_time = find(dimidx==5); % 5 is the index for time
                        % reorganize start and count arrays
                        read_start(do_time)=start.(cc_names{do_time});
                        read_count(do_time)=count.(cc_names{do_time});
532
                        read_stride(do_time)=stride.(cc_names{do_time});
533
                        % search for the non_restrictive variable
534 535 536 537
                        %                         cc=1
                        %                         while ~(length( start.(cc_names{cc}))==1);cc=cc+1;end
                        % esto esta mal.... tengo que incluir otra opcion por si tenemos una
                        % variable de dos dimensiones donde los dos son restrictivas....
538 539
                        cc=find(~do_restrict);
                        if isempty(cc);cc=length(cc_names);end
540 541
                        read_start(cc)=start.(cc_names{cc});
                        read_count(cc)=count.(cc_names{cc});
542
                        read_stride(cc)=stride.(cc_names{cc});
543 544
                        do_other = setdiff(dimidx,[dimidx(cc),5]) ; % one of these is also restrictive...
                        do_other=find(dimidx==do_other);
545

546 547 548
                        for cc=1:length(start.(cc_names{do_other}))
                            read_start(do_other)=start.(cc_names{do_other})(cc);
                            read_count(do_other)=count.(cc_names{do_other})(cc);
549 550
                            read_stride(do_other)=stride.(cc_names{do_other});
                            var_dump=netcdf.getVar(nc,varID,read_start,read_count,read_stride,'double');
551 552 553 554 555 556 557 558 559
                            switch do_other
                                case 1
                                    eval([varnames{aa},'(cc,:,:)=var_dump;'])
                                case 2
                                    eval([varnames{aa},'(:,cc,:)=var_dump;'])
                                case 3
                                    eval([varnames{aa},'(:,:,cc)=var_dump;'])
                            end
                            clear var_dump
560
                        end
561 562 563
                    else % time is not one of them so we need to restrict both variables...
                        % in this case it doesn't really matter
                        % which one we restrict firts...
564

565 566 567 568
                        for kk=1:length(start.(cc_names{1}))
                            % reorganize start and count arrays
                            read_start(1)=start.(cc_names{1})(kk);
                            read_count(1)=count.(cc_names{1})(kk);
569
                            read_stride(1)=stride.(cc_names{1});
570 571 572
                            for cc=1:length(start.(cc_names{2}))
                                read_start(2)=start.(cc_names{2})(cc);
                                read_count(2)=count.(cc_names{2})(cc);
573 574
                                read_stride(2)=stride.(cc_names{2});
                                var_dump=netcdf.getVar(nc,varID,read_start,read_count,read_stride,'double');
Pierre Cazenave's avatar
Pierre Cazenave committed
575

576 577 578
                                eval([varnames{aa},'(kk,cc)=var_dump;'])
                                clear var_dump
                            end
579

580
                        end
581

582
                    end
583

584 585
                    eval(['selection.',varnames{aa},'.start=start;'])
                    eval(['selection.',varnames{aa},'.count=count;'])
Pierre Cazenave's avatar
Pierre Cazenave committed
586

587 588 589 590 591 592 593 594 595 596 597
                case 3 % three dimension to restrict!
                    % but the variable can have more than 2 dimensions
                    switch dimens
                        % initialise variable
                        case 2
                            rr=[min(sum(count.(cc_names{1})),dimLength(1)) min(sum(count.(cc_names{2})),dimLength(2))];
                        case 3
                            rr=[min(sum(count.(cc_names{1})),dimLength(1)),...
                                min(sum(count.(cc_names{2})),dimLength(2)),...
                                min(sum(count.(cc_names{3})),dimLength(3))];
                    end
Pierre Cazenave's avatar
Pierre Cazenave committed
598

599 600 601 602 603 604 605 606
                    eval([varnames{aa},'=nan(rr);'])
                    % check if time is one of them
                    if isempty(find(dimidx==5));disp('This won''t work, try again');return;end
                    do_time = find(dimidx==5); % 5 is the index for time
                    % reorganize start and count arrays
                    read_start(do_time)=start.(cc_names{do_time});
                    read_count(do_time)=count.(cc_names{do_time});
                    read_stride(do_time)=stride.(cc_names{do_time});
Pierre Cazenave's avatar
Pierre Cazenave committed
607

608 609 610 611 612 613 614 615 616 617 618
                    % search for the non_restrictive variable
                    %                         cc=1
                    %                         while ~(length( start.(cc_names{cc}))==1);cc=cc+1;end
                    % esto esta mal.... tengo que incluir otra opcion por si tenemos una
                    % variable de dos dimensiones donde los dos son restrictivas....
                    [~,do_other] = setdiff(dimidx,[dimidx(do_time)]) ; % these are also restrictive and are not time...
                    if length(count.(cc_names{do_other(1)})) <length(count.(cc_names{do_other(2)}))
                        do_one=do_other(1);do_two=do_other(2);
                    else
                        do_one=do_other(2);do_two=do_other(1);
                    end
Pierre Cazenave's avatar
Pierre Cazenave committed
619

620 621 622 623
                    for cc=1:length(start.(cc_names{do_one}))
                        read_start(do_one)=start.(cc_names{do_one})(cc);
                        read_count(do_one)=count.(cc_names{do_one})(cc);
                        read_stride(do_one)=stride.(cc_names{do_one});
Pierre Cazenave's avatar
Pierre Cazenave committed
624

625 626 627 628
                        for pp=1:length(start.(cc_names{do_two}))
                            read_start(do_two)=start.(cc_names{do_two})(pp);
                            read_count(do_two)=count.(cc_names{do_two})(pp);
                            read_stride(do_two)=stride.(cc_names{do_two});
Pierre Cazenave's avatar
Pierre Cazenave committed
629

630 631 632 633 634 635 636 637
                            var_dump=netcdf.getVar(nc,varID,read_start,read_count,read_stride,'double');
                            eval([varnames{aa},'(pp,cc,:)=var_dump;'])
                        end
                        clear var_dump
                    end
                    eval(['selection.',varnames{aa},'.start=start;'])
                    eval(['selection.',varnames{aa},'.count=count;'])
                case 0 % there are NO dimensions to restrict and 3 dimensions haven't been coded yet!!
Pierre Cazenave's avatar
Pierre Cazenave committed
638

639 640 641
                    for nn=1:length(cc_names)
                        read_start(nn)=start.(cc_names{nn});
                        read_count(nn)=count.(cc_names{nn});
642
                         read_stride(nn)=stride.(cc_names{nn});
643
                    end
644
                    eval([varnames{aa},'=netcdf.getVar(nc,varID,read_start,read_count,read_stride,''double'');'])
645 646
                    eval(['selection.',varnames{aa},'.start=start;'])
                    eval(['selection.',varnames{aa},'.count=count;'])
647

648 649
            end
    end
650
    eval(['data.(varnames{aa}) = ',varnames{aa},';'])
651
    eval(['clear ',varnames{aa}])
652
%     toc
653
end
654 655

%--------------------------------------------------------------------------
656
% Tidy up, finish and return data
657
%--------------------------------------------------------------------------
658 659

netcdf.close(nc)
660 661 662

if ftbverbose
    fprintf('end   : %s \n', subname)
Pierre Cazenave's avatar
Pierre Cazenave committed
663
end