% created by Nils Hase 2016 (nilshase@math.uni-bremen.de)

% This script will divide the US into several regions just based on
% geographical ideas

function [regions,region_name] = get_regions(plotter)

load('synthetic_methane_data.mat','land_grids')

% Split up using the get_index function
regions = {};
act_set = [get_index(-130,-104,32,40);get_index(-130,-112,41,41);get_index(-103,-103,37,40)];
regions{1} = act_set;
act_set = [get_index(-135,-103,49,55)];
regions{2} = act_set;
act_set = setdiff(get_index(-130,-105,40,50),[regions{1};regions{2}]);
regions{3} = act_set;
act_set = setdiff([get_index(-91,-75,25,35)],[get_index(-91,-91,35,35);get_index(-91,-91,29,30)]);
regions{4} = act_set;
act_set = [get_index(-90,-68,36,42);get_index(-86,-68,43,44);get_index(-91,-91,38,41);get_index(-81,-68,45,45);get_index(-77,-68,46,46);get_index(-74,-68,47,47)];
regions{5} = act_set;
act_set = [get_index(-120,-107,25,31);get_index(-106,-106,25,30);get_index(-105,-105,25,29);get_index(-104,-104,25,28);get_index(-103,-102,25,29);get_index(-101,-101,25,27);get_index(-100,-100,25,26);get_index(-99,-98,25,25)];
regions{6} = act_set;
act_set = setdiff(get_index(-110,-85,25,39),vertcat(regions{1:6}));
regions{7} = act_set;
act_set = setdiff(get_index(-105,-87,40,48),vertcat(regions{1:7}));
regions{8} = act_set;
act_set = setdiff(get_index(-80,-50,40,55),vertcat(regions{1:8}));
regions{9} = act_set;
act_set = setdiff((1:1469)',vertcat(regions{1:9}));
regions{10} = act_set;

region_name = {'SW','Can W','NW','SE','NE','Mex','S','N','Can E','Can'};

resort = [4,5,7,8,1,3,6,9,10,2];

% Plot the regions
color = [1.7,2.3,3.2,3.9,5.3,6.5,7.2,8,9,10];
xx = [-115,-120,-117,-87,-84,-110,-98,-98,-75,-95];
yy = [37,52,45,33,39.5,28,34,44.5,52,52];

if plotter == 1
    f = zeros(length(land_grids),1);

    for k = 1:length(regions);
        f(regions{k}) = color(k);
    end

    map_result_less(f);
    caxis([0 10])
    colormap jet
    for k = 1:length(regions)
        text(xx(k),yy(k),region_name{k},'fontsize',12,'fontweight','bold');
    end
    colorbar off
end

region_name = region_name(resort);
regions = regions(resort);

end
