%% PLOT 3D topography and land use and land cover data.

clc 
clear 
close all 

%% read in data 
addpath('/Users/yang560/Dropbox/PNNL/Matlab_functions/utilities/')

file = '../12_SGP_POLAR/zhao_dem/geo_em.d01.nc.1km';
lat2d     = ncread(file, 'XLAT_M');
lon2d     = ncread(file,'XLONG_M');
lulc = ncread(file, 'LU_INDEX');
hgt  = ncread(file, 'HGT_M');

[aspect,slope,gradN,gradE] = gradientm(lat2d,lon2d,hgt);

lat2d = double (lat2d);
lon2d = double(lon2d);

lulc = double (lulc);
hgt  = double (hgt);
%% plot 3d figure 

stamp_lat = [36.881, 35.615, 36.841, 36.605, 36.431, 37.1509, 36.819, 36.9255, 37.0694, 35.8615, 36.1166, 36.3106,  35.8797, 36.37354, 36.31937 ,37.133];
stamp_lon = [-98.285,-96.065 ,-96.427 ,-97.485 , -98.284, -98.362, -97.8199, -97.0817,  -96.7606, -97.0695, -97.5112, -97.928, -98.1728, -97.06905, -96.76203,  -97.266];

lhsh_lat = [36.819,37.0694 , 35.8615, 36.1166, 35.615, 37.150, 36.9255,36.3106 , 35.8797];
lhsh_lon = [-97.8199, -96.7606,-97.0695,-97.5112,-96.065, -98.362, -97.0817,-97.928, -98.1728];

make_it_tight = true;
subplot = @(m,n,p) subtightplot (m, n, p, [0.065 0.045], [0.12 0.025], [0.055 0.13]);
if ~make_it_tight,  clear subplot;  end

figure(1)
set(gcf,'PaperPositionMode', 'auto')
set(gcf,'Position', [100 100 1000 600])
%------------------------------------------------------
subplot(2,2,1)
surf(lon2d, lat2d, hgt); 
shading interp; view(20,50); view(2);
xlabel('Longitude')
ylabel('Latitude')
zlabel('Surface Elevation (m)')
colormap (jet)
h = colorbar; 
% h.Label.String = ['DEM (m)'];
title(h, 'DEM(m)')

title_left('a) DEM ')
zlim([0 800]);
caxis ([100 700])
ylim([35.1, 38.1])
xlim([-99.3,-95.6])
% text(-96, 37.5, 780, '1km')
%------------------------------------------------------
subplot(2,2,2)
surf(lon2d, lat2d, hgt,aspect); 
shading interp; view(2); 
xlabel('Longitude')
% ylabel('Latitude')
zlabel('Surface Elevation (m)')
colormap (jet)
h = colorbar;
% h.Label.String = ['aspect angle'];

title(h, 'aspect angle')
title_left('b) aspect angle')
zlim([0 800]);
ylim([35.1, 38.1])
xlim([-99.3,-95.6])
caxis ([0 360])
% text(-96, 37.5, 780, '1km')
%------------------------------------------------------
subplot(2,2,3)
surf(lon2d, lat2d, hgt,slope); 
shading interp; view(2); 
xlabel('Longitude')
ylabel('Latitude')
zlabel('Surface Elevation (m)')
colormap (jet)
h = colorbar;
% h.Label.String = ['slope angle'];
title(h, 'slope angle')
title_left('c) slope angle')
zlim([0 800]);
ylim([35.1, 38.1])
xlim([-99.3,-95.6])
caxis([0, 3])
%------------------------------------------------------
subplot(2,2,4)
% surf(lon2d, lat2d, hgt,lulc); 
surf(lon2d, lat2d, lulc); hold on;
alpha 0.7;
plot(stamp_lon, stamp_lat, 'KO','Linewidth', 3, 'MarkerSize', 10); hold on;
plot(lhsh_lon, lhsh_lat, 'KX', 'Linewidth', 3,'MarkerSize', 10); hold on;

shading interp; view(2); 
xlabel('Longitude')
% ylabel('Latitude')
zlabel('Surface Elevation (m)')
colormap (jet(17))
h = colorbar;
% h.Label.String = ['land ocver type'];
title(h, 'land cover type')
caxis([.5,17.5]) % value from the exact bottom to the exact top of the colorbar

h.Ticks = linspace(1,17,17);
h.TickLabels = {'Evergreen Needle Forest','Evergreen Broadleaf Forest','Deciduous Needleleaf Forest',...
                'Deciduous Broadleaf Forest','Mixed Forests','Closed Shrublands',...
    'Open Shrublands','Woody Savannas','Savannas',...
    'Grasslands','Permanent wetlands','Croplands',...
    'Urban and Built-Up','cropland/natural vegetation mosaic','Snow and Ice',...
    'Barren or Sparsely Vegetated','Water'};

title_left('d) land cover type')
zlim([0 800]);
ylim([35.1, 38.1])
xlim([-99.3,-95.6])

% ----- add locations of LH/SH + STAMP polygons-----

print('-depsc', 'dem_1km_4plt');


