-
Notifications
You must be signed in to change notification settings - Fork 110
/
aggregate.m
95 lines (84 loc) · 2.57 KB
/
aggregate.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
function C = aggregate(DEMlowres,DEMhighres,aggfun)
%AGGREGATE resampling a GRIDobj using aggregation
%
% Syntax
%
% C = aggregate(B,A)
% C = aggregate(B,xy)
% C = aggregate(B,xyz)
% C = aggregate(...,aggfun)
%
% Description
%
% This function resamples the grid A to C to match the extent and
% resolution of grid B. B must spatially overlap with A and must have a
% coarser resolution. By default, aggregate uses the mean to calculate
% new values in grid C, but aggregate takes any other function that takes
% a vector and returns a scalar (e.g. median, std, ...).
%
% Values to be aggregated can also be supplied as list of coordinates
% (and attributes). This is particularly useful if point density for
% each pixel is greater than one.
%
% Input arguments
%
% A high resolution GRIDobj
% B low resolution GRIDobj. A and B must have the same
% coordinate system
% xy two column matrix with coordinates. aggfun is @numel.
% xyz three column matrix with coordinates and third column being
% some attribute (e.g. elevation, ...)
% aggfun anonymous function that aggregate values. The function must
% take a vector and return a scalar. The default is @mean.
% Other possible functions are @numel to obtain counts,
% @median, @std, ...
%
% Output arguments
%
% C GRIDobj with same extent and resolution as B and aggregated
% values derived from A
%
%
% See also: GRIDobj/resample, accumarray
%
% Author: Wolfgang Schwanghart (w.schwanghart[at]geo.uni-potsdam.de)
% Date: 12. July, 2018
if nargin == 2
aggfun = @mean;
end
if isa(DEMhighres,'GRIDobj')
% A is a high res GRIDobj
[x,y] = getcoordinates(DEMhighres);
IX = zeros(DEMhighres.size);
[X,Y] = getcoordinates(DEMlowres);
sy = size(y);
for r = 1:numel(x)
IX(:,r) = coord2ind(X,Y,repmat(x(r),sy),y);
end
Z = DEMhighres.Z(:);
IX = IX(:);
inan = isnan(IX);
IX(inan) = [];
Z(inan) = [];
else
% A is a list of coordinates and attributes
xyz = DEMhighres;
IX = coord2ind(DEMlowres,xyz(:,1),xyz(:,2));
inan = isnan(IX);
xyz(inan,:) = [];
IX(inan) = [];
if size(xyz,2) == 3
Z = xyz(:,end);
else
Z = true; %true(size(IX));
aggfun = @numel;
end
end
if isempty(Z)
C = GRIDobj(DEMlowres);
return
end
z = accumarray(IX,Z,[prod(DEMlowres.size) 1],aggfun);
C = DEMlowres;
C.Z = reshape(z,DEMlowres.size);
% C.Z(C.Z == 0) = nan;