%% Examples using lidar.nc
url = 'D:\OpenEarthRawData\lidar\oregon_6.nc';

%% read shorelines 0f 1930 and 2002
tic
shoreEast = nc_varget(url, 'shore_east');
shoreNorth = nc_varget(url, 'shore_north');
[lon,lat] = convertCoordinates(shoreEast,shoreNorth,'CS1.code',32610,'CS2.code',4326);
KMLline(lat,lon,'fileName','shoreline_2002.kml','lineColor',[.5 1 1], 'lineWidth',[2]);
shoreEast = nc_varget(url, 'shore_east_1930');
shoreNorth = nc_varget(url, 'shore_north_1930');
[lon,lat] = convertCoordinates(shoreEast,shoreNorth,'CS1.code',32610,'CS2.code',4326);
KMLline(lat,lon,'fileName','shoreline_1930.kml','lineColor',[1 0.5 1], 'lineWidth',[2]);
KMLmerge_files('fileName','shorelines.kmz','sourceFiles',{'shoreline_1930.kml' 'shoreline_2002.kml'})
delete('shoreline_1930.kml');delete('shoreline_2002.kml');
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');
lon = nc_varget(url, 'lon',[time_index-1, id_index-1, 0], [1, 1, length(crossShoreCoordinate)]);
lat = nc_varget(url, 'lat',[time_index-1, id_index-1, 0], [1, 1, length(crossShoreCoordinate)]);
z = nc_varget(url, 'altitude', [time_index-1, id_index-1, 0], [1, 1, length(crossShoreCoordinate)]);
nan_id = lon == -9999 | lat == -9999 | z == -9999 | z == nan;
lon = lon(~nan_id);lat = lat(~nan_id);z = z(~nan_id);
toc

%% make kml file of transect
filename = 'Transect 200';
KMLline(lat,lon,z+10,'fileName',[filename '_1.kml'],'fillColor',jet(5),'fillAlpha',[1 0 1 0 1],'lineColor',[1 1 1],'lineWidth',[1],'lineAlpha',[0.5]);
KMLline([lat(1) lat(1)]',[lon(1) lon(1)]',[0 z(1)+35]','fileName',[filename '_2.kml'],'lineColor',[1 1 1]);
KMLtext(lat(1),lon(1),'Transect 200',z(1)+40,'fileName',[filename '_3.kml'])
KMLmerge_files('fileName',[filename '.kml'],'sourceFiles',{[filename '_1.kml'] [filename '_2.kml'] [filename '_3.kml']});delete([filename '*.kml']);

%% make multiple transect kml files
tic
time_index = find(time == time(1));
crossShoreCoordinate = nc_varget(url, 'cross_shore');
transects = [150:350];

for j = 1 : length(transects)
    id_index = find(id == transects(j));
    lon = nc_varget(url, 'lon',[time_index-1, id_index-1, 0], [1, 1, length(crossShoreCoordinate)]);
    lat = nc_varget(url, 'lat',[time_index-1, id_index-1, 0], [1, 1, length(crossShoreCoordinate)]);
    z = nc_varget(url, 'altitude', [time_index-1, id_index-1, 0], [1, 1, length(crossShoreCoordinate)]);
    nan_id = lon == -9999 | lat == -9999 | z == -9999 | z == nan;
    lon = lon(~nan_id);lat = lat(~nan_id);z = z(~nan_id);
    filename = ['Transect ' num2str(transects(j),'%05.0f')];
    KMLline(lat,lon,z+10,'fileName',[filename '_1.kml'],'fillColor',jet(5),'fillAlpha',[1 0 1 0 1],'lineColor',[1 1 1],'lineWidth',[1],'lineAlpha',[0.5]);
    KMLline([lat(1) lat(1)]',[lon(1) lon(1)]',[0 z(1)+35]','fileName',[filename '_2.kml'],'lineColor',[1 1 1]);
    KMLtext(lat(1),lon(1),filename,z(1)+40,'fileName',[filename '_3.kml']);
    KMLmerge_files('fileName',[filename '.kml'],'sourceFiles',{[filename '_1.kml'] [filename '_2.kml'] [filename '_3.kml']});delete([filename '_*.kml']);
    sourceFiles{j}=[filename '.kml'];
end
KMLmerge_files('fileName',['Transects_' num2str(transects(1),'%05.0f') '-' num2str(transects(end),'%05.0f') '.kmz'],'sourceFiles',sourceFiles);delete(['Transect *.kml']);
toc

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

%% make KML trisurf plot

% filter strangely high data points first
tri = delaunay(lon,lat);
z(tri(any((abs(diff(z(tri),1,2))>15),2),:)) = nan;

% normal triangles
tri = delaunay(lon(~isnan(z)),lat(~isnan(z)));
trisurf(tri,lon(~isnan(z)),lat(~isnan(z)),z(~isnan(z)));shading interp

% filtered triangles
tri2 = delaunay_simplified(lon(~isnan(z)),lat(~isnan(z)),z(~isnan(z)),.25,'maxSize',30000);
trisurf(tri2,lon(~isnan(z)),lat(~isnan(z)),z(~isnan(z)));shading interp

% check triangular structure
trimesh(tri,lon(~isnan(z)),lat(~isnan(z)))

% show in GE 
KMLtrisurf(tri2,lat(~isnan(z)),lon(~isnan(z)),z(~isnan(z)),'fileName','trisurf.kmz')

%% make KML surf plot
% interpolate to grid
[G.lat,G.lon] = meshgrid(linspace(min(lat),max(lat),150),linspace(min(lon),max(lon),200));
G.z = griddata(lat(~isnan(z)),lon(~isnan(z)),z(~isnan(z)),G.lat,G.lon)

% surf(G.lon,G.lat,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

filename = 'surf.kml';
KMLsurf(G.lat,G.lon,G.z,'fileName',[filename '_1.kml'],'zScaleFun',@(z)(z+20)*4,'colorMap',@(m)colormap_cpt('bathymetry_vaklodingen',m),'colorSteps',200,'cLim',[-50 25]);
KMLline(G.lat,G.lon,G.z ,'fileName',[filename '_2.kml'],'zScaleFun',@(z)(z+20)*4,'fill',false);
KMLline(G.lat',G.lon',G.z','fileName',[filename '_3.kml'],'zScaleFun',@(z)(z+20)*4,'fill',false);
KMLmerge_files('fileName',['Surf.kmz'],'sourceFiles',{[filename '_1.kml'],[filename '_2.kml'],[filename '_3.kml']});delete([filename '_*.kml']);

% %% Make rendered png
% tri = delaunay(lon(~isnan(z)),lat(~isnan(z)));
% h = trisurf(tri,lon(~isnan(z)),lat(~isnan(z)),z(~isnan(z)));shading interp
% colormap(colormap_cpt('bathymetry_vaklodingen',200));clim([-50 25]);
% lightangle(90,60);shading interp;material([.7 .3 0.3]);lighting phong
% KMLfig2png(h,'fileName',['Rendered_png.kml'],'levels',[-3 3],'alpha',0.7);

% plot(linspace(-50,25,200),colormap_cpt('bathymetry_vaklodingen',200))
% tri2 = delaunay_simplified(lon(~isnan(z)),lat(~isnan(z)),z(~isnan(z)),2,'maxSize',10000);
% subplot(2,1,1)
% trisurf(tri,lon(~isnan(z)),lat(~isnan(z)),z(~isnan(z)));shading interp
% subplot(2,1,2)
% KMLcontour3(G.lat,G.lon,G.z,'zScaleFun',zScaleFun,'colorMap',@(m)colormap_cpt('bathymetry_vaklodingen',m),'colorSteps',200,'cLim',[-50 25])
