close all;
clear all;
clc;

load conus;
S = shaperead('landareas','UseGeoCoords',true);
radius = 6370;

% set plotting range
xmax = -68;
xmin = -125.3;
ymax = 49.5;
ymin = 24.5;

% xmaxb = xmax+1;
% xminb = xmin-1;
% ymaxb = ymax+1;
% yminb = ymin-1;


[clat,clon,rlat,rlon,subNo]=textread('box.txt','%f %f %f %f %d');
numbox=numel(clat);

filename = strcat('scaler_',num2str(clat(1)),'_',num2str(clon(1)),'.dat');
[lat, lon, elev, scaler] = textread(char(filename), '%f %f %f %f'); 
for k=2:numbox
    
    % read data
    filename = strcat('scaler_',num2str(clat(k)),'_',num2str(clon(k)),'.dat');
    [lat1, lon1, elev1, scaler1] = textread(char(filename), '%f %f %f %f');  
    lat=[lat;lat1];
    lon=[lon;lon1];
    elev=[elev;elev1];
    scaler=[scaler;scaler1];  
    clear lat1 lon1 elev1 scaler1
end


outfigdir = './';          % figure output directory, needs / on end

xi = xmin:0.05:xmax;
yi = ymin:0.05:ymax;
[XI,YI] = meshgrid(xi,yi);

% xib = xminb:0.04:xmaxb;
% yib = yminb:0.04:ymaxb;
% [XIB,YIB] = meshgrid(xib,yib);

% ZI = griddata(lon,lat,elev,XIB,YIB,'cubic');
ZZV = griddata(lon,lat,scaler,XI,YI,'cubic');
elev_scl=[50,500,1500,2500];

% set scaler for analysis 
cmap=zeros(25,3);
for i=1:26
    cmap(i,:)=(26-i)*0.04;
end

v1=0;
v2=1;
h=0.1;
t=(v2-v1)/24;
v=v1:t:v2;
scl = v1:t:v2;

figure(3)
box off;
axis off;
% axesm('MapProjection','lambertstd','Grid','on','MapLatLimit',[ymin ymax],'MapLonLimit',[xmin xmax],'MeridianLabel','on','ParallelLabel','on','MlineLocation',[xmin (xmin+xmax)/2 xmax] ,'MLabelRound',-1,'PlineLocation',[ymin (ymin+ymax)/2 ymax],'PLabelRound',-1);
axesm('MapProjection','lambertstd','Grid','on','MapLatLimit',[ymin ymax],'MapLonLimit',[xmin xmax]);
contourfm(YI,XI,ZZV,v);
shading flat;
colormap(cmap);
cbh=colorbar;
caxis([v1 v2]);
set(cbh,'YTick',v1:h:v2);
% [C,h]=contourm(YIB,XIB,ZI,elev_scl,'k--');
% clabel(C,h,elev_scl,'FontName','Times-Roman','FontSize',6,'Color','k','Rotation',0);
numbox=numel(subNo);

% for k=1:numbox
%     ymin_box(k)=clat(k)-rlat(k);
%     ymax_box(k)=clat(k)+rlat(k);
%     xmin_box(k)=clon(k)-rlon(k);
%     xmax_box(k)=clon(k)+rlon(k);
% %     txt=strcat(num2str(sstvt(k)),'/',num2str(rawsNo(k)));
%     if (subNo(k)==20)
%         plotm([ymin_box(k) ymin_box(k)],[xmin_box(k) xmax_box(k)],'k','LineWidth',1);
%         plotm([ymax_box(k) ymax_box(k)],[xmax_box(k) xmin_box(k)],'k','LineWidth',1);
% %         plotm([ymin_box(k) ymax_box(k)],[xmax_box(k) xmax_box(k)],'k','LineWidth',1);
%         plotm([ymax_box(k) ymin_box(k)],[xmin_box(k) xmin_box(k)],'k','LineWidth',1);
% %         textm((clat(k)+0.6),(clon(k)-0.6),txt,'FontName','Times-Roman','FontSize',8,'FontWeight','bold','Color','k','HorizontalAlignment','center','Clipping','on');
%     else
%         plotm([ymin_box(k) ymin_box(k)],[xmin_box(k) xmax_box(k)],'k','LineWidth',1);
%         plotm([ymax_box(k) ymax_box(k)],[xmax_box(k) xmin_box(k)],'k','LineWidth',1);
%         plotm([ymin_box(k) ymax_box(k)],[xmax_box(k) xmax_box(k)],'k','LineWidth',1);
%         plotm([ymax_box(k) ymin_box(k)],[xmin_box(k) xmin_box(k)],'k','LineWidth',1);
% %         textm(clat(k),clon(k),txt,'FontName','Times-Roman','FontSize',8,'FontWeight','bold','Color','k','HorizontalAlignment','center','Clipping','on');
%     end
% 
% end
nums=numel(lon);

for i=1:nums
    if (lon(i)<-110)
        geoshow(statelat,statelon,'DisplayType','line','Color',[0.3,0.3,0.3],'LineStyle','-','LineWidth',0.1);
        geoshow([S.Lat], [S.Lon],'DisplayType','line','color',[0.3,0.3,0.3],'LineStyle','-','LineWidth',0.1);
    else
        geoshow(statelat,statelon,'DisplayType','line','Color',[0.3,0.3,0.3],'LineStyle','-','LineWidth',0.1);
        geoshow(uslat,uslon,'DisplayType','line','Color',[0.3,0.3,0.3],'LineStyle','-','LineWidth',0.1);
        geoshow([S.Lat], [S.Lon],'DisplayType','line','color',[0.3,0.3,0.3],'LineStyle','-','LineWidth',0.1);
    end
end


% geoshow(statelat,statelon,'DisplayType','line','Color','r','LineStyle','-','LineWidth',0.1);
% geoshow(uslat,uslon,'DisplayType','line','Color','r','LineStyle','-','LineWidth',0.1);
% geoshow([S.Lat], [S.Lon],'DisplayType','line','color','r','LineStyle','-','LineWidth',0.1);
filename = strcat(outfigdir,'shaded_map_gray.eps');
saveas(gcf,filename,'epsc2');