%% Examples using lidar.nc
addpath d:\sonnevi\Projecten\10002364_USGS\code_amy\
url = 'D:\oregon_6.nc';

%% Lookup netcdf info
info = nc_info(url)

%% read and plot 'contours' 1.4 s
tic
contours = nc_varget(url, 'contour');
UCIT_plotLandboundary('Lidar Data US',0,0); hold on;
transectcontours = repmat(nan,3*length(contours),2);
transectcontours(1:3:3*length(contours),1) = contours(1:length(contours),1);transectcontours(2:3:3*length(contours),1) = contours(1:length(contours),2);transectcontours(1:3:3*length(contours),2) = contours(1:length(contours),4);transectcontours(2:3:3*length(contours),2) = contours(1:length(contours),3);
plot(transectcontours(:,1),transectcontours(:,2),'b');
axis equal;
toc

%% read data of single transect 0.2 s
tic
id = nc_varget(url, 'id');
time = nc_varget(url, 'time');
id_index = find(id == 200);
time_index = find(time == time(1));
crossShoreCoordinate = nc_varget(url, 'cross_shore');
d.xe = nc_varget(url, 'cross_shore_distance',[time_index-1, id_index-1, 0], [1, 1, length(crossShoreCoordinate)]);
d.ze = nc_varget(url, 'altitude', [time_index-1, id_index-1, 0], [1, 1, length(crossShoreCoordinate)]);
toc

%% Plot the transect 0.17 s
tic
figure;
plot(d.xe,d.ze,'xb');hold on;
xlim([min(d.xe(d.xe~=-9999)) max(d.xe(d.xe~=-9999))]);
toc

%% read data and metadata of single transect 1.7 s
tic
d = readLidarDataNetcdf(url, 200, time(1));
toc

%% Add USGS metadata 0.14s
tic
regression = nc_varget(url, 'regression',[time_index-1, id_index-1, 0], [1, 1, length(crossShoreCoordinate)]);
plot(d.xe(~isnan(regression)),d.ze(~isnan(regression)),'or');
shorePos = nc_varget(url, 'shorepos', [id_index-1], [1]);
MHW = nc_varget(url, 'mean_high_water', [id_index-1], [1]);
plot(shorePos, MHW,'om','markersize',10);
line([min(d.xe(d.xe~=-9999)) max(d.xe(d.xe~=-9999))],[MHW MHW],'color','k');
toc

% %% plot id's 60 s
% tic
% ids = nc_varget(url, 'id');
% h_i = cellfun(@text,{contours(:,1)},{contours(:,2)}, {num2str(ids)},'UniformOutput',false);
% toc

%% load slightly filtered patch of raw data of multiple transects in UTM 2
tic
x = nc_varget(url, 'x', [0,150,0], [1,200,400],[1 1 10]);
y = nc_varget(url, 'y', [0,150,0], [1,200,400],[1 1 10]);
z = nc_varget(url, 'altitude', [0,150,0], [1,200,400],[1 1 10]);
id = x == -9999 | y == -9999 | z == -9999 | z == nan;
x = x(~id);y = y(~id);z = z(~id);
toc

%% load raw data of all transects (with thinning of 10) 162 s
tic
crossShoreCoordinate = nc_varget(url, 'cross_shore');
x = nc_varget(url, 'x', [0,0,0], [2,251,462], [1 100 10]);
y = nc_varget(url, 'y', [0,0,0], [2,251,462], [1 100 10]);
z = nc_varget(url, 'altitude', [0,0,0], [2,251,462], [1 100 10]);
toc


%% plot them 'USGS style' (all transects: 43 s)
tic
load colormapMHWjump20
UCIT_plotLandboundary('Lidar Data US',1,1); hold on;
cdots_amy(x,y,z,(19 * .5/8),(44 * .5/8),cmapMHWjump20); axis equal;
toc

%% add shoreline position
tic
shoreEast = nc_varget(url, 'shore_east');
shoreNorth = nc_varget(url, 'shore_north');
plot(shoreEast,shoreNorth,'ow','marker','o','markerfacecolor','w','markeredgecolor','k','markersize',6);
toc

%% use different plotmethod
tic
x = x(z < 50); y = y(z < 50); z = z(z < 50); 
plotc(x,y,z);view(3);
toc

%% make triangular grid and plot
tri = delaunay(x,y);
trisurf(tri,x(~isnan(z)),y(~isnan(z)),z(~isnan(z)));
shading interp; view(2);colorbar;axis equal;colormap(colormap_cpt('bathymetry_vaklodingen',200));clim([-50 25]);

% filter triangles
% nan_id = any((abs(diff(z(tri),1,2))>15),2);
% z(tri(nan_id,:)) = nan;
diff_array = abs(diff(z(tri),1,2))>5;
z(tri(diff_array(:,1)&sum(diff_array,2)~=2,1)) = nan;
z(tri(sum(diff_array,2)==2,2)) = nan;
z(tri(diff_array(:,2)&sum(diff_array,2)~=2,3)) = nan;

% trimesh(tri,x(~isnan(z)),y(~isnan(z)));

tri2 = delaunay_simplified(x(~isnan(z)),y(~isnan(z)),z(~isnan(z)),.25,'maxSize',30000);
trisurf(tri2,x(~isnan(z)),y(~isnan(z)),z(~isnan(z)));
shading interp; view(2);colorbar;axis equal;colormap(colormap_cpt('bathymetry_vaklodingen',200));clim([-50 25]);

%% interpolate to rectangular grid and plot
tic
[G.x,G.y] = meshgrid(linspace(min(x),max(x),1000),linspace(min(y),max(y),3000));
G.z = griddata(x(~isnan(z)),y(~isnan(z)),z(~isnan(z)),G.x,G.y)

surf(G.x,G.y,G.z);shading interp; axis equal;view(2);
colormap(colormap_cpt('bathymetry_vaklodingen',200));clim([-50 25]);
lightangle(90,60);shading interp;material([.7 .3 0.3]);lighting phong
toc