% OPENS A LAND COVER FILE IN NETCDF AND CREATES THE FORT.13 ATTRIBUTES

filelocation = 'C:\Users\jgonzalez\Downloads\prlandcover_ncint.nc';

% GET COORDINATE LIMITS

north = ncreadatt(filelocation,'GDAL_Geographics','Northernmost_Northing');
south = ncreadatt(filelocation,'GDAL_Geographics','Southernmost_Northing');
east = ncreadatt(filelocation,'GDAL_Geographics','Easternmost_Easting');
west = ncreadatt(filelocation,'GDAL_Geographics','Westernmost_Easting');

dx = 0.000305637;
dy = 0.000426472;

nx = 8427;
ny = 4613;

% MAKE MESH FOR DATA

x = linspace(east,west,nx);
y = linspace(south,north,ny);

[XX,YY] = meshgrid(x,y);
XX = fliplr(XX);
YY = fliplr(YY);

% GET LANDCOVER DATA

landcoverdata = ncread(filelocation,'Band1');
landcoverdata = rot90(landcoverdata,1);

%landcoverdata(find(isnan(landcoverdata) == 1)) = 0;

% ASSIGN Z0 TO EACH LAND USE CATEGORY, FROM BUNYA 2010

% 11 = open water
landcoverdata(find(landcoverdata == 11)) = 0.001;

% 21 = low residential
landcoverdata(find(landcoverdata == 21)) = 0.330;

% 22 = medium residential
landcoverdata(find(landcoverdata == 22)) = 0.500;

% 23 = commercial
landcoverdata(find(landcoverdata == 23)) = 0.390;

% 24 = high residential
landcoverdata(find(landcoverdata == 24)) = 0.500;

% 31 = bare rock
landcoverdata(find(landcoverdata == 31)) = 0.090;

% 32 = gravel pit
landcoverdata(find(landcoverdata == 32)) = 0.180;

% 41 = deciduous forest
landcoverdata(find(landcoverdata == 41)) = 0.650;

% 42 = evergreen forest
landcoverdata(find(landcoverdata == 42)) = 0.720;

% 43 = mixed forest
landcoverdata(find(landcoverdata == 43)) = 0.710;

% 52 = shrub
landcoverdata(find(landcoverdata == 52)) = 0.120;

% 71 = grassland
landcoverdata(find(landcoverdata == 71)) = 0.040;

% 81 = pasture
landcoverdata(find(landcoverdata == 81)) = 0.060;

% 82 = cultivated crops
landcoverdata(find(landcoverdata == 82)) = 0.060;

% 90 = woody wetlands
landcoverdata(find(landcoverdata == 90)) = 0.550;

% 91 = palaustrine forested wetland
landcoverdata(find(landcoverdata == 91)) = 0.550;

% 92 = herbaceous wetland
landcoverdata(find(landcoverdata == 92)) = 0.110;

% 93 = estuarine forested wetland
landcoverdata(find(landcoverdata == 93)) = 0.550;

% 94 = estuarine shrub wetland
landcoverdata(find(landcoverdata == 94)) = 0.120;

% 95 = emergent herbaceous wetland
landcoverdata(find(landcoverdata == 95)) = 0.110;

% 96 = palaustrine emergent wetland
landcoverdata(find(landcoverdata == 96)) = 0.110;

% 97 = estuarine emergent wetland
landcoverdata(find(landcoverdata == 97)) = 0.110;

% anything else
landcoverdata(find(landcoverdata > 1)) = 0.001;

%%

f1 = figure(1)
pcolor(landcoverdata)
shading flat;
colorbar

caxis([0 0.8]);
title('PRVI roughness length');
c = colorbar;
ylabel(c,'z0 (m)')

set(f1, 'Units', 'inches');
set(f1,'position',[0 0 8 4]);
set(f1, 'PaperUnits', 'inches', 'PaperSize', [8 4], 'PaperPosition',[0 0 8 4], 'PaperPositionMode', 'manual');
print(f1,'-r150','-dpng','prvi_roughness1.png')


f2 = figure(2)
pcolor(landcoverdata(1500:2100,500:7000))
shading flat;
colorbar

caxis([0 0.8]);
title('PRVI roughness length');
c = colorbar;
ylabel(c,'z0 (m)')

set(f2, 'Units', 'inches');
set(f2,'position',[0 0 10 2]);
set(f2, 'PaperUnits', 'inches', 'PaperSize', [10 2], 'PaperPosition',[0 0 10 2], 'PaperPositionMode', 'manual');
print(f2,'-r150','-dpng','prvi_roughness2.png')


%% PROCESS ST. CROIX LAND COVER DATA 

location = 'I:\GRIDS\landcover\usvi\stcroix_wgs84.nc';

% land cover data
landcoverdata = ncread(location,'z');

% dimensions x,y
dim = ncread(location,'dimension');

% range in x,y
x_range = ncread(location,'x_range');
y_range = ncread(location,'y_range');

% ASSIGN Z0 TO EACH LAND USE CATEGORY, FROM BUNYA 2010

% 0,21 = open water
landcoverdata(find(landcoverdata == 0)) = 0.001;
landcoverdata(find(landcoverdata == 21)) = 0.001;

% 5 = low residential
landcoverdata(find(landcoverdata == 5)) = 0.330;

% 6 = cultivated crops
landcoverdata(find(landcoverdata == 6)) = 0.060;

% 7 = pasture
landcoverdata(find(landcoverdata == 7)) = 0.060;

% 8 = grassland / herbaceous
landcoverdata(find(landcoverdata == 8)) = 0.040;

% 9 = deciduous forest
landcoverdata(find(landcoverdata == 9)) = 0.650;

% 10 = evergreen forest
landcoverdata(find(landcoverdata == 10)) = 0.720;

% 11 = mixed forest
landcoverdata(find(landcoverdata == 11)) = 0.710;

% 12 = shrub
landcoverdata(find(landcoverdata == 12)) = 0.120;

% 13 = palaustrine forested wetland
landcoverdata(find(landcoverdata == 13)) = 0.550;

% 14 = palaustrine shrub wetland
landcoverdata(find(landcoverdata == 13)) = 0.120;

% 15 = palaustrine emergent wetland
landcoverdata(find(landcoverdata == 15)) = 0.110;

% 16 = estuarine forested wetland
landcoverdata(find(landcoverdata == 16)) = 0.550;

% 17 = estuarine shrub wetland
landcoverdata(find(landcoverdata == 17)) = 0.120;

% 18 = estuarine emergent wetland
landcoverdata(find(landcoverdata == 18)) = 0.110;

% anything else
landcoverdata(find(landcoverdata > 1)) = 0.001;

% DO PLOT

z2 = rot90( reshape(landcoverdata,dim(1),dim(2)) , 1);

f1 = figure(1)
pcolor(z2)
shading flat;
colorbar

caxis([0 0.8]);
title('St. Croix roughness length');
c = colorbar;
ylabel(c,'z0 (m)')

set(f1, 'Units', 'inches');
set(f1,'position',[0 0 8 4]);
set(f1, 'PaperUnits', 'inches', 'PaperSize', [8 4], 'PaperPosition',[0 0 8 4], 'PaperPositionMode', 'manual');
print(f1,'-r150','-dpng','stcroix_roughness1.png')


%% PROCESS ST. THOMAS LAND COVER DATA 

location = 'I:\GRIDS\landcover\usvi\stthomas_wgs84.nc';

% land cover data
landcoverdata = ncread(location,'z');

% dimensions x,y
dim = ncread(location,'dimension');

% range in x,y
x_range = ncread(location,'x_range');
y_range = ncread(location,'y_range');

% ASSIGN Z0 TO EACH LAND USE CATEGORY, FROM BUNYA 2010

% 0,21 = open water
landcoverdata(find(landcoverdata == 0)) = 0.001;
landcoverdata(find(landcoverdata == 21)) = 0.001;

% 5 = low residential
landcoverdata(find(landcoverdata == 5)) = 0.330;

% 6 = cultivated crops
landcoverdata(find(landcoverdata == 6)) = 0.060;

% 7 = pasture
landcoverdata(find(landcoverdata == 7)) = 0.060;

% 8 = grassland / herbaceous
landcoverdata(find(landcoverdata == 8)) = 0.040;

% 9 = deciduous forest
landcoverdata(find(landcoverdata == 9)) = 0.650;

% 10 = evergreen forest
landcoverdata(find(landcoverdata == 10)) = 0.720;

% 11 = mixed forest
landcoverdata(find(landcoverdata == 11)) = 0.710;

% 12 = shrub
landcoverdata(find(landcoverdata == 12)) = 0.120;

% 13 = palaustrine forested wetland
landcoverdata(find(landcoverdata == 13)) = 0.550;

% 14 = palaustrine shrub wetland
landcoverdata(find(landcoverdata == 13)) = 0.120;

% 15 = palaustrine emergent wetland
landcoverdata(find(landcoverdata == 15)) = 0.110;

% 16 = estuarine forested wetland
landcoverdata(find(landcoverdata == 16)) = 0.550;

% 17 = estuarine shrub wetland
landcoverdata(find(landcoverdata == 17)) = 0.120;

% 18 = estuarine emergent wetland
landcoverdata(find(landcoverdata == 18)) = 0.110;

% anything else
landcoverdata(find(landcoverdata > 1)) = 0.001;

z2 = rot90( reshape(landcoverdata,dim(1),dim(2)) , 1);

f1 = figure(1)
pcolor(z2)
shading flat;
colorbar

caxis([0 0.8]);
title('St. Thomas roughness length');
c = colorbar;
ylabel(c,'z0 (m)')

set(f1, 'Units', 'inches');
set(f1,'position',[0 0 8 4]);
set(f1, 'PaperUnits', 'inches', 'PaperSize', [8 4], 'PaperPosition',[0 0 8 4], 'PaperPositionMode', 'manual');
print(f1,'-r150','-dpng','stthomas_roughness1.png')

%% PROCESS ST. JOHN LAND COVER DATA 

location = 'I:\GRIDS\landcover\usvi\stjohn_wgs84.nc';

% land cover data
landcoverdata = ncread(location,'z');

% dimensions x,y
dim = ncread(location,'dimension');

% range in x,y
x_range = ncread(location,'x_range');
y_range = ncread(location,'y_range');

% ASSIGN Z0 TO EACH LAND USE CATEGORY, FROM BUNYA 2010

% 0,21 = open water
landcoverdata(find(landcoverdata == 0)) = 0.001;
landcoverdata(find(landcoverdata == 21)) = 0.001;

% 5 = low residential
landcoverdata(find(landcoverdata == 5)) = 0.330;

% 6 = cultivated crops
landcoverdata(find(landcoverdata == 6)) = 0.060;

% 7 = pasture
landcoverdata(find(landcoverdata == 7)) = 0.060;

% 8 = grassland / herbaceous
landcoverdata(find(landcoverdata == 8)) = 0.040;

% 9 = deciduous forest
landcoverdata(find(landcoverdata == 9)) = 0.650;

% 10 = evergreen forest
landcoverdata(find(landcoverdata == 10)) = 0.720;

% 11 = mixed forest
landcoverdata(find(landcoverdata == 11)) = 0.710;

% 12 = shrub
landcoverdata(find(landcoverdata == 12)) = 0.120;

% 13 = palaustrine forested wetland
landcoverdata(find(landcoverdata == 13)) = 0.550;

% 14 = palaustrine shrub wetland
landcoverdata(find(landcoverdata == 13)) = 0.120;

% 15 = palaustrine emergent wetland
landcoverdata(find(landcoverdata == 15)) = 0.110;

% 16 = estuarine forested wetland
landcoverdata(find(landcoverdata == 16)) = 0.550;

% 17 = estuarine shrub wetland
landcoverdata(find(landcoverdata == 17)) = 0.120;

% 18 = estuarine emergent wetland
landcoverdata(find(landcoverdata == 18)) = 0.110;

% anything else
landcoverdata(find(landcoverdata > 1)) = 0.001;

z2 = rot90( reshape(landcoverdata,dim(1),dim(2)) , 1);

f1 = figure(1)
pcolor(z2)
shading flat;
colorbar

caxis([0 0.8]);
title('St. John roughness length');
c = colorbar;
ylabel(c,'z0 (m)')

set(f1, 'Units', 'inches');
set(f1,'position',[0 0 8 4]);
set(f1, 'PaperUnits', 'inches', 'PaperSize', [8 4], 'PaperPosition',[0 0 8 4], 'PaperPositionMode', 'manual');
print(f1,'-r150','-dpng','stjohn_roughness1.png')
