%% Script Description 

clear all 
close all

set(0,'defaultfigureposition',[500 500 500 420])
format shortg

%load('mod_nlcd2011_lnd.mat')

%% Determin the values and ncread the files into matlab variables.

% Variables that won't be modified.
lat = ncread('../output/nlcd2011-dec31/geo_em.d03.nc','XLAT_M');
lon = ncread('../output/nlcd2011-dec31/geo_em.d03.nc','XLONG_M');
hgt = ncread('../output/nlcd2011-dec31/geo_em.d03.nc','HGT_M');
flg = ncread('../output/nlcd2011-dec31/geo_em.d03.nc','LANDMASK');

% Variables to modify.
fin_lnd =  ncread('../output/nlcd2011-dec31/geo_em.d03.nc','LU_INDEX');
lnd2011 = ncread('../output/nlcd2011-dec31/geo_em_new.d03.nc','LU_INDEX');
lnd2006 = ncread('../meso3/wrf-nlcd40/WPS/geo_em.d03.nc','LU_INDEX');
lndUSGS = ncread('../output/usgs-dec31/geo_em.d03.nc','LU_INDEX');

fin_flg = ncread('../output/nlcd2011-dec31/geo_em.d03.nc','LANDMASK');
flg2011 = ncread('../output/nlcd2011-dec31/geo_em_new.d03.nc','LANDMASK');
flg2006 = ncread('../meso3/wrf-nlcd40/WPS/geo_em.d03.nc','LANDMASK');
flgUSGS = ncread('../output/usgs-dec31/geo_em.d03.nc','LANDMASK');

fin_lnd = double(fin_lnd);
lnd2011 = double(lnd2011);
lnd2006 = double(lnd2006);
lndUSGS = double(lndUSGS);

lat = double(lat(:,:,1));
lon = double(lon(:,:,1));

lats = lat';
lons = lon';

hgt = double(hgt);
flg = double(flg);

latlim = [min(lat(:)) max(lat(:))];
lonlim = [min(lon(:)) max(lon(:))];

%%

id  = {'ISFS1',   'ISFS2',   'ISFS3',   'ISFS4',   'ISFS5',   'ISFS6',   'ISFS7',   'ISS-N',  'ISS-S'};
lon_st = [-112.1031, -111.9136, -111.8448, -111.9946, -111.8045, -112.0606, -111.9729, -111.925, -111.9252];
lat_st = [40.7785,   40.7374,   40.6980,    40.6810,   40.6127,   40.5847,  40.5140,   40.61,    40.6006];

% Find closest grid point to station locations

for y = 1:length(lons(:,1))-1
   for x = 1:length(lons(1,:))-1
      for i = 1:length(lon_st)
         if (lons(y,x) < lon_st(i) & lon_st(i) < lons(y,x+1)) & (lats(y,x) < lat_st(i) & lat_st(i) < lats(y+1,x))
            disp('yay')
            x_grid(i) = x;
            y_grid(i) = y;
         end
      end
   end
end

lon_grid = [lons(y_grid(1),(x_grid(1))), lons(y_grid(2),(x_grid(2))), lons(y_grid(3),(x_grid(3))), lons(y_grid(4),(x_grid(4))), lons(y_grid(5),(x_grid(5))), lons(y_grid(6),(x_grid(6))), lons(y_grid(7),(x_grid(7)))];
lat_grid = [lats(y_grid(1),(x_grid(1))), lats(y_grid(2),(x_grid(2))), lats(y_grid(3),(x_grid(3))), lats(y_grid(4),(x_grid(4))), lats(y_grid(5),(x_grid(5))), lats(y_grid(6),(x_grid(6))), lats(y_grid(7),(x_grid(7)))];


%% Set boundaries for land use plotting.

x_min = 150;    % Boundaries for modification to snow.
x_max = 425;
y_min = 180;
y_max = 425;

lx_min = 160; 
lx_max = 240;
ly_min = 290;
ly_max = 405;

lake_lnd2011 = lnd2011(lx_min:lx_max,ly_min:ly_max);
lake_lnd2006 = lnd2006(lx_min-5:lx_max-5,ly_min-5:ly_max-5);
lake_lndUSGS = lndUSGS(lx_min:lx_max,ly_min:ly_max);
lake_mod = fin_lnd(lx_min:lx_max,ly_min:ly_max);

vx_min = 212; 
vx_max = 240;
vy_min = 290;
vy_max = 322;

valley_lnd2011 = lnd2011(vx_min:vx_max,vy_min:vy_max);
valley_lnd2006 = lnd2006(vx_min-5:vx_max-5,vy_min-5:vy_max-5);
valley_lndUSGS = lndUSGS(vx_min:vx_max,vy_min:vy_max);
valley_mod = fin_lnd(vx_min:vx_max,vy_min:vy_max);

%% Find all unique instances of land use type in domain for landuse frac.

lake_lnd2011_unique = unique(lake_lnd2011);
lake_lnd2006_unique = unique(lake_lnd2006);
lake_lndUSGS_unique = unique(lake_lndUSGS);
lake_mod_unique = unique(lake_mod);

valley_lnd2011_unique = unique(valley_lnd2011);
valley_lnd2006_unique = unique(valley_lnd2006);
valley_lndUSGS_unique = unique(valley_lndUSGS);
valley_mod_unique = unique(valley_mod);

nlcd_lnd_id = {'Evergreen Needleleaf Forest', ...  
          'Evergreen Broadleaf Forest', ...
          'Deciduous Needleleaf Forest', ... 
          'Deciduous Broadleaf Forest', ... 
          'Mixed Forest', ... 
          'Closed Shrubland', ... 
          'Open Shrubland', ... 
          'Woody Savanna', ... 
          'Savanna', ... 
          'Grassland', ... 
          'Permanent Wetland', ... 
          'Cropland', ... 
          'Urban and Built-Up', ...  
          'Cropland / Natural Veg. Mosaic', ... 
          'Permanent Snow', ... 
          'Barren / Sparsely Vegetated', ... 
          'IGBP Water', ... 
          'Unclassified', ...
          'Fill Value', ...
          'Unclassified', ...
          'Open Water', ...
          'Perennial Ice/Snow', ... 
          'Developed Open Space', ...  
          'Developed Low Intensity', ...  
          'Developed Medium Intensity', ...  
          'Developed High Intensity', ...  
          'Barren Land', ... 
          'Deciduous Forest', ... 
          'Evergreen Forest', ...
          'Mixed Forest', ... 
          'Dwarf Scrub', ... 
          'Shrub/Scrub', ... 
          'Grassland/Herbaceous', ... 
          'Sedge/Herbaceous', ... 
          'Lichens', ... 
          'Moss', ... 
          'Pasture/Hay', ... 
          'Cultivated Crops', ... 
          'Woody Wetland', ... 
          'Emergent Herbaceous Wetland' 
          };
      
usgs_lnd_id = {'Urban and Built-Up Land', ...   
     'Dryland Cropland and Pasture', ...  
     'Irrigated Cropland and Pasture', ...  
     'Mixed Dryland/Irrigated Cropland and Pasture', ...  
     'Cropland/Grassland Mosaic', ... 
     'Cropland/Woodland Mosaic', ...  
     'Grassland', ...  
     'Shrubland', ...  
     'Mixed Shrubland/Grassland', ... 
     'Savanna', ...  
     'Deciduous Broadleaf Forest', ...  
     'Deciduous Needleleaf Forest', ...  
     'Evergreen Broadleaf Forest', ... 
     'Evergreen Needleleaf Forest', ...   
     'Mixed Forest', ...  
     'Water Bodies', ...  
     'Herbaceous Wetland', ...  
     'Wooded Wetland', ...  
     'Barren or Sparsely Vegetated', ...  
     'Herbaceous Tundra', ...  
     'Wooded Tundra', ...  
     'Mixed Tundra', ...  
     'Bare Ground Tundra', ...  
     'Snow or Ice', ...  
     'Playa', ...  
     'Lava', ...    
     'White Sand'
     };
 
%% New Land Use Plots
 
usgs_lnd_id_10 = {'Water', ...   
    'Wetland', ...
    'Devloped Urban', ...    
    'Evergreen Forest', ... 
    'Deciduous Forest', ...    
    'Irrigated Cropland and Pasture', ...
    'Cropland/Grassland Mosaic', ... 
    'Grassland', ...  
    'Shrubland', ...
    'Barren Land', ...
    'High Albedo Surface', ...  
    };

cmap_usgs = colormap([0.0 0.0 1.0; ...
    0.0 1.0 1.0; ...
    0.4 0.4 0.4; ...
    0.0 0.4 0.2; ...
    0.2 0.6 0.2; ...
    0.3 0.7 0.0; ...
    0.8 1.0 0.2; ...
    0.0 1.0 0.0; ...
    0.8 0.4 0.2; ...
    0.6 0.4 0.0; ...   
    1.0 1.0 1.0]);

nlcd_lnd_id_10 = {'Water', ... 
    'Wetland', ...
    'Developed High Intensity', ...
    'Developed Medium Intensity', ... 
    'Developed Low Intensity', ...
    'Developed Open Space', ...  
    'Evergreen Forest', ...
    'Deciduous Forest', ... 
    'Cultivated Crops', ...
    'Pasture/Hay', ...
    'Grassland', ...
    'Shrubland', ... 
    'Barren Land' 
    };

cmap_2011 = colormap([0.0 0.0 1.0; ...
    0.0 1.0 1.0; ...
    0.3 0.3 0.3; ...
    0.4 0.4 0.4; ...
    0.5 0.5 0.5; ...
    0.6 0.6 0.6; ...
    0.0 0.4 0.2; ...
    0.2 0.6 0.2; ...
    0.3 0.7 0.0; ...
    0.8 1.0 0.2; ...
    0.0 1.0 0.0; ...
    0.8 0.4 0.2; ...
    0.6 0.4 0.0
    ]);
 
%lndUSGS_10 = zeros(size(lndUSGS));
lndUSGS_id_10 = [16,21,1,14,11,3,5,7,8,19,24];
maxalb_usgs =   [0.7,0.55,0.46,0.52,0.58,0.66,0.68,0.7,0.6,0.75,0.82]

for x = 1:length(lndUSGS(:,1))
    for y = 1:length(lndUSGS)
        for i = 0:length(usgs_lnd_id_10)-1
            if lndUSGS(x,y) == lndUSGS_id_10(i+1)
                lndUSGS_10(x,y) = i;
                albUSGS_10(x,y) = maxalb_usgs(i+1);
            end
        end
    end
end     

lnd2011_id_10 = [17,40,26,25,24,23,29,28,38,37,33,32,27];
maxalb_nlcd =   [0.7,0.68,0.7,0.75,0.75,0.75,0.42,0.56,0.75,0.75,0.7,0.8,0.8]

for x = 1:length(lnd2011(:,1))
    for y = 1:length(lnd2011)
        for i = 0:length(nlcd_lnd_id_10)-1
            if lnd2011(x,y) == lnd2011_id_10(i+1)
               lnd2011_10(x,y) = i;
               albNLCD_10(x,y) = maxalb_nlcd(i+1);
            end
            if fin_lnd(x,y) == lnd2011_id_10(i+1)
               fin_lnd_10(x,y) = i;
            end   
            if lnd2006(x,y) == lnd2011_id_10(i+1)
               lnd2006_10(x,y) = i;
            end  
        end
    end
end    

station_land_use = [fin_lnd_10((x_grid(1)-1):(x_grid(1)+1),(y_grid(1)-1):(y_grid(1)+1)); ...
    fin_lnd_10((x_grid(2)-1):(x_grid(2)+1),(y_grid(2)-1):(y_grid(2)+1)); ...
    fin_lnd_10((x_grid(3)-1):(x_grid(3)+1),(y_grid(3)-1):(y_grid(3)+1)); ...
    fin_lnd_10((x_grid(4)-1):(x_grid(4)+1),(y_grid(4)-1):(y_grid(4)+1)); ...
    fin_lnd_10((x_grid(5)-1):(x_grid(5)+1),(y_grid(5)-1):(y_grid(5)+1)); ...
    fin_lnd_10((x_grid(6)-1):(x_grid(6)+1),(y_grid(6)-1):(y_grid(6)+1)); ...
    fin_lnd_10((x_grid(7)-1):(x_grid(7)+1),(y_grid(7)-1):(y_grid(7)+1))];

stn_unique = unique(station_land_use);

for i = 1:length(stn_unique)
    station_land_use(find(station_land_use == stn_unique(i))) = i-1 
end

station_land_use = reshape(station_land_use',[3,3,7]);


%%

close all

figure('Position', [100, 100, 650, 800])

subplot(2,1,1)
title('NLCD 2011 Land Use')
ax1 = usamap([lat(vx_min,vy_min),lat(vx_max,vy_max)],[lon(vx_min,vy_max),lon(vx_max,vy_min)]);
surfm(lat,lon,lnd2011_10)
colormap(ax1,cmap_2011)
h1 = colorbar;
caxis([0 13])
set(h1,'YTick',[0.5:1:length(nlcd_lnd_id_10)+1],'YTickLabel',nlcd_lnd_id_10)
mlabel('off')
hold on
mod = plotm(lat_grid, lon_grid, 'w.', 'MarkerSize', 20)
textm(lat_grid,lon_grid+0.01,{'ISFS1','ISFS2','ISFS3','ISFS4','ISFS5','ISFS6','ISFS7'},'FontSize',12)
hold on
iss = plotm(lat_st(8:9), lon_st(8:9), 'k.', 'MarkerSize', 20)
textm(lat_st(8)+0.01,lon_st(8)+0.01,'ISS-N','FontSize',12)
textm(lat_st(9)-0.01,lon_st(9)+0.01,'ISS-S','FontSize',12)
hold on
kslc = plotm(40.77,-111.97,'k.','MarkerSize', 20)
textm(40.77, -111.97+0.01,'KSLC','FontSize',12)
hold on
time_height = plotm(lat(225,302),lon(225,302),'k*')
hold on
v_sub = plotm(lat(vx_min:vx_max,vy_min),lon(vx_min:vx_max,vy_min),'k-','LineWidth',2)
hold on
v_sub = plotm(lat(vx_min:vx_max,vy_max),lon(vx_min:vx_max,vy_max),'k-','LineWidth',2)
hold on
v_sub = plotm(lat(vx_min,vy_min:vy_max),lon(vx_min,vy_min:vy_max),'k-','LineWidth',2)
hold on
v_sub = plotm(lat(vx_max,vy_min:vy_max),lon(vx_max,vy_min:vy_max),'k-','LineWidth',2)

setm(gca,'mlinelocation',0.2)

fig_nudge(ax1,-0.1,0)

fl = figm_label('a','ul',14,0.04,0.05,4500,4500);
set(fl,'color','k');

subplot(2,1,2)
title('USGS Land Use')
ax2 = usamap([lat(vx_min,vy_min),lat(vx_max,vy_max)],[lon(vx_min,vy_max),lon(vx_max,vy_min)]);
surfm(lat,lon,lndUSGS_10)
colormap(ax2,cmap_usgs)
h1 = colorbar;
caxis([0 11])
set(h1,'YTick',[0.5:1:length(usgs_lnd_id_10)+1],'YTickLabel',usgs_lnd_id_10)
hold on
mod = plotm(lat_grid,lon_grid,'w.','MarkerSize', 20)
textm(lat_grid,lon_grid+0.01,{'ISFS1','ISFS2','ISFS3','ISFS4','ISFS5','ISFS6','ISFS7'},'FontSize',12)
hold on
iss = plotm(lat_st(8:9),lon_st(8:9),'k.','MarkerSize', 20)
textm(lat_st(8)+0.01,lon_st(8)+0.01,'ISS-N','FontSize',12)
textm(lat_st(9)-0.01,lon_st(9)+0.01,'ISS-S','FontSize',12)
hold on
kslc = plotm(40.77,-111.97,'k.','MarkerSize', 20)
textm(40.77, -111.97+0.01,'KSLC','FontSize',12)
hold on
time_height = plotm(lat(225,302),lon(225,302),'k*')
hold on
v_sub = plotm(lat(vx_min:vx_max,vy_min),lon(vx_min:vx_max,vy_min),'k-','LineWidth',2)
hold on
v_sub = plotm(lat(vx_min:vx_max,vy_max),lon(vx_min:vx_max,vy_max),'k-','LineWidth',2)
hold on
v_sub = plotm(lat(vx_min,vy_min:vy_max),lon(vx_min,vy_min:vy_max),'k-','LineWidth',2)
hold on
v_sub = plotm(lat(vx_max,vy_min:vy_max),lon(vx_max,vy_min:vy_max),'k-','LineWidth',2)

setm(gca,'mlinelocation',0.2)
setm(gca,'mlabellocation',0.2)

fl = figm_label('b','ul',14,0.04,0.06,4500,4500);
set(fl,'color','k');

fig_nudge(ax2,-0.1,0.08)

%%

figure(1)

for i = 1:7
   subplot(2,4,i)
   ax(i) = imagesc(station_land_use(:,:,i));
   hold on
   plot(2,2,'k.','MarkerSize',50)
   title(strcat('ISFS ',int2str(i)))
   set(gca,'Ydir','normal')
   set(gca,'YTick',[1.5,2.5,3.5],'YTickLabel',[1.3,2.7,4])
   set(gca,'XTick',[1.5,2.5,3.5],'XTickLabel',[1.3,2.7,4])
   xlabel('Distance (km)')
   if i == 1
      ylabel('Distance (km)')
   elseif i == 5
      ylabel('Distance (km)')
   end
end

h=colorbar;
colormap(cmap_2011(stn_unique+1,:))
caxis([0 length(stn_unique)])
set(h,'Position', [.7425 .106 .028 .34]);
set(h,'YTick',0.5:1:length(stn_unique),'YTickLabel',nlcd_lnd_id_10(stn_unique+1))


%%

figure(1)
title('USGS Land Use')
ax1 = usamap([lat(lx_min,ly_min),lat(lx_max,ly_max)],[lon(lx_min,ly_max),lon(lx_max,ly_min)]);
surfm(lat,lon,lndUSGS_10)
colormap(cmap_usgs)
h1 = colorbar;
caxis([0 11])
set(h1,'YTick',[0.5:1:length(usgs_lnd_id_10)+1],'YTickLabel',usgs_lnd_id_10)

figure(2)
title('USGS MAXALB')
ax2 = usamap([lat(lx_min,ly_min),lat(lx_max,ly_max)],[lon(lx_min,ly_max),lon(lx_max,ly_min)]);
surfm(lat,lon,albUSGS_10)
colormap(jet)
h2 = colorbar;
caxis([0.4 0.85])

figure(3)
title('NLCD 2011 Land Use')
ax = usamap([lat(lx_min,ly_min),lat(lx_max,ly_max)],[lon(lx_min,ly_max),lon(lx_max,ly_min)]);
surfm(lat,lon,lnd2011_10)
colormap(cmap_2011)
h1 = colorbar;
caxis([0 13])
set(h1,'YTick',[0.5:1:length(nlcd_lnd_id_10)+1],'YTickLabel',nlcd_lnd_id_10)

figure(4)
title('NLCD 2011 MAXALB')
ax2 = usamap([lat(lx_min,ly_min),lat(lx_max,ly_max)],[lon(lx_min,ly_max),lon(lx_max,ly_min)]);
surfm(lat,lon,albNLCD_10)
colormap(jet)
h2 = colorbar;
caxis([0.4 0.85])

figure(4)
title('NLCD 2011 Land Use Modified')
ax = usamap([lat(lx_min,ly_min),lat(lx_max,ly_max)],[lon(lx_min,ly_max),lon(lx_max,ly_min)]);
surfm(lat,lon,fin_lnd_10)
colormap(cmap_2011)
h1 = colorbar;
caxis([0 13])
set(h1,'YTick',[0.5:1:length(nlcd_lnd_id_10)+1],'YTickLabel',nlcd_lnd_id_10)

%%
close all

figure('Position', [100, 100, 900, 800])

subplot(2,2,1)
title('USGS')
ax1 = usamap([lat(lx_min,ly_min),lat(lx_max,ly_max)],[lon(lx_min,ly_max),lon(lx_max,ly_min)]);
surfm(lat,lon,lndUSGS_10)
caxis([0 11])
h1 = colorbar;
colormap(ax1,cmap_usgs)
set(h1,'YTick',[0.5:1:length(usgs_lnd_id_10)+1],'YTickLabel',usgs_lnd_id_10)
mlabel('off')

fl = figm_label('a','ul',14,0.04,0.05,16000,16000);
set(fl,'color','k');

subplot(2,2,2)
title('NLCD 2006')
ax2 = usamap([lat(lx_min,ly_min),lat(lx_max,ly_max)],[lon(lx_min,ly_max),lon(lx_max,ly_min)]);
surfm(lat,lon,lnd2006_10)
caxis([0 13])
colormap(ax2,cmap_2011)
plabel('off')
mlabel('off')

fig_nudge(ax2,0.01,0)

fl = figm_label('b','ul',14,0.04,0.05,16000,16000);
set(fl,'color','k');

subplot(2,2,3)
title('NLCD 2011')
ax3 = usamap([lat(lx_min,ly_min),lat(lx_max,ly_max)],[lon(lx_min,ly_max),lon(lx_max,ly_min)]);
surfm(lat,lon,lnd2011_10)
caxis([0 13])
h3 = colorbar;
colormap(ax3,cmap_2011)
set(h3,'YTick',[0.5:1:length(nlcd_lnd_id_10)+1],'YTickLabel',nlcd_lnd_id_10)

fig_nudge(ax3,-0.037,0.08)

fl = figm_label('c','ul',14,0.04,0.05,16000,16000);
set(fl,'color','k');

subplot(2,2,4)
title('NLCD 2011 Modified')
ax4 = usamap([lat(lx_min,ly_min),lat(lx_max,ly_max)],[lon(lx_min,ly_max),lon(lx_max,ly_min)]);
surfm(lat,lon,fin_lnd_10)
caxis([0 13])
colormap(ax4, cmap_2011)
plabel('off')

fig_nudge(ax4,0.01,0.08)

fl = figm_label('d','ul',14,0.04,0.05,16000,16000);
set(fl,'color','k');

%% Plot lake domain land cover.

figure(1)
subplot(2,2,1)
title('USGS Land Use')
ax = usamap([lat(lx_min,ly_min),lat(lx_max,ly_max)],[lon(lx_min,ly_max),lon(lx_max,ly_min)]);
surfm(lat,lon,lndUSGS)
h1 = colorbar;
set(h1,'YTick',lake_lndUSGS_unique+0.5,'YTickLabel',usgs_lnd_id(lake_lndUSGS_unique))
contourcmap('parula',1)

subplot(2,2,2)
title('NLCD 2006 Land Use')
ax = usamap([lat(lx_min-5,ly_min-5),lat(lx_max-5,ly_max-5)],[lon(lx_min-5,ly_max-5),lon(lx_max-5,ly_min)]);
surfm(lat,lon,lnd2006)
h2 = colorbar;
set(h2,'YTick',lake_lnd2006_unique+0.5,'YTickLabel',nlcd_lnd_id(lake_lnd2006_unique))
contourcmap('parula',1)

subplot(2,2,3)
title('NLCD 2011 Land Use')
ax = usamap([lat(lx_min,ly_min),lat(lx_max,ly_max)],[lon(lx_min,ly_max),lon(lx_max,ly_min)]);
surfm(lat,lon,lnd2011)
h3 = colorbar;
set(h3,'YTick',lake_lnd2011_unique+0.5,'YTickLabel',nlcd_lnd_id(lake_lnd2011_unique))
contourcmap('parula',1)


subplot(2,2,4)
title('NLCD 2011 Land Use Modified')
ax = usamap([lat(lx_min,ly_min),lat(lx_max,ly_max)],[lon(lx_min,ly_max),lon(lx_max,ly_min)]);
surfm(lat,lon,fin_lnd)
h4 = colorbar;
set(h4,'YTick',lake_mod_unique+0.5,'YTickLabel',nlcd_lnd_id(lake_mod_unique))
contourcmap('parula',1)


%% Plot valley domain land cover.

figure(2)
subplot(2,2,1)
title('USGS Land Use')
ax = usamap([lat(vx_min,vy_min),lat(vx_max,vy_max)],[lon(vx_min,vy_max),lon(vx_max,vy_min)]);
surfm(lat,lon,lndUSGS)
h1 = colorbar;
set(h1,'YTick',valley_lndUSGS_unique+0.5,'YTickLabel',usgs_lnd_id(valley_lndUSGS_unique))
contourcmap('parula',1)

subplot(2,2,2)
title('NLCD 2006 Land Use')
ax = usamap([lat(vx_min-5,vy_min-5),lat(vx_max-5,vy_max-5)],[lon(vx_min-5,vy_max-5),lon(vx_max-5,vy_min-5)]);
surfm(lat,lon,lnd2006)
h2 = colorbar;
set(h2,'YTick',valley_lnd2006_unique+0.5,'YTickLabel',nlcd_lnd_id(valley_lnd2006_unique))
contourcmap('parula',1)

subplot(2,2,3)
title('NLCD 2011 Land Use')
ax = usamap([lat(vx_min,vy_min),lat(vx_max,vy_max)],[lon(vx_min,vy_max),lon(vx_max,vy_min)]);
surfm(lat,lon,lnd2011)
h3 = colorbar;
set(h3,'YTick',valley_lnd2011_unique+0.5,'YTickLabel',nlcd_lnd_id(valley_lnd2011_unique))
contourcmap('parula',1)

subplot(2,2,4)
title('NLCD 2011 Land Use Modified')
ax = usamap([lat(vx_min,vy_min),lat(vx_max,vy_max)],[lon(vx_min,vy_max),lon(vx_max,vy_min)]);
surfm(lat,lon,fin_lnd)
h4 = colorbar;
set(h4,'YTick',valley_mod_unique+0.5,'YTickLabel',nlcd_lnd_id(valley_mod_unique))
contourcmap('parula',1)

%% Big Maps

figure(3)
title('NLCD 2011 Terrain Height (m)')
ax = usamap([lat(x_min,y_min),lat(x_max,y_max)],[lon(x_min,y_max),lon(x_max,y_min)]);
%surfm(lat(lx_min:lx_max,ly_min:ly_max),lon(lx_min:lx_max,ly_min:ly_max),lake_mod)
contourfm(lat,lon,hgt,[1000:125:4000])
h4 = colorbar;
% set(h4,'YTick',lake_mod_unique+0.5,'YTickLabel',nlcd_lnd_id(lake_mod_unique))
contourcmap('parula',1)
hold on
contourm(lat,lon,flg,[1],'w')
% hold on
% plotm(lat(365,240),lon(365,240),'w*')

c = colorcube(100);

figure(4)
title('NLCD 2011 Land Use Modified')
ax = usamap([lat(vx_min,vy_min),lat(vx_max,vy_max)],[lon(vx_min,vy_max),lon(vx_max,vy_min)]);
surfm(lat(vx_min:vx_max,vy_min:vy_max),lon(vx_min:vx_max,vy_min:vy_max),valley_mod)
h4 = colorbar;
colormap(c(min(valley_mod_unique):max(valley_mod_unique)-1,:))
set(h4,'YTick',valley_mod_unique+0.5,'YTickLabel',nlcd_lnd_id(valley_mod_unique))

figure(5)
title('NLCD 2011 Land Use Modified')
ax = usamap([lat(vx_min,vy_min),lat(vx_max,vy_max)],[lon(vx_min,vy_max),lon(vx_max,vy_min)]);
surfm(lat(vx_min:vx_max,vy_min:vy_max),lon(vx_min:vx_max,vy_min:vy_max),hgt(vx_min:vx_max,vy_min:vy_max))
hold on
colorbar
plotm(lat(212:240,307),lon(212:240,307),'w-','linewidth',3)
hold on
plotm(lat(230,290:322),lon(230,290:322),'w-','linewidth',3)

figure(6)
title('NLCD 2011 Land Use')
ax = usamap([lat(lx_min,ly_min),lat(lx_max,ly_max)],[lon(lx_min,ly_max),lon(lx_max,ly_min)]);
surfm(lat,lon,lnd2011)
h3 = colorbar;
set(h3,'YTick',lake_lnd2011_unique+0.5,'YTickLabel',nlcd_lnd_id(lake_lnd2011_unique))
contourcmap('parula',1)

figure(7)
title('NLCD 2011 Land Use')
ax = usamap([lat(lx_min,ly_min),lat(lx_max,ly_max)],[lon(lx_min,ly_max),lon(lx_max,ly_min)]);
surfm(lat,lon,fin_lnd)
h3 = colorbar;
set(h3,'YTick',lake_lnd2011_unique+0.5,'YTickLabel',nlcd_lnd_id(lake_lnd2011_unique))
contourcmap('parula',1)
hold on
contourm(lat,lon,hgt,[0:200:5000],'w')
plotm(lat(219,343),lon(219,343),'r.','MarkerSize',15)

flg_diff = flgUSGS-fin_flg;

figure(8)
subplot(1,2,1)
title('Landmask Difference, NLCD 2011 - USGS')
ax = usamap([lat(lx_min,ly_min),lat(lx_max,ly_max)],[lon(lx_min,ly_max),lon(lx_max,ly_min)]);
surfm(lat,lon,flg_diff)
contourcmap('gray',[-1:1])
colorbar()
hold on
contourm(lat,lon,flg,'r','LevelStep',1,'LineWidth',1.5)
contourm(lat,lon,hgt,'k','LevelStep',125,'LineWidth',0.5)

%% Percentage of each land use for each domain/dataset.

for i = 1:length(lake_lnd2011_unique)
   per_lake_lnd2011(i) = length(find(lake_lnd2011==lake_lnd2011_unique(i)))/length(lake_lnd2011(:))
end

for i = 1:length(lake_lnd2006_unique)
   per_lake_lnd2006(i) = length(find(lake_lnd2006==lake_lnd2006_unique(i)))/length(lake_lnd2006(:))
end

for i = 1:length(lake_lndUSGS_unique)
   per_lake_lndUSGS(i) = length(find(lake_lndUSGS==lake_lndUSGS_unique(i)))/length(lake_lndUSGS(:))
end

for i = 1:length(lake_mod_unique)
   per_lake_mod(i) = length(find(lake_mod==lake_mod_unique(i)))/length(lake_mod(:))
end

for i = 1:length(valley_lnd2011_unique)
   per_valley_lnd2011(i) = length(find(valley_lnd2011==valley_lnd2011_unique(i)))/length(valley_lnd2011(:))
end

for i = 1:length(valley_lnd2006_unique)
   per_valley_lnd2006(i) = length(find(valley_lnd2006==valley_lnd2006_unique(i)))/length(valley_lnd2006(:))
end

for i = 1:length(valley_lndUSGS_unique)
   per_valley_lndUSGS(i) = length(find(valley_lndUSGS==valley_lndUSGS_unique(i)))/length(valley_lndUSGS(:))
end

for i = 1:length(valley_mod_unique)
   per_valley_mod(i) = length(find(valley_mod==valley_mod_unique(i)))/length(valley_mod(:))
end

