% Script to apply the physically based method from Pope et al. (2016) to % the red band of Sentinel-2 to derive lake depths for a series of % satellite imagery. Uses the same parameters as for Landsat-8 for % Sentinel-2 given the close overlap of the bands for the two satellites % This script should only be run following the validation of the % physically based method applied to Sentinel-2 against the same method % applied to Landsat-8. This is done within the script called % S2_L8_depth_physical_comp.m There are other scripts that do similar % things but use empirical (depth vs. reflectance) relationships instead of % the physical ones that are dealt with here. % Andrew Williamson, University of Cambridge, September 2017 %% Specify input directories and parameters main_dir = '/Users/Andrew/Documents/MATLAB/PhD_GrIS_SGLs_S2/'; S2_folder = '/Users/Andrew/Documents/MATLAB/PhD_GrIS_SGLs_S2/outputs_S2_2017'; S2_bands_folder = [S2_folder '/bands']; S2_NDWI_folder = [S2_folder '/NDWI_masks']; S2_mask_folder = [main_dir '/masks']; S2_out_depth_folder = [S2_folder '/depths_physical']; S2_out_info_folder = [S2_folder '/lake_info_physical']; S2_pix_siz = 0.01; S2_pix_siz_m = 10; % specifying S2 resolution g_red = 0.7507; % g value from Pope et al. (2016) for Landsat-8 S2_scalar = 10000; %% Make directories for saving files mkdir(S2_out_depth_folder); mkdir(S2_out_info_folder); %% Load data into workspace cd(S2_bands_folder); S2_red_list = dir('*_red.mat'); cd(S2_NDWI_folder); S2_NDWI_list = dir('*_NDWI_SGL_mask.mat'); cd(S2_mask_folder); GrIS_mask = imread('S2_mask_10m.tif'); num_img = length(S2_red_list); string = num2str(num_img); disp(['Calculating lake depths for ' string ' Sentinel-2 images.']) disp('Lake depths calculated with physically based methods of Pope et al. (2016).') for i = 1:num_img % get data and write image details to command line dates = {S2_red_list(:).name}; dates = dates'; % transpose matrix as it's nicer this way current_date = dates{i,1}; year = current_date(1,4:7); month = current_date(1,8:9); day = current_date(1,10:11); curr_date_str = [day '/' month '/' year]; % format date string nicely curr_imnum_str = num2str(i); disp(' ') disp(['Now analysing image from ' curr_date_str ', image number '... curr_imnum_str ' of ' string '.']) % new variable for short name string for use later short_dt_str = [year month day]; % get the filenames for loading reflectance data and masks red_fname = S2_red_list(i).name; NDWI_fname = S2_NDWI_list(i).name; % load the data into the workspace (note that the red_band % and NDWI_defined_SGLs names below are based on how the original % variables were saved to the .mat files - change if needed) cd(S2_bands_folder); disp('Loading red-band surface reflectance data.') load(red_fname); % loads variable with name red_band cd(S2_NDWI_folder); disp('Loading NDWI-derived lake area mask.') load(NDWI_fname); lake_mask = NDWI_defined_SGLs; % convert surface reflectance to double format for analysis and mask % no-data values (zeroes in S2), also mask the GrIS edge red_band = double(red_band); red_band = red_band ./ S2_scalar; red_band(red_band == 0) = NaN; % get darkest pixel value from the red band rinf_red = min(min(red_band)); rinf_red_str = num2str(rinf_red,'%.3f'); disp(['Darkest red band pixel value is ' rinf_red_str '.']); if rinf_red > 0.1 disp('Warning: darkest red band value is high - manually check.'); end red_band(GrIS_mask == 0) = NaN; lake_mask(GrIS_mask == 0) = NaN; % get a lake-number-labelled NDWI mask lake_mask(isnan(lake_mask)) = 0; [lake_mask_labelled,lakenum] = bwlabel(lake_mask,8); % default 8 connectivity % set up blank arrays for depth information lake_depths = zeros(size(lake_mask)); lakeinfo = zeros(lakenum,8); % set up structuring element e = strel('disk',1); % to dilate lake mask by 1 pixel for ad calcs % do the lake depth calculations for each lake in the array disp(['Calculating depths for all ' num2str(lakenum) ' lakes on the image.']); for n = 1:lakenum % defines number of times to loop % show the lake number being worked on at the command line disp(['Now working on lake number ' num2str(n) ' of a total ' num2str(lakenum) ' lakes.']) lmask = zeros(size(lake_mask)); % creates lake mask with size as above tmpref_red = red_band; % tmpref is equal to values of red band % Get the individual lake mask lmask(lake_mask_labelled == n) = 1; % put lake area information into the lakeinfo array lakeinfo(n,1) = n; % # of lake to 1st col of 'lakeinfo' array lakeinfo(n,2) = sum(sum(lmask)); % sum of pixels within lake lakeinfo(n,7) = (sum(sum(lmask)) * S2_pix_siz^2); % lake area in km^2 % Dilate lake mask by one ring of pixels to find bare-ice albedo lmaskdil = imdilate(lmask,e); % e is defined above - structuring element ledge = lmaskdil - lmask; % just get the edge pixels edgepix_red = (red_band(ledge == 1)); % find the red edge pixel refl ad_red = nanmean(edgepix_red); % take the mean of these values (ignore NaNs) lakeinfo(n,5) = ad_red; % put this in lake info file % mask reflectance arrays tmpref_red(lmask == 0) = NaN; % calculate lake depths with methods of Sneed and Hamilton (2007); % see details in Pope and others (2016) or Williamson and others (2017) tmpdepth_red = (log(ad_red - rinf_red) - log(tmpref_red - rinf_red)) / (-g_red); % add the two depth values together and take the average of them tmpdepth = tmpdepth_red; tmpdepth = -tmpdepth; % to get thee back to the correct +ve values % (this line may cause issues as it was added after the testing and % without re-running, so if it does, be sure to remove it and then just % take negatives of all lake depth arrays in the lake tracking / % drainage analysis) tmpdepth(lmask == 0) = 0; lakeinfo(n,3) = sum(sum(tmpdepth)); % sum of all pixel depths lakeinfo(n,4) = max(max(tmpdepth)); % max lake depth lakeinfo(n,6) = mean(mean(tmpdepth)); % mean lake depth lakeinfo(n,8) = (sum(sum(tmpdepth)) * (S2_pix_siz_m)^2); % volume (m^3) % add calculated depths to lake depth array lake_depths = lake_depths + tmpdepth; end % end of lake-depth calculation loop % mask the lake depth arrays off the ice-sheet edge with NaNs lake_depths(GrIS_mask == 0) = NaN; % get new variable names from above new_depth_var_name = ['S2_' short_dt_str '_depths']; new_info_var_name = ['S2_' short_dt_str '_lake_info']; % save the output files cd(S2_out_depth_folder); save([new_depth_var_name '.mat'],'lake_depths') cd(S2_out_info_folder); save([new_info_var_name '.mat'],'lakeinfo') % display progress update to command line disp(' '); progress = 100/num_img * i; % percentage based on number of images disp(['Lake depth calculations now ' num2str(progress,'%.1f') '% complete.']) end cd(main_dir)