Menu
▾
▴
[Octave-cvsupdate] SF.net SVN: octave:[10663] trunk/octave-forge/extra
|
From: <aba...@us...> - 2012-06-22 15:04:20
|
Revision: 10663
http://octave.svn.sourceforge.net/octave/?rev=10663&view=rev
Author: abarth93
Date: 2012-06-22 15:04:09 +0000 (Fri, 22 Jun 2012)
Log Message:
-----------
initial import of ncArray
Added Paths:
-----------
trunk/octave-forge/extra/ncArray/
trunk/octave-forge/extra/ncArray/DESCRIPTION
trunk/octave-forge/extra/ncArray/inst/
trunk/octave-forge/extra/ncArray/inst/cached_decompress.m
trunk/octave-forge/extra/ncArray/inst/ncCatArray.m
trunk/octave-forge/extra/ncArray/inst/ncCatData.m
trunk/octave-forge/extra/ncArray/inst/nccoord.m
trunk/octave-forge/extra/ncArray/inst/test_ncarray.m
Added: trunk/octave-forge/extra/ncArray/DESCRIPTION
===================================================================
--- trunk/octave-forge/extra/ncArray/DESCRIPTION (rev 0)
+++ trunk/octave-forge/extra/ncArray/DESCRIPTION 2012-06-22 15:04:09 UTC (rev 10663)
@@ -0,0 +1,11 @@
+Name: ncArray
+Version: 1.0.0
+Date: 2012-06-22
+Author: Alexander Barth <bar...@gm...>
+Maintainer: Alexander Barth <bar...@gm...>
+Title: ncArray
+Description: Subsetting and concatenating of NetCDF files
+Depends: octave (>= 3.4.0)
+Autoload: yes
+License: GPL version 2 or later
+Url: http://octave.sf.net
Added: trunk/octave-forge/extra/ncArray/inst/cached_decompress.m
===================================================================
--- trunk/octave-forge/extra/ncArray/inst/cached_decompress.m (rev 0)
+++ trunk/octave-forge/extra/ncArray/inst/cached_decompress.m 2012-06-22 15:04:09 UTC (rev 10663)
@@ -0,0 +1,118 @@
+% [fname,success]=cached_decompress(filename)
+%
+% Decompress a file if it is not already in cache
+%
+% Input:
+% filename: name of the file which is possibly compressed
+%
+% Output:
+% fname: the filename of the uncompressed file
+%
+% Example:
+%
+%
+% Alexander Barth, 2012-06-13
+%
+
+
+function [fname]=cached_decompress(url)
+global CASHED_GUNZIP_DIR
+global CASHED_GUNZIP_LOG_FID
+
+
+cache_dir = CASHED_GUNZIP_DIR;
+if isempty(cache_dir)
+% cache_dir = fullfile(getenv('HOME'),'tmp','Cache');
+ cache_dir = fullfile(getenv('HOME'),'tmp','Cache');
+end
+
+
+if endswith(url,'.gz') || endswith(url,'.bz2')
+ if exist(cache_dir,'dir') ~= 7
+ error(['cache directory for compressed files does not exist. '...
+ 'Please create the directory %s or change le value of the '...
+ 'global variable CASHED_GUNZIP_DIR'],cache_dir);
+ end
+else
+ fname = url;
+ return
+end
+
+% where to print logs? default to screen
+
+fid=CASHED_GUNZIP_LOG_FID;
+
+if (isempty(fid))
+ fid=1;
+end
+
+% form filename for cache
+
+fname = url;
+fname = strrep(fname,'/','_SLASH_');
+fname = strrep(fname,'*','_STAR_');
+fname = fullfile(cache_dir,fname);
+
+% test if in cache
+
+if exist(fname,'file') ~= 2
+ if endswith(url,'.gz')
+ syscmd('gunzip --stdout "%s" > "%s"',url,fname);
+ else
+ syscmd('bunzip2 --stdout "%s" > "%s"',url,fname);
+ end
+else
+% fprintf(fid,'retrieve from cache %s\n',url);
+end
+
+% check cache size
+
+d=dir(cache_dir);
+cashe_size = sum([d.bytes]);
+max_cache_size = 1e10;
+
+if (cashe_size > max_cache_size)
+
+ % look for oldest files
+ fdate = zeros(1,length(d));
+ for i=1:length(d);
+ fdate(i) = datenum(d(i).date);
+ end
+
+ [fdate,index] = sort(fdate,'descend');
+ d=d(index);
+
+ cum_size = cumsum([d(:).bytes]);
+ todelete = find(cum_size > max_cache_size);
+
+ for i=todelete
+ if (d(i).isdir == 0)
+ fprintf(fid,'clean cashe: delete %s\n', d(i).name);
+ delete(fullfile(cache_dir,d(i).name));
+ end
+ end
+end
+end
+
+function t = endswith(s,ext)
+
+if length(ext) <= length(s)
+ t = strcmp(s(end-length(ext)+1:end),ext);
+else
+ t = 0;
+end
+end
+
+
+function syscmd(varargin)
+
+cmd = sprintf(varargin{:});
+%disp(cmd)
+status = 0;
+[status, output] = system(cmd);
+
+disp(output);
+if status ~= 0
+ error(['command "' cmd '" failed: ' output]);
+end
+end
\ No newline at end of file
Added: trunk/octave-forge/extra/ncArray/inst/ncCatArray.m
===================================================================
--- trunk/octave-forge/extra/ncArray/inst/ncCatArray.m (rev 0)
+++ trunk/octave-forge/extra/ncArray/inst/ncCatArray.m 2012-06-22 15:04:09 UTC (rev 10663)
@@ -0,0 +1,81 @@
+% C = ncCatArray(dim,filenames,varname)
+% C = ncCatArray(dim,pattern,varname)
+% C = ncCatArray(dim,filenamefun,varname,range)
+%
+% create a concatenated array from variables (varname) in a list of
+% netcdf files along dimension dim.Individual elements can be accessed by
+% subscribs, e.g. C(2,3) and the corrsponding subset of the appropriate file is loaded
+%
+% This list of netcdf files can be specified as a cell array (filenames),
+% shell wildcard pattern (e.g. file_*.nc) or a function handle
+% filenamefun. In this later case, this i-th filename is
+% filenamefun(range(i)).
+%
+% Example:
+%
+% data = ncCatArray(3,{'file-20120708.nc','file-20120709.nc'},'SST')
+%
+% data = ncCatArray(3,'file-*.nc','SST')
+%
+% data = ncCatArray(3,@(t) ['file-' datestr(t,'yyyymmdd') '.nc'],...
+% datenum(2012,07,08):datenum(2012,07,09));
+%
+% Note: in Octave the glob function is used to determine files matching the
+% shell wildcard pattern, while in Matlab rdir is used. The function rdir
+% is available from Matlab exchange under BSD license
+% (http://www.mathworks.com/matlabcentral/fileexchange/19550).
+
+% Author: Alexander Barth (bar...@gm...)
+%
+function ncCA = ncCatArray(dim,pattern,varname,range)
+
+if iscell(pattern)
+ filenames = pattern;
+
+elseif ischar(pattern)
+ try
+ filenames = glob(pattern);
+ catch
+ try
+ d = rdir(pattern);
+ filenames = {d(:).name};
+ catch
+ error(['The function rdir or glob (octave) is not available. '...
+ 'rdir can be installed from '...
+ 'http://www.mathworks.com/matlabcentral/fileexchange/19550']);
+ end
+ end
+elseif isa(pattern, 'function_handle')
+ filenames = cell(1,length(range));
+
+ for i=1:length(range)
+ filenames{i} = pattern(range(i));
+ end
+end
+
+arrays = cell(1,length(filenames));
+
+for i=1:length(filenames)
+ arrays{i} = ncArray(filenames{i},varname);
+end
+
+
+ncCA = CatArray(dim,arrays);
+
+% Copyright (C) 2012 Alexander Barth <bar...@gm...>
+%
+% This program is free software; you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation; either version 2 of the License, or
+% (at your option) any later version.
+%
+% This program is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with this program; If not, see <http://www.gnu.org/licenses/>.
+
+
+
Added: trunk/octave-forge/extra/ncArray/inst/ncCatData.m
===================================================================
--- trunk/octave-forge/extra/ncArray/inst/ncCatData.m (rev 0)
+++ trunk/octave-forge/extra/ncArray/inst/ncCatData.m 2012-06-22 15:04:09 UTC (rev 10663)
@@ -0,0 +1,111 @@
+% C = ncCatArray(dim,filenames,varname)
+% C = ncCatArray(dim,pattern,varname)
+% C = ncCatArray(dim,filenamefun,varname,range)
+%
+% create a concatenated array from variables (varname) in a list of
+% netcdf files along dimension dim.Individual elements can be accessed by
+% subscribs, e.g. C(2,3) and the corrsponding subset of the appropriate file is loaded
+
+% This list of netcdf files can be specified as a cell array (filenames),
+% shell wildcard pattern (e.g. file_*.nc) or a function handle
+% filenamefun. In this later case, this i-th filename is
+% filenamefun(range(i)).
+%
+% Example:
+%
+% data = ncCatArray(3,{'file-20120708.nc','file-20120709.nc'},'SST')
+%
+% data = ncCatArray(3,'file-*.nc','SST')
+%
+% data = ncCatArray(3,@(t) ['file-' datestr(t,'yyyymmdd') '.nc'],...
+% datenum(2012,07,08):datenum(2012,07,09));
+%
+% Note: in Octave the glob function is used to determine files matching the
+% shell wildcard pattern, while in Matlab rdir is used. The function rdir
+% is available from Matlab exchange under BSD license
+% (http://www.mathworks.com/matlabcentral/fileexchange/19550).
+
+% Author: Alexander Barth (bar...@gm...)
+%
+function data = ncCatData(dim,pattern,varname,range)
+
+catdimname = '_cat_dim';
+
+if iscell(pattern)
+ filenames = pattern;
+
+elseif ischar(pattern)
+ try
+ filenames = glob(pattern);
+ catch
+ try
+ d = rdir(pattern);
+ filenames = {d(:).name};
+ catch
+ error(['The function rdir or glob (octave) is not available. '...
+ 'rdir can be installed from '...
+ 'http://www.mathworks.com/matlabcentral/fileexchange/19550']);
+ end
+ end
+elseif isa(pattern, 'function_handle')
+ filenames = cell(1,length(range));
+
+ for i=1:length(range)
+ filenames{i} = pattern(range(i));
+ end
+end
+
+if nargin == 3
+ range = 1:length(filenames);
+end
+
+var = arr(dim,filenames,varname);
+
+[dims,coord] = nccoord(cached_decompress(filenames{1}),varname);
+
+if dim > length(dims)
+ % concatenate is new dimension
+ dims{dim} = catdimname;
+ coord(dim).dims = {catdimname};
+ coord(dim).val = range;
+end
+
+
+for i=1:length(coord)
+ % coordinates do also depend on the dimension only which we concatenate
+ coord(i).val = arr(dim,filenames,coord(i).name);
+ if dim > length(coord(i).dims)
+ coord(i).dims{dim} = catdimname;
+ end
+end
+
+data = ncData(var,dims,coord);
+
+end
+
+
+function CA = arr(dim,filenames,varname)
+arrays = cell(1,length(filenames));
+for i=1:length(filenames)
+ arrays{i} = ncArray(filenames{i},varname);
+end
+
+CA = CatArray(dim,arrays);
+end
+% Copyright (C) 2012 Alexander Barth <bar...@gm...>
+%
+% This program is free software; you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation; either version 2 of the License, or
+% (at your option) any later version.
+%
+% This program is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with this program; If not, see <http://www.gnu.org/licenses/>.
+
+
+
Added: trunk/octave-forge/extra/ncArray/inst/nccoord.m
===================================================================
--- trunk/octave-forge/extra/ncArray/inst/nccoord.m (rev 0)
+++ trunk/octave-forge/extra/ncArray/inst/nccoord.m 2012-06-22 15:04:09 UTC (rev 10663)
@@ -0,0 +1,54 @@
+% coord = nccoord(filename,varname)
+% get coordinates of a variables using CF convention
+
+function [dims,coord] = nccoord(filename,varname)
+
+finfo = ncinfo(filename);
+vinfo = ncinfo(filename,varname);
+
+% determine coordinates
+% using CF convention
+
+dims = vinfo.Dimensions;
+
+% create empty coord array with the fields name and dims
+coord = struct('name',{},'dims',{});
+
+% check the coordinate attribute
+index = strmatch('coordinates',{vinfo.Attributes(:).Name});
+if ~isempty(index)
+ tmp = strsplit(vinfo.Attributes(index).Value,' ');
+
+ for i=1:length(tmp)
+ coord = addcoord(coord,tmp{i},finfo);
+ end
+end
+
+% check for coordinate dimensions
+for i=1:length(dims)
+ % check if variable with the same name than the dimension exist
+ index = strmatch(dims{i},{finfo.Variables(:).Name});
+ if ~isempty(index)
+ coord = addcoord(coord,dims{i},finfo);
+ end
+end
+
+
+end
+
+function coord = addcoord(coord,name,finfo)
+
+% check if coordinate is aleady in the list
+if isempty(strmatch(name,{coord(:).name}))
+
+ % check if name is variable
+ index = strmatch(name,{finfo.Variables(:).Name});
+ if ~isempty(index)
+ c.name = name;
+ c.dims = finfo.Variables(index).Dimensions;
+
+ coord(end+1) = c;
+ end
+end
+
+end
\ No newline at end of file
Added: trunk/octave-forge/extra/ncArray/inst/test_ncarray.m
===================================================================
--- trunk/octave-forge/extra/ncArray/inst/test_ncarray.m (rev 0)
+++ trunk/octave-forge/extra/ncArray/inst/test_ncarray.m 2012-06-22 15:04:09 UTC (rev 10663)
@@ -0,0 +1,233 @@
+function test_ncarray()
+% test ncArray, ncCatArray, ncData and ncCatData
+
+varname = 'SST';
+
+tmpdir = tempname;
+mkdir(tmpdir);
+
+tmpfname = tempname(tmpdir);
+for i = 1:3
+ files{i} = fullfile(tmpdir,sprintf('file%d.nc',i));
+ ncarray_example_file(files{i},randn(220,144));
+end
+
+
+filename = files{1};
+
+
+% test ncread/ncwrite
+
+copyfile(files{1},tmpfname);
+SST_ref = ncread(files{1},'SST');
+ncwrite(tmpfname,'SST',zeros(size(SST_ref)));
+test = ncread(tmpfname,'SST');
+
+assert(all(test(:) == 0))
+
+ncwrite(tmpfname,'SST',SST_ref);
+test = ncread(tmpfname,'SST');
+assert(isequalwithequalnans(test,SST_ref))
+
+
+%%% test ncArray
+
+% reading
+
+copyfile(files{2},tmpfname);
+SST = ncArray(tmpfname,varname);
+test = SST(:,:,:);
+SST_ref = ncread(tmpfname,varname);
+
+assert(isequalwithequalnans(test,SST_ref))
+assert(isempty(SST(:,:,:,[])));
+assert(isequalwithequalnans(SST_ref, SST(:,:,:,1)))
+
+ind = floor(numel(SST_ref) * rand(100,1))+1;
+assert(isequalwithequalnans(SST(ind),SST_ref(ind)))
+
+% writing
+
+r = round(randn(size(SST)));
+SST(:,:,:) = r;
+SST_ref = ncread(tmpfname,varname);
+assert(isequalwithequalnans(r,SST_ref));
+
+SST(:,:,:) = 3 * r;
+SST_ref = ncread(tmpfname,varname);
+assert(isequalwithequalnans(3 * r,SST_ref));
+
+
+%%% CatArray
+
+% reading
+
+CA = CatArray(3,{...
+ ncArray(filename,varname),...
+ ncArray(files{2},varname),...
+ ncArray(files{3},varname)...
+ });
+
+assert(isequalwithequalnans(size(CA),[220 144 3]))
+
+SST_ref = ncread(filename,'SST');
+tmp2 = CA(:,:,1);
+assert(isequalwithequalnans(SST_ref,tmp2))
+
+
+
+SST_test = CA(:,:,2);
+SST_ref = ncread(files{2},'SST');
+
+
+
+assert(isequalwithequalnans(SST_test,SST_ref))
+
+CA2 = CatArray(4,{...
+ ncArray(files{1},varname),...
+ ncArray(files{2},varname),...
+ ncArray(files{3},varname)...
+ });
+
+SST_test = CA2(:,:,:,2);
+
+assert(isequalwithequalnans(SST_test,SST_ref))
+
+CA2 = ncCatArray(3,{...
+ files{1},...
+ files{2},...
+ files{3}},...
+ varname);
+
+SST_test = CA2(:,:,2);
+assert(isequalwithequalnans(SST_test,SST_ref))
+
+CA2 = ncCatArray(3,fullfile(tmpdir,'file*nc'),varname);
+SST_test = CA2(:,:,2);
+assert(isequalwithequalnans(SST_test,SST_ref))
+
+
+CA2 = ncCatArray(3,...
+ @(i) fullfile(tmpdir,sprintf('file%d.nc',i)),...
+ varname,...
+ 1:3);
+
+SST_test = CA2(:,:,2);
+assert(isequalwithequalnans(SST_test,SST_ref))
+
+SST_ref = cat(3,...
+ ncread(files{1},'SST'),...
+ ncread(files{2},'SST'),...
+ ncread(files{3},'SST'));
+
+
+assert(isequalwithequalnans(CA2(:,:,:),SST_ref))
+
+assert(isequalwithequalnans(CA2(:,:,1),SST_ref(:,:,1)))
+assert(isequalwithequalnans(CA2(3:5:50,3:5:100,1),SST_ref(3:5:50,3:5:100,1)))
+assert(isequalwithequalnans(CA2(3:5:50,3:5:100,2),SST_ref(3:5:50,3:5:100,2)))
+assert(isequalwithequalnans(CA2(3:5:50,3:5:100,3),SST_ref(3:5:50,3:5:100,3)))
+assert(isequalwithequalnans(CA2(3:5:50,3:5:100,end),SST_ref(3:5:50,3:5:100,end)))
+assert(isequalwithequalnans(CA2(50,100,1:3),SST_ref(50,100,1:3)))
+assert(isequalwithequalnans(CA2(3:5:50,3:5:100,1:2:3),SST_ref(3:5:50,3:5:100,1:2:3)))
+assert(isequalwithequalnans(CA2(3:5:50,3:5:end,1:2:3),SST_ref(3:5:50,3:5:end,1:2:3)))
+assert(isequalwithequalnans(CA2(3:5:50,3:5:end,:),SST_ref(3:5:50,3:5:end,:)))
+ind = floor(numel(SST_ref) * rand(100,1))+1;
+assert(isequalwithequalnans(CA2(ind),SST_ref(ind)))
+
+% writing
+
+for i=1:3
+ list{i} = tempname;
+ copyfile(filename,list{i});
+end
+
+CA2 = ncCatArray(3,list,varname);
+r = round(randn(size(CA2)));
+CA2(:,:,:) = r;
+
+check = ncread(list{2},varname);
+assert(isequalwithequalnans(check,r(:,:,2)))
+
+r2 = round(randn(size(CA2)));
+r(3:5:50,3:5:end,:) = r2(3:5:50,3:5:end,:);
+CA2(3:5:50,3:5:end,:) = r2(3:5:50,3:5:end,:);
+assert(isequalwithequalnans(CA2(:,:,:),r))
+
+r(end-1,3:5:end,1:2:3) = 2*r2(end-1,3:5:end,1:2:3);
+CA2(end-1,3:5:end,1:2:3) = 2*r2(end-1,3:5:end,1:2:3);
+assert(isequalwithequalnans(CA2(:,:,:),r))
+
+
+
+
+if 1
+ % test ncData (constructor: ncData(var,dims,coord)
+
+ SST = ncArray(filename,varname);
+ SST_ref = ncread(filename,varname);
+ lon_ref = ncread(filename,'lon');
+
+ coord(1).val = ncArray(filename,'lon');
+ coord(1).dims = {'x','y'};
+
+ coord(2).val = ncArray(filename,'lat');
+ coord(2).dims = {'x','y'};
+
+ coord(3).val = ncArray(filename,'time');
+ coord(3).dims = {'time'};
+
+ data = ncData(SST,{'x','y','time'},coord);
+
+ [x,y,t] = data(:,:,:).coord;
+
+ assert(isequalwithequalnans(data(:,:,:),SST_ref))
+ assert(isequalwithequalnans(x,lon_ref))
+
+ assert(isequalwithequalnans(data(),SST_ref))
+ [x,y,t] = data().coord;
+ assert(isequalwithequalnans(x,lon_ref))
+
+ assert(isequalwithequalnans(data(1:3:end,:,:),SST_ref(1:3:end,:,:)))
+ [x,y,t] = data(1:3:end,:,:).coord;
+ assert(isequalwithequalnans(x,lon_ref(1:3:end,:)))
+
+ % test ncData (constructor: ncData(filename,varname)
+ SST = ncData(filename,varname);
+ [x,y,t] = data(:,:,:).coord;
+
+ assert(isequalwithequalnans(data(:,:,:),SST_ref))
+ assert(isequalwithequalnans(x,lon_ref))
+
+ assert(isequalwithequalnans(data(),SST_ref))
+ [x,y,t] = data().coord;
+ assert(isequalwithequalnans(x,lon_ref))
+
+ assert(isequalwithequalnans(data(1:3:end,:,:),SST_ref(1:3:end,:,:)))
+ [x,y,t] = data(1:3:end,:,:).coord;
+ assert(isequalwithequalnans(x,lon_ref(1:3:end,:)))
+
+
+ assert(strcmp(SST.units,'degC'))
+ assert(strcmp(SST.('units'),'degC'))
+
+end
+
+
+% read compressed data
+zname = [tmpfname '.gz'];
+system(['gzip --stdout ' tmpfname ' > ' zname]);
+
+SST = ncData(zname,'SST');
+SST_ref = ncread(tmpfname,'SST');
+assert(isequalwithequalnans(SST(),SST_ref))
+
+
+
+CA2 = ncCatData(3,fullfile(tmpdir,'file*nc'),varname);
+SST_test = CA2(:,:,2);
+SST_ref = ncread(files{2},'SST');
+assert(isequalwithequalnans(SST_test,SST_ref))
+
+disp('All tests passed.')
+
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
