%% Unpack Data

% A script to segment panchromatic imagery and pool maps stored in the
% 'FinalProducts' directory into buffered one kilometer square tiles in the
% 'Workspace' directory. This script needs to be run first if you plan to
% experiment with training a UNet or performing histogram matching.

% author:   cabolt
% date:     11.12.2021, La Mesilla NM

clear; clc; close all

%% Loop through the survey areas

sites   = dir(fullfile('FinalProducts', 'imagery'));
sites   = {sites(~ismember({sites.name}, {'.', '..'})).name};

for k = 1:numel(sites)
    % generate a list of dates from which data is available at site 
    dates   = dir(fullfile('FinalProducts', 'imagery', sites{k}, '*.tif'));
    dates   = cellfun(@(c) c(1:end-4), {dates.name}, 'uniformoutput', false);
    
    % find the site bounding box in utm coordinates
    info    = geotiffinfo(fullfile('FinalProducts', 'imagery', ...
              sites{k}, [dates{1} '.tif']));
    bbox    = info.BoundingBox;
    e0      = (bbox(1,1)/1000) : (bbox(2,1)/1000 - 1);
    e0      = arrayfun(@(a) num2str(a), e0, 'uniformoutput', false);
    n0      = (bbox(1,2)/1000) :  (bbox(2,2)/1000 - 1);
    n0      = arrayfun(@(a) num2str(a), fliplr(n0), 'uniformoutput', false);
    
    for kk = 1:numel(dates)
        % segment images into buffered 1sqkm tiles, one date at a time
        wvimg   = imread(fullfile('FinalProducts', 'imagery', sites{k}, ...
                  [dates{kk} '.tif']));
        wvimg   = padarray(wvimg, [28 28], 'symmetric');
        pools   = imread(fullfile('FinalProducts', 'pools_raster', sites{k}, ...
                  [dates{kk} '.tif']));
        pools   = padarray(pools, [28 28], 'symmetric');
        
        for i = 1:numel(n0)
            for j = 1:numel(e0)
                sqkm    = [e0{j} '_' n0{i}];
                if ~isfolder(fullfile('Workspace', 'imagery', sites{k}, sqkm))
                    mkdir(fullfile('Workspace', 'imagery', sites{k}, sqkm));
                end
                if ~isfolder(fullfile('Workspace', 'results', 'post_CRF', sites{k}, sqkm))
                    mkdir(fullfile('Workspace', 'results', 'post_CRF', sites{k}, sqkm));
                end
                
                wvimgTile   = wvimg((i-1)*2000+1:i*2000+56, ...         % segmented satellite imagery has 14 meter buffer
                                    (j-1)*2000+1:j*2000+56);
                poolsTile   = pools((i-1)*2000+29:i*2000+28, ...        % segmented pools map has no buffer
                                    (j-1)*2000+29:j*2000+28);
                imwrite(wvimgTile, fullfile('Workspace', 'imagery', ...
                    sites{k}, sqkm, [dates{kk} '.tif']));
                imwrite(poolsTile, fullfile('Workspace', 'results', 'post_CRF', ...
                    sites{k}, sqkm, [dates{kk} '.tif']));
            end
        end
    end
    
    disp(['done unpacking imagery from ' num2str(k) ' of ' num2str(numel(sites)) ...
        ' sites! (' sites{k} ')']); 
end
