-
Notifications
You must be signed in to change notification settings - Fork 0
/
escape_nc_search.m
58 lines (47 loc) · 2.15 KB
/
escape_nc_search.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
function [indices, timestamps] = escape_nc_search(ncfile,LWC_threshold)
% Search for cloud passes based on liquid water content.
% Outputs list of indices for the cloud pass, as well as the Start and
% End time of segments in UTC
% Setting the duration threshold (in seconds) will limit cloud passes to sections
% that are above the lwc threshold for the specified duration or
% longer.
% Define LWC threshold for cloud
% LWC_threshold = 0.005; % g/m3
duration_threshold = 10; % seconds
% Get data from the netCDF file
time = ncread(ncfile,'Time');
cdplwc = ncread(ncfile,'lwc_cdp_sp_rt'); % cdp liquid water content
nevzlwc = ncread(ncfile, 'lwc_nevz_sp_rt'); % nevzerov liquid water content
nevztwc = ncread(ncfile, 'twc_nevz_sp_rt'); % nevzerov total water content
% Which lwc value to use for the search
LWC = nevzlwc;
% Find logical vector where lwc > threshold
binaryVector = LWC > LWC_threshold;
% Label each region with a label - an "ID" number.
[labeledVector, numRegions] = bwlabel(binaryVector);
% Measure lengths of each region and the indexes
measurements = regionprops(labeledVector, LWC, 'Area', 'PixelValues', 'PixelIdxList');
% Find regions where the area (length) are 3 or greater and
% put the values into a cell of a cell array
indices=[];
n=1;
for k = 1 : numRegions
if measurements(k).Area >= duration_threshold
% Area (length) is duration_threshold or greater, so store the values.
out{n} = measurements(k).PixelValues;
indices{n} = measurements(k).PixelIdxList;
n=n+1;
end
end
% Display the regions that meet the criteria
% celldisp(out);
for p = 1 : length(indices)
i_start = indices{p}(1);
i_end = indices{p}(end);
s_start = datetime(time(i_start), 'convertfrom', 'posixtime', 'Format', 'MM/dd/yy HH:mm:ss.SSS');
s_end = datetime(time(i_end), 'convertfrom', 'posixtime', 'Format', 'MM/dd/yy HH:mm:ss.SSS');
s_start.Format = 'HH:mm:ss.SSS';
s_end.Format = 'HH:mm:ss.SSS';
timestamps{p} = [s_start, s_end];
end
end