Commit eee55d6e authored by Rafael Laboissiere's avatar Rafael Laboissiere

Update upstream source from tag 'upstream/1.0.4'

Update to upstream version '1.0.4'
with Debian dir 0e6f0902150e36dd1e926aef92e5fee40d77c7b6
parents ddd02460 ec22a8f1
/debian/files
/debian/octave-ncarray.debhelper.log
/debian/octave-ncarray.substvars
/debian/octave-ncarray/
/global-list
/local-list
Name: ncarray
Version: 1.0.3
Date: 2014-05-04
Version: 1.0.4
Date: 2016-11-29
Author: Alexander Barth <barth.alexander@gmail.com>
Maintainer: Alexander Barth <barth.alexander@gmail.com>
Title: ncarray
......
Summary of important user-visible changes for ncarray 1.0.4:
------------------------------------------------------------
** Avoid obsolete strmatch function
** Use a tolerance in test script (bug #49391)
Summary of important user-visible changes for ncarray 1.0.3:
------------------------------------------------------------
......
......@@ -32,7 +32,7 @@ end
function x = funelem(x,nm)
% make sure x is not an ncArray
x = full(x);
mask = isnan(x) || isnan(nm);
mask = isnan(x) | isnan(nm);
diff = zeros(size(x));
diff(mask) = x(mask) - nm(mask);
x = {diff.^2, ~mask};
......
......@@ -2,17 +2,20 @@ function B = subsref(self,idx)
% handle case with a single subscript
% for example CA(234)
if strcmp(idx.type,'()') && length(idx.subs) == 1
% indices of elements in B
ind = idx.subs{1};
if strcmp(ind,':')
ind = [1:numel(self)]';
end
% output array
B = zeros(size(ind));
B(:) = NaN;
if self.overlap
% transform index to subscript
% transform linear index ind to subscript subs
subs = cell(1,self.nd);
[subs{:}] = ind2sub(size(self),ind);
......@@ -26,22 +29,7 @@ if strcmp(idx.type,'()') && length(idx.subs) == 1
for i=1:length(ind)
idxe.subs = mat2cell(subs(i,:),1,ones(1,self.nd));
%B(i) = subsref_canonical(self,idxe);
[idx_global,idx_local,sz] = idx_global_local_(self,idxe);
tmp = zeros(sz);
tmp(:) = NaN;
for j=1:self.na
% get subset from j-th array
subset = subsref(self.arrays{j},idx_local{j});
% set subset in global array B
tmp = subsasgn(tmp,idx_global{j},subset);
end
B(i) = tmp;
B(i) = subsref_canonical(self,idxe);
end
else
% assume all arrays does not overlap
......
......@@ -10,8 +10,13 @@ fprintf('Coordinates:\n')
for i = 1:length(c)
tmp = sprintf('%dx',size(c(i).val));
fprintf(' Name: "%s" standard name: "%s" size %s\n',...
c(i).name,c(i).standard_name,tmp(1:end-1));
stdname = c(i).standard_name;
if isempty(stdname)
stdname = '<unset>';
end
fprintf(' Name: %15s; standard name: %25s; size %10s\n',...
c(i).name,stdname,tmp(1:end-1));
end
......
......@@ -2,6 +2,7 @@
% Load a subset of a variable based on range of coordiante variables.
% The names of the coordinates (coord_name1, coord_name2,...) coorespond to the standard_name attribute.
% Only 1-dimensional coordinates are currently supported.
% Time units are converted to "datenum".
%
%
% Example
......@@ -15,6 +16,18 @@ for i = 1:length(c)
c(i).v = full(c(i).val);
% per default take all data along a dimension
c(i).index = ':';
% convert time units
if ~isempty(strfind(c(i).units,'since'))
[t0,f] = nctimeunits(c(i).units);
c(i).v = f*double(c(i).v) + t0;
end
% change vertical axis to positive up
if strcmp(c(i).positive,'down')
c(i).v = -double(c(i).v);
end
c(i).sub = c(i).v;
end
......@@ -38,7 +51,11 @@ for i = 1:2:length(varargin)
if numel(range) == 1
dist = abs(c(j).v - range);
[mindist,i] = min(dist);
%i
%mindist
%c(j).v(i)
%datevec(c(j).v(i))
else
i = find(range(1) < c(j).v & c(j).v < range(end));
i = min(i):max(i);
......@@ -49,6 +66,7 @@ for i = 1:2:length(varargin)
end
idx = substruct('()',{c.index});
%idx
data = subsref (self,idx);
varargout = {data,c.sub};
......
......@@ -38,7 +38,7 @@
% Web: http://modb.oce.ulg.ac.be/mediawiki/index.php/ncArray
% hidded constructor signature:
% data = ncArray(filename,varname)
% data = ncArray(var,dims,coord);
% is used to create data with coordinate values by ncCatArray
function retval = ncArray(varargin)
......
......@@ -33,7 +33,6 @@ if length(idx) == 2 && strcmp(idx(2).type,'.') && strcmp(idx(2).subs,'coord')
j = sort(j);
idx_c.type = '()';
idx_c.subs = idx(1).subs(j);
varargout{i} = subsref(self.coord(i).val,idx_c);
end
else
......
......@@ -37,14 +37,19 @@ else
if length(tmp) == 1
start(i) = tmp;
else
% check if indexes are at regular intervals
test = tmp(1):tmp(2)-tmp(1):tmp(end);
if size(tmp,2) == 1
% tmp is a row vector, make test also a row vector
test = test';
end
if all(tmp == test)
start(i) = tmp(1);
stride(i) = tmp(2)-tmp(1);
count(i) = (tmp(end)-tmp(1))/stride(i) +1;
else
error('indeces');
error('indeces not at regular intervals');
end
end
end
......
......@@ -52,7 +52,10 @@ if strcmp(idx.type,'()')
elseif strcmp(idx.type,'.')
% load attribute
name = idx.subs;
index = strmatch(name,{self.vinfo.Attributes(:).Name});
% strmatch is obsolete
% index = strmatch(name,{self.vinfo.Attributes(:).Name});
index = find(strcmp(name,{self.vinfo.Attributes(:).Name}));
if isempty(index)
error('variable %s has no attribute called %s',self.varname,name);
......
% Decompress a file using a cache.
%
% [fname,success]=cached_decompress(filename)
% [fname,success] = cached_decompress(filename)
%
% Input:
% filename: name of the file which is possibly compressed
......@@ -14,15 +14,21 @@
% CACHED_DECOMPRESS_LOG_FID (default 1): file id for log message
% CACHED_DECOMPRESS_MAX_SIZE (default 1e10): maximum size of cache in bytes.
% Alexander Barth, 2012-06-13
%
function [fname]=cached_decompress(url)
function fname = cached_decompress(url)
global CACHED_DECOMPRESS_DIR
global CACHED_DECOMPRESS_LOG_FID
global CACHED_DECOMPRESS_MAX_SIZE
if startswith(url,'http:') || ...
~(endswith(url,'.gz') || endswith(url,'.bz2') || endswith(url,'.xz'))
% opendap url or not compressed file
fname = url;
return
end
cache_dir = CACHED_DECOMPRESS_DIR;
if isempty(cache_dir)
......@@ -32,13 +38,6 @@ if isempty(cache_dir)
fprintf('creating directory %s for temporary files.\n',cache_dir);
end
if beginswith(url,'http:') || ...
~(endswith(url,'.gz') || endswith(url,'.bz2') || endswith(url,'.xz'))
% opendap url or not compressed file
fname = url;
return
end
%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 '...
......@@ -47,10 +46,10 @@ end
% where to print logs? default to screen
fid=CACHED_DECOMPRESS_LOG_FID;
fid = CACHED_DECOMPRESS_LOG_FID;
if (isempty(fid))
fid=1;
fid = 1;
end
% form filename for cache
......@@ -58,22 +57,35 @@ end
fname = url;
fname = strrep(fname,'/','_SLASH_');
fname = strrep(fname,'*','_STAR_');
fname = strrep(fname,'\','_BSLASH_');
fname = strrep(fname,':','_COLON_');
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);
cmd = 'gunzip --stdout -';
elseif endswith(url,'.bz2')
syscmd('bunzip2 --stdout "%s" > "%s"',url,fname);
cmd = 'bunzip2 --stdout -';
elseif endswith(url,'.xz')
cmd = 'unxz --stdout -';
else
cmd = 'cat';
end
if startswith(url,'ftp://')
syscmd('curl --silent "%s" | %s > "%s"',url,cmd,fname);
else
syscmd('unxz --stdout "%s" > "%s"',url,fname);
syscmd('%s < "%s" > "%s"',cmd,url,fname);
end
else
% fprintf(fid,'retrieve from cache %s\n',url);
end
% check cache size
d=dir(cache_dir);
......@@ -107,23 +119,23 @@ if (cashe_size > max_cache_size)
end
end
function t = beginswith(s,pre)
if length(pre) <= length(s)
t = strcmp(s(1:length(pre)),pre);
else
function t = startswith(s,ext)
if length(ext) <= length(s)
t = strcmp(s(1:length(ext)),ext);
else
t = 0;
end
end
end
function t = endswith(s,ext)
if length(ext) <= length(s)
if length(ext) <= length(s)
t = strcmp(s(end-length(ext)+1:end),ext);
else
else
t = 0;
end
end
end
......@@ -141,7 +153,7 @@ end
end
% Copyright (C) 2012-2013 Alexander Barth <barth.alexander@gmail.com>
% Copyright (C) 2012-2013, 2015 Alexander Barth <barth.alexander@gmail.com>
%
% 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
......
......@@ -80,6 +80,11 @@ elseif ischar(pattern)
'http://www.mathworks.com/matlabcentral/fileexchange/19550']);
end
end
if isempty(filenames)
error('ncArray:nomatch','no file found matching %s',pattern);
end
elseif isa(pattern, 'function_handle')
filenames = cell(1,length(range));
......@@ -120,12 +125,24 @@ if (dim > length(dims)) && ~isempty(coord)
dims{dim} = catdimname;
coord(dim).dims = {catdimname};
coord(dim).val = range;
coord(dim).name = catdimname;
end
for i=1:length(coord)
%test if value is already defined, if yes do nothing
if isfield(coord(i),'val')
if ~isempty(coord(i).val)
continue
end
end
% the number of the dimension might be different
% find in coord(i).dims the index of the dimension called dims{dim}
% for example we concatenate over time, then two situations can arrise
% the coordinate variable lon can dependent on time (dimc is not empty)
% or it does not depdent on time (dimc is empty)
dimc = find(strcmp(coord(i).dims,dims{dim}));
if isempty(dimc)
......@@ -133,12 +150,13 @@ for i=1:length(coord)
coord(i).val = ncBaseArray(filenames{1},coord(i).name,'vinfo',vinfo);
else
% coordinates do also depend on the dimension over which we concatenate
%i,coord(i).name,dimc,dims{dim}
coord(i).val = arr(dimc,filenames,coord(i).name,finfos);
end
if dim > length(coord(i).dims)
coord(i).dims{dim} = catdimname;
end
%if dim > length(coord(i).dims)
% coord(i).dims{dim} = catdimname;
%end
end
data = ncArray(var,dims,coord);
......
% Tutorial for using ncArray
% It is advised to run this script in an empty directory.
% It will delete and overwrite files named file1.nc, file2.nc and file3.nc.
% size of the example data (2x3)
n = 3;
m = 2;
% create 3 files (file1.nc, file2.nc,...) with a 2x3 variable called SST
data = zeros(n,m);
disp('create example files: file1.nc, file2.nc, file3.nc')
for i = 1:3
data(:) = i;
files{i} = sprintf('file%d.nc',i);
delete(files{i});
ncarray_example_file(files{i},data);
end
% Using ncArray
SST = ncArray('file1.nc','SST');
disp('load the entire file')
data = SST(:,:,:);
disp('get the attribute units')
units = SST.units;
disp('load a particular value');
data = SST(3,2,1);
% Using ncCatArray
disp('concatenate the files over the 3rd dimension (here time)')
SST = ncCatArray(3,{'file1.nc','file2.nc','file3.nc'},'SST');
% or just
% SST = ncCatArray(3,'file*.nc','SST');
disp('load all 3 files');
data = SST(:,:,:);
disp('load a particular value (1,2,1) of the 3rd file');
data = SST(1,2,3);
function ncarray_example_file(filename,data)
% create an example NetCDF file with the name filename and given data
if ~isempty(which('nccreate')) && ~isempty(which('ncwriteatt')) && ...
~isempty(which('ncwrite'))
% use matlab netcdf high level interface
function ncarray_example_file(filename,data)
%dtype = 'single';
dtype = 'double';
sz = size(data);
% Variables
nccreate(filename,'lon','Format','classic','Datatype',dtype,...
'Dimensions',{'x',220, 'y',144});
'Dimensions',{'x',sz(1), 'y',sz(2)});
ncwriteatt(filename,'lon','long_name','Longitude')
ncwriteatt(filename,'lon','units','degrees_east')
nccreate(filename,'lat','Datatype',dtype,'Dimensions',{'x',220, 'y',144});
nccreate(filename,'lat','Datatype',dtype,'Dimensions',{'x',sz(1), 'y',sz(2)});
ncwriteatt(filename,'lat','long_name','Latitude')
ncwriteatt(filename,'lat','units','degrees_north')
......@@ -22,7 +20,7 @@ if ~isempty(which('nccreate')) && ~isempty(which('ncwriteatt')) && ...
ncwriteatt(filename,'time','units','days since 1858-11-17 00:00:00 GMT')
nccreate(filename,'SST','Datatype',dtype,'Dimensions',...
{'x',220, 'y',144, 'time',1});
{'x',sz(1), 'y',sz(2), 'time',1});
ncwriteatt(filename,'SST','missing_value',single(9999))
ncwriteatt(filename,'SST','_FillValue',single(9999))
ncwriteatt(filename,'SST','units','degC')
......@@ -30,56 +28,19 @@ if ~isempty(which('nccreate')) && ~isempty(which('ncwriteatt')) && ...
ncwriteatt(filename,'SST','coordinates','lat lon')
ncwrite(filename,'SST',data);
else
% use octcdf interface
nc = netcdf(filename,'c');
% dimensions
nc('x') = size(data,1);
nc('y') = size(data,2);
nc('time') = size(data,3);
% variables
nc{'lon'} = ncdouble('y','x'); % 31680 elements
nc{'lon'}.long_name = ncchar('Longitude');
nc{'lon'}.units = ncchar('degrees_east');
nc{'lat'} = ncdouble('y','x'); % 31680 elements
nc{'lat'}.long_name = ncchar('Latitude');
nc{'lat'}.units = ncchar('degrees_north');
nc{'time'} = ncdouble('time'); % 1 elements
nc{'time'}.long_name = ncchar('Time');
nc{'time'}.units = ncchar('days since 1858-11-17 00:00:00 GMT');
nc{'SST'} = ncdouble('time','y','x'); % 31680 elements
nc{'SST'}.missing_value = ncdouble(9999);
nc{'SST'}.FillValue_ = ncdouble(9999);
nc{'SST'}.units = ncchar('degC');
nc{'SST'}.long_name = ncchar('Sea Surface Temperature');
nc{'SST'}.coordinates = ncchar('lat lon');
% global attributes
nc{'SST'}(:) = permute(data,[3 2 1]);
close(nc)
end
% Copyright (C) 2012,2013 Alexander Barth <barth.alexander@gmail.com>
%
% 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/>.
% Copyright (C) 2012,2013,2015 Alexander Barth <barth.alexander@gmail.com>
%
% 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/>.
......@@ -12,6 +12,9 @@
% coord is an empty structure if no coordinate information have been
% found.
% TODO: use a predictable order for coord:
% lon, lat, depth, time, ensemble,...
% Author: Alexander Barth (barth.alexander@gmail.com)
function [dims,coord] = nccoord(filename,varname)
......@@ -32,7 +35,8 @@ vinfo = finfo.Variables(index);
dims = {vinfo.Dimensions(:).Name};
% create empty coord array with the fields name and dims
coord = struct('name',{},'dims',{},'standard_name',{},'units',{});
coord = struct('name',{},'dims',{},'standard_name',{},...
'units',{},'positive',{});
% check the coordinate attribute
if ~isempty(vinfo.Attributes)
......@@ -73,6 +77,7 @@ if isempty(find(strcmp(name,{coord(:).name}),1))
c.dims = {d(:).Name};
c.standard_name = [];
c.units = [];
c.positive = [];
% get standard_name attribute if present
i = find(strcmp('standard_name',{vinfo.Attributes(:).Name}));
......@@ -86,6 +91,12 @@ if isempty(find(strcmp(name,{coord(:).name}),1))
c.units = vinfo.Attributes(i).Value;
end
% get positive attribute if present
i = find(strcmp('positive',{vinfo.Attributes(:).Name}));
if ~isempty(i)
c.positive = vinfo.Attributes(i).Value;
end
coord(end+1) = c;
end
end
......@@ -93,7 +104,7 @@ end
end
% Copyright (C) 2012, 2013 Alexander Barth <barth.alexander@gmail.com>
% Copyright (C) 2012, 2013, 2015 Alexander Barth <barth.alexander@gmail.com>
%
% 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
......
......@@ -7,7 +7,7 @@
function t = ncreadtime(filename,varname)
t = ncread(filename,varname);
t = double(ncread(filename,varname));
units = ncreadatt(filename,varname,'units');
[t0,f] = nctimeunits(units);
......
......@@ -5,31 +5,65 @@
% Parse the netCDF time unit u and returns the time offset (days since 31
% December 1 BC, as datenum) and scaling factor f (in days).
% See the netCDF CF convention for the structure of the time units.
% http://cfconventions.org/Data/cf-conventions/cf-conventions-1.6/build/cf-conventions.html#time-coordinate
% Also: http://www.unidata.ucar.edu/software/thredds/current/netcdf-java/CDM/CalendarDateTime.html
function [t0,f] = nctimeunits(u)
if strfind(u,'seconds')
% years in days for udunits
% http://cfconventions.org/Data/cf-conventions/cf-conventions-1.6/build/cf-conventions.html#time-coordinate
year_in_days = 365.242198781;
l = strfind(u,'since');
if length(l) ~= 1
error(['time units sould expect one "since": "' u '"']);
end
period = strtrim(lower(u(1:l-1)));
reference_date = strtrim(u(l+6:end));
if strcmp(period,'millisec') || strcmp(period,'msec')
f = 1/(24*60*60*1000);
elseif strcmp(period,'second') || strcmp(period,'seconds') ...
|| strcmp(period,'s') || strcmp(period,'sec')
f = 1/(24*60*60);
elseif strfind(u,'hours')
elseif strcmp(period,'minute') || strcmp(period,'minutes') ...
|| strcmp(period,'min')
f = 1/(24*60);
elseif strcmp(period,'hour') || strcmp(period,'hours') ...
|| strcmp(period,'hr')
f = 1/24;
elseif strfind(u,'days')
elseif strcmp(period,'day') || strcmp(period,'days')
f = 1;
elseif strcmp(period,'week') || strcmp(period,'weeks')
f = 1/(24*60*60*7);
elseif strcmp(period,'year') || strcmp(period,'years') ...
strcmp(period,'yr')
f = year_in_days;
elseif strcmp(period,'month') || strcmp(period,'months') ...
strcmp(period,'mon')
f = year_in_days/12;
else
error(['unknown units "' u '"']);
error(['unknown units "' period '"']);
end
l = strfind(u,'since')+6;
try
t0 = datenum(u(l:end),'yyyy-mm-dd HH:MM:SS');
catch
if strcmp(reference_date,'1900-01-01 00:00:0.0')
t0 = datenum(1900,1,1);
else
try
t0 = datenum(u(l:end),'yyyy-mm-ddTHH:MM:SS');
t0 = datenum(reference_date,'yyyy-mm-dd HH:MM:SS');
catch
try
t0 = datenum(u(l:end),'yyyy-mm-dd');
t0 = datenum(reference_date,'yyyy-mm-ddTHH:MM:SS');
catch
error(['date format is not recogized ' u(l:end)])
try
t0 = datenum(reference_date,'yyyy-mm-dd');
catch
error(['date format is not recogized ' reference_date])
end
end
end
end
% this function is necessary because of the limitation of matlab
function assertAlmostEqual(observed,expected)
% tolerance for testing
tol = 1e-10;
% for compatibility with matlab which does not have
% assert (OBSERVED, EXPECTED, TOL)
assert(max(abs(observed(:) - expected(:))) < tol)
% for octave prior to 3.8.0
if isempty(which('isequaln'))
isequaln = @(x,y) isequalwithequalnans(x,y);
end
% check also NaNS
assert(isequal(isnan(observed),isnan(expected)))
end
This diff is collapsed.
......@@ -24,11 +24,11 @@ end
data = ncCatArray(3,files,varname);
reddata = nanmean(data,3);
reddataref = nanmean(dataref,3);
assert(isequaln(reddata, reddataref))
assertAlmostEqual(reddata, reddataref)
reddata = nansum(data,3);
reddataref = nansum(dataref,3);
assert(isequaln(reddata, reddataref))
assertAlmostEqual(reddata, reddataref)
reddata = nanvar(data,[],3);
reddataref = nanvar(dataref,[],3);
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment