% takes ca. 45 min for 800 png over 38 tiles

OPT.rijkswaterstaat_dataset_name = 'DienstZeeland'; % jarkus, % vaklodingen

close all
%% get overview
nc_dir  = ['D:\opendap.deltares.nl\thredds\rijkswaterstaat\',OPT.rijkswaterstaat_dataset_name];
[CATALOG.urlPath, CATALOG.x, CATALOG.y,bb.x,bb.y] = grid_orth_getmapinfofromdataset(nc_dir);

CATALOG.urlPath = cellfun(@(x) strrep(x,['http://opendap.deltares.nl/thredds/dodsC/opendap/rijkswaterstaat/',OPT.rijkswaterstaat_dataset_name],nc_dir),CATALOG.urlPath,'UniformOutput',0);

%% get ldb for reference
L.x = nc_varget('http://opendap.deltares.nl/thredds/dodsC/opendap/deltares/landboundaries/holland_fillable.nc','x');
L.y = nc_varget('http://opendap.deltares.nl/thredds/dodsC/opendap/deltares/landboundaries/holland_fillable.nc','y');

tic

%% TEST 1: plot per tile per availble timestep

   AX = subplot_meshgrid(2,1,.05,.05,[nan nan],nan);
   
   for i=1:length(CATALOG.urlPath)
   
      ncfile = [CATALOG.urlPath{i}];
      x = nc_varget (ncfile,'x');
      y = nc_varget (ncfile,'y');
      t = nc_cf_time(ncfile,'time');
   
      BB.x = [x(1) x(end) x(end) x(1) x(1)];
      BB.y = [y(end) y(end) y(1) y(1) y(end)];
   
      for it=1:length(t)
      
         %% get z
         z = nc_varget (ncfile,'z',[(it-1) 0 0],[1 -1 -1]);
        
         %% plot overview
         axes(AX(1));cla
         plot(poly_join(bb.x),poly_join(bb.y),'color',[.5 .5 .5])
         hold on
         plot(BB.x,BB.y,'r','linewidth',2)
         plot(L.x,L.y,'k')
         axis equal
   
         dx = range(poly_join(bb.x))./50;
         dy = range(poly_join(bb.y))./50;
         axis([min(poly_join(bb.x))-dx max(poly_join(bb.x))+dx ...
               min(poly_join(bb.y))-dy max(poly_join(bb.y))+dy])
   
         axis off
         
         %% plot data
         axes(AX(2));cla   
         pcolorcorcen(x,y,z);
         title([mktex(filename(ncfile)),'  -  ',datestr(t(it),'yyyy-mm-dd'),' (',num2str(it),'/',num2str(length(t)),')'])
         axis equal
         grid on
         axis([x(1) x(end) y(end) y(1)])
         set(gca,'xtick',[x(  1):2.5e3:x(end)])
         set(gca,'ytick',[y(end):2.5e3:y(  1)])
         tickmap('xy','texttype','text')
         clim([-25 0])
         colorbarwithhtext('z [m]','horiz')
         hold on
         plot(L.x,L.y,'k')
        
         print2screensizeoverwrite([filepathstrname(CATALOG.urlPath{i}),'_',num2str(it),'_',datestr(t(it),'yyyy-mm-dd'),'.png'])
         
      end
      
   end
   toc
   tic

%% TEST 2: plot per period for all tiles (UCIT in batch mode for whole area instead of polygon)
   
   OPT.axis = [min(poly_join(bb.x)) max(poly_join(bb.x)) ...
               min(poly_join(bb.y)) max(poly_join(bb.y))] - 10; % go from corners to centers
   dyear = 10;  % every year make a plot that fills up to 10 years back
   year0 = 1955;
   year1 = 2012;
   for year = year1:-1:year0;
   
       %% get data
       [X, Y, Z, Ztime] = grid_orth_getDataFromNetCDFGrids(CATALOG.urlPath,...
       OPT.axis(1),OPT.axis(2),OPT.axis(3),OPT.axis(4),... % [-20000 90000 362500 425000]
            'starttime',datenum(year,1,1),...
         'datathinning',5,...
       'searchinterval',-dyear*365.25);
   
      if all(isnan(Ztime(:)))
          warning(['all nan for period ',num2str(year-dyear),' - ',num2str(year)])
      else
   
       %% plot data
         close all
         AX= subplot_meshgrid(2,1,[.03 .01 .1],.05,[nan .05],[nan]);
         axes(AX(1))
         pcolorcorcen(X,Y,Z)
         hold on
         plot(poly_join(bb.x),poly_join(bb.y),'color',[.5 .5 .5])
         plot(L.x,L.y,'k')
         axis equal
   
         dx = range(poly_join(bb.x))./50;
         dy = range(poly_join(bb.y))./50;
         axis([min(poly_join(bb.x))-dx max(poly_join(bb.x))+dx ...
               min(poly_join(bb.y))-dy max(poly_join(bb.y))+dy])
   
         grid on
         tickmap('xy','texttype','text')
         title([OPT.rijkswaterstaat_dataset_name,'  -  ',datestr(datenum(year-dyear,1,1),'yyyy'),' - ',datestr(datenum(year,1,1),'yyyy')])
         clim([-25 0])
         colorbarwithvtext('z [m]','position',get(AX(2),'position').*[1 4 1 .6])
         delete(AX(2))
   
         print2screensizeoverwrite([fileparts(CATALOG.urlPath{1}),filesep,OPT.rijkswaterstaat_dataset_name,'_',num2str(year-dyear),' - ',num2str(year)])
       
       %% plot timestamps
         close all
         AX= subplot_meshgrid(2,1,[.03 .01 .1],.05,[nan .05],[nan]);
         axes(AX(1))
         
         % (code from grid_orth_plotDataInPolygon
         v = unique(Ztime(find(~isnan(Ztime))));
         nv = length(v);
         
         ZtimeIndex=Ztime;
	     for iv=1:nv
	     mask = (Ztime==v(iv));
	     ZtimeIndex(mask)=iv;
	     end
         pcolorcorcen(X,Y,ZtimeIndex)
         hold on
         plot(poly_join(bb.x),poly_join(bb.y),'color',[.5 .5 .5])
         plot(L.x,L.y,'k')
         axis equal
   
         dx = range(poly_join(bb.x))./50;
         dy = range(poly_join(bb.y))./50;
         axis([min(poly_join(bb.x))-dx max(poly_join(bb.x))+dx ...
               min(poly_join(bb.y))-dy max(poly_join(bb.y))+dy])
   
         grid on
         tickmap('xy','texttype','text')
         title([OPT.rijkswaterstaat_dataset_name,'  -  ',datestr(datenum(year-dyear,1,1),'yyyy'),' - ',datestr(datenum(year,1,1),'yyyy')])

         %% find unqiue date values

         if nv > 0
         caxis   ([1-.5 nv+.5])
         colormap(jet(nv));
         [ax,c1] =  colorbarwithvtext('',1:nv+1,'position',get(AX(2),'position').*[1 4 1 .6]); %#ok<NASGU>
         set(ax,'yticklabel',datestr(v,29))
         end

         delete(AX(2))
   
         print2screensizeoverwrite([fileparts(CATALOG.urlPath{1}),filesep,OPT.rijkswaterstaat_dataset_name,'_',num2str(year-dyear),' - ',num2str(year),'_datestamp'])
      end
   end
   
   toc