function out=hedgeMakeStationMap(stationBatch,boundaries,cMap,displayBool)
% hedgeMakeStationMap reads in model data calibrated to regional weather
% stations and makes a map of the contiguous 48 states using that data.
% Used to make Figure 6C.
%=========================================================================
% Usage:
% hedgeMakeStationMap(stationBatch,boundaries,cMap,displayBool)
% Inputs:
%       * stationBatch = cell array of calibrated station data (output of
%            hedgeGeography.m
%       * boundaries = state boundary image, var boundaries in
%            hedgeStateBoundaries.mat
%       * cMap = colorMap for the background interpolation and data points,
%            posneg2 in hedgeColorMaps.mat
%       * displayBool = 1x4 array indicating which display options to use
%            in the final figure, as follows:
%                 [X - - -] = display the interpolated background map
%                 [- X - -] = display state boundaries
%                 [- - X -] = display black points at station locations on the pixel map
%                 [- - - X] = superimpose vector-based data points on the map
% Output:
%       * out = matrix of the background image BvH values

% set the latitude and longitude resolution, and the bounds of each
latRes=.1;
latRange=(28:latRes:53)*pi()/180;
longRes=.1;
longRange=(-130:longRes:-65)*pi()/180;

% extract data on stations from station data cell array
latLongs=cell2mat(stationBatch(:,4:5))*pi()/180;
BvHs=cell2mat(stationBatch(:,3));

%initialize a variable to store the distance of each pixel to all stations
distsToStats=zeros(size(BvHs));

%initialize the map
llMap=zeros(length(latRange),length(longRange));
llMap2=llMap;

latTick=1; % loop parameter

% if showing the background interpolated map
if displayBool(1)==1
    %for all lat and long pixels
    for i=latRange;
        longTick=1;
        for j=longRange;
            
            %for each pixel, calculate the great circle distance to each station
            for k=1:size(latLongs,1)
                a=sin((i-latLongs(k,1))/2)^2 + cos(i)*cos(latLongs(k,1))*sin((j-latLongs(k,2))/2)^2;
                c=2*atan2(sqrt(a),sqrt(1-a));
                distsToStats(k)=c;
            end
            
            %scale and sum the distances by an exponential term, to weight
            %BvH values of close stations more highly
            statWeights=1./(distsToStats).^3;       %change this exponent to change interpolation
            pxColor=sum(BvHs.*statWeights)/sum(statWeights);
            
            %save this pixel value
            llMap(latTick,longTick)=pxColor;
            
            %increment the ticker
            longTick=longTick+1;
        end
        latTick=latTick+1;
    end
    llMap2=llMap; %save into working map
end

% if showing the state boundaries
if displayBool(2)==1
    if size(boundaries,1)>1 % you can pass a scalar instead of the boundary display bool if you don't want to show it
        llMap2(flipud(boundaries))=Inf;
    end
end

% if showing pixels at station locations
if displayBool(3)==1
    
    %for all stations
    for i=1:size(latLongs,1)
        % get the station lat long
        lati=latLongs(i,1);
        longi=latLongs(i,2);
        
        % find the closest pixel
        sRow=(latRange-lati).^2==min((latRange-lati).^2);
        sCol=(longRange-longi).^2==min((longRange-longi).^2);
        
        % set those pixels to black/extreme color
        llMap2(sRow,sCol)=-Inf;
    end
end

%set the scale of the colormap
colorScale=8;
colorOffset=129;

%start the figure
figure;
hold on;

% if displaying something on the pixel map background
if sum(displayBool(1:3))>0;
    image((llMap2)*colorScale+colorOffset); %display the image with the right color
    colormap(cMap)
end

% if displaying vector-based versions of the station positions on top of
% the pixel map
if displayBool(4)==1
    %start a vector for the colors of the points
    BvHs=round(BvHs*colorScale+colorOffset);
    BvHs(BvHs>256)=256;
    BvHs(BvHs<1)=1;
    %get the pixel positions of these lat longs
    latLongs=latLongs/(pi()/180);
    latLongs(:,1)=(latLongs(:,1)-28)*10+1;  % these values would need to be adjusted if the bounds of the map are changed
    latLongs(:,2)=(latLongs(:,2)-(-130))*10+1;
    %draw the points
    scatter(latLongs(:,2),latLongs(:,1),40,cMap(BvHs,:),'filled','LineWidth', 1, 'MarkerEdgeColor', 'k')
end;

%frame the image
xlim([1 651]);
ylim([1 251]);

%save the BvH map data as output
out=flipud(llMap);