% (c) Rebecca Dell, University of Cambridge, June 2020. % Adapted from the work of Williamson et al. (2018) (doi:10.5194/tc-2018-56). % Who calculated depths following Pope et al. (2016) (doi:10.5194/tc-10-15-2016). % Function originally also used by Banwell et al. (2014) (doi:10.3189/2014AoG66A049), Arnold et al. % (2014) (doi:10.5194/tc-8-1149-2014), and Williamson et al. (2017) (doi:10.1016/j.rse.2017.04.032). function [lake_depths_c] = L8_depth_mask(rinf_red,rinf_pan,B8,B4,combined_matrix,filled_combined) e1 = strel('disk',1); % to dilate lake mask by 1 pixel e2 = strel('disk',2); % to dilate lake mask by 2 pixels g_red = 0.7507; % value from Pope and others (2016) g_pan = 0.3817; % value from Pope and others (2016) [combined_mask_labelled,lakenum_c] = bwlabel(combined_matrix,8); lake_depths_c = zeros(size(rinf_red)); for c = 1:lakenum_c % defines number of times to loop lmask_c = zeros(size(B4)); % creates lake mask with size as above tmpref_pan_c = B8; % tmpref is equal to values of pan band tmpref_red_c = B4; % tmpref is equal to values of red band % Get the individual lake mask lmask_c(combined_mask_labelled == c) = 1; % Dilate lake mask by one ring of pixels to find bare-ice albedo lmaskdil_c = imdilate(lmask_c,e1); % e is defined above - structuring element lmaskdil_c2 = imdilate(lmask_c,e2); % e is defined above - structuring element ledge_c = lmaskdil_c2 - lmaskdil_c; % just get the edge pixels edgepix_pan_c = (tmpref_pan_c(ledge_c == 1)); % find the pan edge pixel refl edgepix_red_c = (tmpref_red_c(ledge_c == 1)); % find the red edge pixel refl ad_pan_c = nansum(nansum(edgepix_pan_c))/sum(sum(~isnan(edgepix_pan_c))); % take the mean of these values (ignore NaNs) ad_red_c = nansum(nansum(edgepix_red_c))/sum(sum(~isnan(edgepix_red_c))); % take the mean of these values (ignore NaNs) % mask reflectance arrays tmpref_pan_c(lmask_c == 0) = NaN; tmpref_red_c(lmask_c == 0) = NaN; % Calculate Depths tmpdepth_pan = (log(ad_pan_c - rinf_pan) - log(tmpref_pan_c - rinf_pan)) / (g_pan); tmpdepth_red = (log(ad_red_c - rinf_red) - log(tmpref_red_c - rinf_red)) / (g_red); tmpdepth_red(tmpdepth_red<0) = NaN; tmpdepth_pan(tmpdepth_pan<0) = NaN; % add the two depth values together and take the average of them tmpdepth_c = tmpdepth_red + tmpdepth_pan; tmpdepth_c = tmpdepth_c./2; % divide by two (i.e. take average) tmpdepth_c(lmask_c == 0) = 0; % if lmask value = 0, then tmpdepth = 0 mean_lake_depth_c = nansum(nansum(tmpdepth_c))/sum(sum(~isnan(tmpdepth_c))); % add calculated depths to lake depth array lake_depths_c = lake_depths_c + tmpdepth_c; %% Calculate depth of filled holes % filled_combined 1's = filled holes holes_mask_c = filled_combined; template_filled_mask_c = zeros(size(filled_combined)); % creates lake mask with size as above specific_combined = zeros(size(filled_combined)); % find specific lake specific_combined(combined_mask_labelled == c) = 2; % Add the holes mask and specific lake mask holes_and_combined = holes_mask_c + specific_combined; % Find where there are holes in the specific lake template_filled_mask_c(holes_and_combined ==3)=1; % assign the specific lake holes the mean lake value template_filled_mask_c(template_filled_mask_c==1)= mean_lake_depth_c; % add the new hole depths to the overall depth mask lake_depths_c = lake_depths_c + template_filled_mask_c; end end