Skip to content

Commit 0d96587

Browse files
cwhansecwhanse
authored
Add pvl_robustfit (#4)
* Add pvl_robustfit Added pvl_robustfit.m to provide an option for pvl_desoto_parameter_estimation and pvl_PVsyst_parameter_estimation that doesn’t rely on the Matlab Statistics Toolbox. Added py_rlm.py which is the Python code used if the Statistics Toolbox is not present. * Huld model parameter estimation Add new function * Correction to pvl_snlinverter change ACPower(ACPower<Ps0)= -1.*abs(Pnt); to ACPower(DCPower<Ps0)= -1.*abs(Pnt); * Clear sky detection algorithm and example data Initial upload * Add code for equally spaced data Execution is much faster with GHI data at equally spaced time stamps * Add pvl_detect_shadows; update help library pvl_detect_shadows identifies shadows from nearby objects such as wires and poles in a time series of GHI data.
1 parent cf226ec commit 0d96587

File tree

152 files changed

+2840
-113
lines changed

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

152 files changed

+2840
-113
lines changed

‎Example Data/PSEL_GHI_2012.mat

1.91 MB
Binary file not shown.

‎Example Data/SMUD22_GHI_data.mat

2.22 MB
Binary file not shown.

‎GHIImage.m

Lines changed: 387 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,387 @@
1+
classdef GHIImage
2+
properties
3+
Value
4+
end
5+
6+
methods
7+
function obj = GHIImage(val)
8+
if nargin>0
9+
if ismatrix(val) && ~isvector(val) &&(isnumeric(val) | islogical(val))
10+
obj.Value = val;
11+
else
12+
error('GHIIMage must be 2D numeric or logical array');
13+
end
14+
end
15+
end
16+
17+
function matrixOut = smooth2a(matrixIn,Nr,Nc)
18+
% Smooths 2D array data. Ignores NaN's.
19+
%
20+
% function matrixOut = smooth2a(matrixIn,Nr,Nc)
21+
%
22+
% This function smooths the data in matrixIn using a mean filter over a
23+
% rectangle of size (2*Nr+1)-by-(2*Nc+1). Basically, you end up replacing
24+
% element "i" by the mean of the rectange centered on "i". Any NaN
25+
% elements are ignored in the averaging. If element "i" is a NaN, then it
26+
% will be preserved as NaN in the output. At the edges of the matrix,
27+
% where you cannot build a full rectangle, as much of the rectangle that
28+
% fits on your matrix is used (similar to the default on Matlab's builtin
29+
% function "smooth").
30+
%
31+
% "matrixIn": original matrix
32+
% "Nr": number of points used to smooth rows
33+
% "Nc": number of points to smooth columns. If not specified, Nc = Nr.
34+
%
35+
% "matrixOut": smoothed version of original matrix
36+
%
37+
%
38+
% Written by Greg Reeves, March 2009.
39+
% Division of Biology
40+
% Caltech
41+
%
42+
% Inspired by "smooth2", written by Kelly Hilands, October 2004
43+
% Applied Research Laboratory
44+
% Penn State University
45+
%
46+
% Developed from code written by Olof Liungman, 1997
47+
% Dept. of Oceanography, Earth Sciences Centre
48+
% G�teborg University, Sweden
49+
% E-mail: [email protected]
50+
51+
%
52+
% Initial error statements and definitions
53+
%
54+
if nargin < 2, error('Not enough input arguments!'), end
55+
56+
N(1) = Nr;
57+
if nargin < 3, N(2) = N(1); else N(2) = Nc; end
58+
59+
if length(N(1)) ~= 1, error('Nr must be a scalar!'), end
60+
if length(N(2)) ~= 1, error('Nc must be a scalar!'), end
61+
62+
%
63+
% Building matrices that will compute running sums. The left-matrix, eL,
64+
% smooths along the rows. The right-matrix, eR, smooths along the
65+
% columns. You end up replacing element "i" by the mean of a (2*Nr+1)-by-
66+
% (2*Nc+1) rectangle centered on element "i".
67+
%
68+
[row,col] = size(matrixIn.Value);
69+
eL = spdiags(ones(row,2*N(1)+1),(-N(1):N(1)),row,row);
70+
eR = spdiags(ones(col,2*N(2)+1),(-N(2):N(2)),col,col);
71+
72+
%
73+
% Setting all "NaN" elements of "matrixIn" to zero so that these will not
74+
% affect the summation. (If this isn't done, any sum that includes a NaN
75+
% will also become NaN.)
76+
%
77+
A = isnan(matrixIn.Value);
78+
matrixIn.Value(A) = 0;
79+
80+
%
81+
% For each element, we have to count how many non-NaN elements went into
82+
% the sums. This is so we can divide by that number to get a mean. We use
83+
% the same matrices to do this (ie, "eL" and "eR").
84+
%
85+
nrmlize = eL*(~A)*eR;
86+
nrmlize(A) = NaN;
87+
88+
%
89+
% Actually taking the mean.
90+
%
91+
matrixOut = eL*matrixIn.Value*eR;
92+
matrixOut = matrixOut./nrmlize;
93+
end
94+
95+
function out = bwareaopen(X,connectedness,threshold)
96+
% BWAREAOPEN: perform area opening of binary image X
97+
%
98+
% Out = BWAREAOPEN(X, Y,SE)
99+
%
100+
% Inputs:
101+
% X: input image array (assumed binary [0 1])
102+
%
103+
% Outputs
104+
% Out: output image array
105+
if ~islogical(X.Value)
106+
error('bwareaopen passed non-binary array')
107+
end
108+
109+
cc = findCC(X,connectedness); % Get connected components
110+
counts = histcounts(cc,'BinMethod','integers');
111+
% integers = 0:65000;
112+
% counts = hist(cc(:),integers);
113+
out = ismember(cc,find(counts(2:end)>threshold));
114+
end
115+
116+
function out = bwclose(img, SE)
117+
% BWCLOSE: morphological close of binary image with structuring element SE
118+
%
119+
% Out = BWCLOSE(In,SE)
120+
%
121+
% Inputs:
122+
% In: input image array (assumed binary [0 1])
123+
% SE: structuring element (assumed binary)
124+
%
125+
% Outputs
126+
% Out: output image array
127+
if ~islogical(img.Value)
128+
error('bwclose passed non-binary array')
129+
end
130+
131+
% Version 1: use convolution
132+
tmp = img.bwdilate(SE);
133+
out = GHIImage(tmp).bwerode(SE);
134+
end
135+
136+
function out = bwdilate(img, SE)
137+
% BWDILATE: Dilate a binary image with structuring element SE
138+
%
139+
% Out = BWDILATE(In,SE)
140+
%
141+
% Inputs:
142+
% In: input image array (assumed binary [0 1])
143+
% SE: structuring element (assumed binary)
144+
%
145+
% Outputs
146+
% Out: output image array
147+
if ~islogical(img.Value)
148+
error('bwdilate passed non-binary array')
149+
end
150+
151+
% Version 1: use convolution
152+
out = conv2(single(img.Value), single(SE),'same') > 0;
153+
end
154+
155+
function out = bwerode(img, SE)
156+
% BWERODE: Erode a binary image with structuring element SE
157+
%
158+
% Out = BWERODE(In,SE)
159+
%
160+
% Inputs:
161+
% In: input image array (assumed binary [0 1])
162+
% SE: structuring element (assumed binary)
163+
%
164+
% Outputs
165+
% Out: output image array
166+
if ~islogical(img.Value)
167+
error('bwerode passed non-binary array')
168+
end
169+
170+
% Version 2: erosion of foreground is dilation of background
171+
out = ~(conv2(single(~img.Value), single(SE), 'same') > 0);
172+
173+
% Version 1: use convolution
174+
% out = conv2(single(img),single(SE),'same') == conv2(ones(size(img)),SE,'same');
175+
% Previously, I used divide, which probably takes longer than array ==
176+
% out = conv2(single(img),single(SE),'same')./conv2(ones(size(img)),SE,'same') == 1;
177+
178+
end
179+
180+
function out = bwhitmiss(X,J,K)
181+
% BWHITMISS: perform hit or miss operation on binary image X
182+
%
183+
% BWHITMISS(X,J) uses the same structuring element in both phases of
184+
% the hit-miss operation: BWHITMISS(X,J) = BWERODE(X,J) & BWERODE(~X,~J)
185+
%
186+
% BWHITMISS(X,J,K) accommodates "don't care" by allowing different
187+
% structuring elements: BWHITMISS(X,J,K) = BWERODE(X,J) & BWERODE(~X,K)
188+
% In this form, K should be the complement of J, with "don't cares" reset
189+
% to 0.
190+
%
191+
% Usage: Out = BWHITMISS(X,J) or Out = BWHITMISS(X,J,K)
192+
%
193+
% Inputs:
194+
% X: input image array (assumed logical)
195+
% J, K: binary structuring elements
196+
%
197+
% Outputs
198+
% Out: output (binary/logical) image array
199+
200+
if ~islogical(X.Value)
201+
error('bwhitmiss passed non-binary array')
202+
end
203+
204+
if nargin < 3
205+
K = ~J;
206+
end
207+
notX = GHIImage(~X.Value);
208+
out = X.bwerode(J) & notX.bwerode(K);
209+
end
210+
211+
function out = bwlengthopen(img,connectedness,threshold)
212+
% BWAREAOPEN: eliminate connected components in img shorter than threshold
213+
%
214+
% Out = BWAREAOPEN(img,Connectedness,MinLength)
215+
%
216+
% Inputs:
217+
% img: input image array (assumed binary [0 1])
218+
% Connectedness: 4 or 8
219+
% MinLength: minimum length of connected component to retain
220+
%
221+
% Outputs
222+
% Out: output image array
223+
224+
if ~islogical(img.Value)
225+
error('bwlengthopen passed non-binary array')
226+
end
227+
228+
if ~ismember([4 8], connectedness)
229+
error('bwlengthopen: connectedness must be 4 or 8')
230+
end
231+
232+
cc = findCC(img,connectedness); % Get connected components
233+
234+
for ii = 1:max(cc(:))
235+
[~,cols] = ind2sub(size(cc),find(cc == ii));
236+
len(ii) = max(cols)-min(cols)+1;
237+
end
238+
239+
out = ismember(cc,find(len>threshold));
240+
end
241+
242+
function out = bwopen(img, SE)
243+
% BWOPEN: morphological open of binary image with structuring element SE
244+
%
245+
% Out = BWOPEN(In,SE)
246+
%
247+
% Inputs:
248+
% In: input image array (assumed binary [0 1])
249+
% SE: structuring element (assumed binary)
250+
%
251+
% Outputs
252+
% Out: output image array
253+
254+
% Version 1: use convolution
255+
tmp = img.bwerode(SE);
256+
257+
out = GHIImage(tmp).bwdilate(SE);
258+
end
259+
260+
function out = graydilate(img, SE)
261+
262+
% img: grayscale image
263+
% SE: structuring element (assumed binary)
264+
265+
[M,N] = size(img.Value);
266+
[SEM, SEN] = size(SE);
267+
268+
halfSEM = floor(SEM/2);
269+
halfSEN = floor(SEN/2);
270+
271+
bufM = M+2*halfSEM;
272+
bufN = N+2*halfSEN;
273+
274+
% Create vector of offsets
275+
offsets = bsxfun(@plus, bufM*(-halfSEN:halfSEN), (-halfSEM:halfSEM)');
276+
offsets = offsets(SE>0);
277+
278+
% Put input image into a buffer (-inf guarantees nonselection by max())
279+
buf = -inf*ones(bufM, bufN);
280+
buf(halfSEM+(1:M),halfSEN+(1:N)) = img.Value;
281+
282+
% Iterate through original image pixels
283+
out = zeros(M, N);
284+
index = offsets + halfSEN * bufM + halfSEM;
285+
for jj = 1:N
286+
for ii = 1:M
287+
index = index + 1;
288+
out(ii,jj) = max(buf(index));
289+
end
290+
% Skip to beginning of next column of original image
291+
index = index + 2*halfSEM;
292+
end
293+
end
294+
295+
function out = grayerode(img, SE)
296+
297+
% img: grayscale image
298+
% SE: structuring element (assumed binary)
299+
300+
[M,N] = size(img.Value);
301+
[SEM, SEN] = size(SE);
302+
303+
halfSEM = floor(SEM/2);
304+
halfSEN = floor(SEN/2);
305+
306+
bufM = M+2*halfSEM;
307+
bufN = N+2*halfSEN;
308+
309+
% Create vector of offsets
310+
offsets = bsxfun(@plus, bufM*(-halfSEN:halfSEN), (-halfSEM:halfSEM)');
311+
offsets = offsets(SE>0);
312+
313+
% Put input image into a buffer (inf guarantees nonselection by min())
314+
buf = inf*ones(bufM, bufN);
315+
buf(halfSEM+(1:M),halfSEN+(1:N)) = img.Value;
316+
317+
% Iterate through original image pixels by columns (reduces multiplies inside loop)
318+
out = zeros(M, N);
319+
index = offsets + halfSEN * bufM + halfSEM;
320+
for jj = 1:N
321+
for ii = 1:M
322+
index = index + 1;
323+
out(ii,jj) = min(buf(index));
324+
end
325+
% Skip to beginning of next column of original image
326+
index = index + 2*halfSEM;
327+
end
328+
end
329+
330+
function out = graygradient(img, SE)
331+
332+
% img: grayscale image
333+
% SE: structuring element (assumed binary)
334+
335+
out = graydilate(img,SE) - grayerode(img,SE);
336+
337+
end
338+
339+
function Connected = findCC(img,connectivity)
340+
%OnePass: One-pass connected component labeling algorithm
341+
% See Wikipedia page "Connected-component labeling"
342+
343+
% Put a buffer of zeros around the input image
344+
[M,N] = size(img.Value);
345+
buf = zeros(M+2,N+2);
346+
buf(2:M+1,2:N+1) = img.Value;
347+
348+
% Initialize
349+
[M,N] = size(buf);
350+
Connected = zeros(M,N);
351+
Mark = 1; % "Value"
352+
Difference = 1; % "Increment"
353+
% Index = [];
354+
Nobj = 0;
355+
if connectivity == 4
356+
Offsets = [-M; -1; M; 1];
357+
elseif connectivity == 8
358+
Offsets = [-M+(-1:1) -1 1 M+(-1:1)]';
359+
else
360+
error('Connectivity should be specified as 4 or 8');
361+
end
362+
363+
% Iterate across rows of original image pixels
364+
for ii = 2:M-1
365+
for jj = 2:N-1
366+
if buf(ii,jj)==1
367+
Nobj = Nobj + 1;
368+
Index = ((jj-1)*M + ii);
369+
Connected(Index) = Mark;
370+
while ~isempty(Index)
371+
buf(Index) = 0;
372+
Neighbors = bsxfun(@plus,Index,Offsets');
373+
Neighbors = unique(Neighbors(:));
374+
Index = Neighbors(buf(Neighbors)==1);
375+
Connected(Index) = Mark;
376+
end
377+
Mark = Mark + Difference;
378+
end
379+
end
380+
end
381+
Connected = Connected(2:M-1,2:N-1);
382+
end
383+
384+
end
385+
end
386+
387+

0 commit comments

Comments
 (0)