%% 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.

lat1 = double(ncread('../output/nlcd2011-dec31/geo_em.d01.nc','XLAT_M'));
lon1 = double(ncread('../output/nlcd2011-dec31/geo_em.d01.nc','XLONG_M'));
hgt1 = double(ncread('../output/nlcd2011-dec31/geo_em.d01.nc','HGT_M'));
flg1 = double(ncread('../output/nlcd2011-dec31/geo_em.d01.nc','LANDMASK'));

lat2 = double(ncread('../output/nlcd2011-dec31/geo_em.d02.nc','XLAT_M'));
lon2 = double(ncread('../output/nlcd2011-dec31/geo_em.d02.nc','XLONG_M'));
hgt2 = double(ncread('../output/nlcd2011-dec31/geo_em.d02.nc','HGT_M'));
flg2 = double(ncread('../output/nlcd2011-dec31/geo_em.d02.nc','LANDMASK'));

lat3 = double(ncread('../output/nlcd2011-dec31/geo_em.d03.nc','XLAT_M'));
lon3 = double(ncread('../output/nlcd2011-dec31/geo_em.d03.nc','XLONG_M'));
hgt3 = double(ncread('../output/nlcd2011-dec31/geo_em.d03.nc','HGT_M'));
flg3 = double(ncread('../output/nlcd2011-dec31/geo_em.d03.nc','LANDMASK'));

s1 = size(lat1);
s2 = size(lat2);
s3 = size(lat3);

 lat1lim = [min(lat1(:)) max(lat1(:))];
 lon1lim = [min(lon1(:)) max(lon1(:))];

%lat1lim = [35 45];
%lon1lim = [-121 -103];

figure('Position',[100 100 500 400])

ax = usamap(lat1lim,lon1lim);
set(ax, 'Visible', 'off')
latlim = getm(ax, 'MapLatLimit');
lonlim = getm(ax, 'MapLonLimit');
states = shaperead('usastatehi','UseGeoCoords', true, 'BoundingBox', [lonlim', latlim']);
geoshow(ax, states, 'FaceColor','none')

hold on
sub1 = plotm(lat1(1:s1(1),1),lon1(1:s1(1),1),'k-','LineWidth',2)
hold on
sub1 = plotm(lat1(1:s1(1),s1(2)),lon1(1:s1(1),s1(2)),'k-','LineWidth',2)
hold on
sub1 = plotm(lat1(1,1:s1(2)),lon1(1,1:s1(2)),'k-','LineWidth',2)
hold on
sub1 = plotm(lat1(s1(1),1:s1(2)),lon1(s1(1),1:s1(2)),'k-','LineWidth',2)
hold on

sub2 = plotm(lat2(1:s2(1),1),lon2(1:s2(1),1),'k-','LineWidth',2)
hold on
sub2 = plotm(lat2(1:s2(1),s2(2)),lon2(1:s2(1),s2(2)),'k-','LineWidth',2)
hold on
sub2 = plotm(lat2(1,1:s2(2)),lon2(1,1:s2(2)),'k-','LineWidth',2)
hold on
sub2 = plotm(lat2(s2(1),1:s2(2)),lon2(s2(1),1:s2(2)),'k-','LineWidth',2)
hold on

sub2 = plotm(lat3(1:s3(1),1),lon3(1:s3(1),1),'k-','LineWidth',2)
hold on
sub2 = plotm(lat3(1:s3(1),s3(2)),lon3(1:s3(1),s3(2)),'k-','LineWidth',2)
hold on
sub2 = plotm(lat3(1,1:s3(2)),lon3(1,1:s3(2)),'k-','LineWidth',2)
hold on
sub2 = plotm(lat3(s3(1),1:s3(2)),lon3(s3(1),1:s3(2)),'k-','LineWidth',2)
hold on

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

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

contourm(lat3(lx_min:lx_max,ly_min:ly_max),lon3(lx_min:lx_max,ly_min:ly_max),flg3(lx_min:lx_max,ly_min:ly_max),'k-')
hold on

v_sub = plotm(lat3(vx_min:vx_max,vy_min),lon3(vx_min:vx_max,vy_min),'k-')
hold on
v_sub = plotm(lat3(vx_min:vx_max,vy_max),lon3(vx_min:vx_max,vy_max),'k-')
hold on
v_sub = plotm(lat3(vx_min,vy_min:vy_max),lon3(vx_min,vy_min:vy_max),'k-')
hold on
v_sub = plotm(lat3(vx_max,vy_min:vy_max),lon3(vx_max,vy_min:vy_max),'k-')

hold on
l_sub = plotm(lat3(lx_min:lx_max,ly_min),lon3(lx_min:lx_max,ly_min),'k-')
hold on
l_sub = plotm(lat3(lx_min:lx_max,ly_max),lon3(lx_min:lx_max,ly_max),'k-')
hold on
l_sub = plotm(lat3(lx_min,ly_min:ly_max),lon3(lx_min,ly_min:ly_max),'k-')
hold on
l_sub = plotm(lat3(lx_max,ly_min:ly_max),lon3(lx_max,ly_min:ly_max),'k-')

textm(50,-127,'(a)','FontSize',14)
textm(max(lat1(:))-0.75,max(lon1(:))-2.5,'D1')
textm(max(lat2(:))-0.75,max(lon2(:))-3.25,'D2')
textm(max(lat3(:))-0.75,max(lon3(:))-3,'D3')