%% script to process JERICO sumemrschool walk and quad data into OGC standard
%  takes order 1 min (matlab scatter print is slowest)

% Running this script requires a check-out of the matlab OpenEarthTools:
% https://publicwiki.deltares.nl/display/OET/Join+OpenEarth
% https://svn.oss.deltares.nl/repos/openearthtools/trunk/matlab/

tic;clear all;close all

OPT.clim = [-4.5 4.5];

process = [1 2 3];

%% load quad
type = 'quad';
load(['..\raw\XYZ_pp_',type,'20140618.mat']);
D(1).X = XYZ.X;
D(1).Y = XYZ.Y;
D(1).Z = XYZ.Z;
D(1).T = XYZ.T;
D(1).name = type;
clear XYZ

%% load walk
type = 'walk';
load(['..\raw\XYZ_pp_',type,'20140618.mat']);
D(2).X = XYZ.X;
D(2).Y = XYZ.Y;
D(2).Z = XYZ.Z;
D(2).T = XYZ.T;
D(2).name = type;
clear XYZ

%% load jetski
type = 'jski';
XYZ = gpsload_corrected(['..\raw\XYZ_jetski_RD_NAP.txt']);
D(3).X = XYZ.x;
D(3).Y = XYZ.y;
D(3).Z = XYZ.z;
D(3).T = [];
D(3).name = type;
clear XYZ

%% convert coordinates
for i=process
[D(i).lon,D(i).lat]=convertCoordinates(D(i).X,D(i).Y,'CS1.code',28992,'CS2.code',4326);
end

%% obtain context data: aerial image towards data
ax = [4.175 4.21 52.045 52.06];
[url,WMS] = wms('server','http://gdsc.nlr.nl/wms/lufo2005?','axis',ax([1 3 2 4]),'styles',1);
if ~exist(['aerial',WMS.ext])
urlwrite(url,['aerial',WMS.ext]); % only once
end
C = nc2struct('D:/opendap.deltares.nl\thredds\dodsC\opendap\deltares\landboundaries\holland_fillable.nc');

%% plot (lat,lon) and context info
close
[IM,map,alpha] = imread(['aerial',WMS.ext]);
image(WMS.x,WMS.y,IM)
colormap(map)
tickmap('ll');grid on;
set(gca,'ydir','normal')
grid on
hold on
for i=process
   plot(D(i).lon,D(i).lat,'DisplayName',D(i).name,'color',[i/3 (3-i)/3 0])
end
axislat
axis(ax); %axis tight
plot(C.lon,C.lat,'color',[.5 .5 .5],'DisplayName','coast')
grid on
legend('location','NorthWest')
print2screensize(['overview.png'])

%% scatter
close
for i=process
   plot(D(i).lon,D(i).lat,'DisplayName',[D(i).name,' track'],'color',[i/3 (3-i)/3 0])
   hold on
   scatter(D(i).lon,D(i).lat,20,D(i).Z,'HandleVisibility','off')
end
axislat
axis(ax); %axis tight
plot(C.lon,C.lat,'color',[.5 .5 .5],'DisplayName','coast')
grid on
clim(OPT.clim)
colorbarwithvtext('DGPS elevation [m]')
legend('location','NorthWest')
%print2screensize('overview_scatter.png')% too slow, Matlab is rubbish

%% plot Google Earth: data towards aerial image
for i=process
KMLline(D(i).lat,D(i).lon,'fileName',[D(i).name,'_line.kml'],...
    'lineColor',[i/3 (3-i)/3 0],...
    'kmlName',D(i).name)
    if ~isempty(D(i).T)
        pngcolorbar = KMLscatter(D(i).lat,D(i).lon,D(i).Z,'fileName',[D(i).name,'_time.kml'],...
        'kmlName',D(i).name,...
        'name',cellstr((num2str(D(i).Z'))),...
        'timeIn',D(i).T,...
        'CBcLim',OPT.clim,...
        'CBcolorSteps',18,...
        'CBcolorTitle',['DGPS elevation [m]'],...
        'colorbar',i==1,... % only one colorbar
        'timeOut',repmat(D(i).T(end),size(D(i).T)));
    else
        pngnames = KMLscatter(D(i).lat,D(i).lon,D(i).Z,'fileName',[D(i).name,'_time.kml'],...
        'kmlName',D(i).name,...
        'CBcolorSteps',18,...
        'CBcLim',OPT.clim,...        
        'name',cellstr((num2str(D(i).Z(:)))),...
        'colorbar',0);
    end
end

%% make overall kml
% copy logo file manually into kmz (is simply zip file)
KMLlogo('jerico_logo.png','fileName','jerico_logo.kml','kmlName','JERICO logo');

KMLmerge_files('sourceFiles',{'walk_line.kml','walk_time.kml',...
                              'quad_line.kml','quad_time.kml',...
                              'jski_line.kml','jski_time.kml','jerico_logo.kml',pngcolorbar{:}},...
                'fileName','JERICO 18th - ShoreMonitoring.kmz',...
                'description','JERICO summerschool Zandmotor practical: https://publicwiki.deltares.nl/display/OET/JERICO, version $HeadURL$ $Id$')
toc