% This script was made as part of the NMDC project "Visualisatie en user interfaces"
% of the Dutch National Model and Data Centre
% http://nmdc.eu/index.php?option=com_content&view=article&id=49&catid=14
%
% It operates on grib files obtained from KNMI. The same data will be made
% available as netCDF too via their OPeNDAP server
% http://geoservices.knmi.nl/thredds/catalog/knmi/opendap/HIRLAM_DEMO/catalog.html

clear all
close all


   OPT.dir = 'D:\PROJECTS\HIRLAM\data\';
   OPT.kml = 'D:\PROJECTS\HIRLAM\data\kml\';
   
   L=dir([OPT.dir,filesep,'*2009*.grib']);
   
   for i=1:length(L)

      FI1 = grib('open',[OPT.dir,filesep,L(i).name]);
      if i==length(L)
      FI2 = FI2;
      else
      FI2 = grib('open',[OPT.dir,filesep,L(i).name]);
      end

      p{i}    = [OPT.kml,filesep,'pc_'  ,filename(FI1.FileName),'.kml'];
      uabs{i} = [OPT.kml,filesep,'U_'  ,filename(FI1.FileName),'.kml'];
      uvc{i}  = [OPT.kml,filesep,'uvc_',filename(FI1.FileName),'.kml'];
      uv{i}   = [OPT.kml,filesep,'uv_' ,filename(FI1.FileName),'.kml'];

     [P ,LAT,LON] = GRIB('read',FI1,1);
      if i>1
      U1    = U2;
      V1    = V2;
      UABS1 = UABS2;
      else
     [U1,LAT,LON] = GRIB('read',FI1,3);
     [V1,LAT,LON] = GRIB('read',FI1,4);
     [~,UABS1] = cart2pol(U1,V1);          
      end
      if i<length(L)
     [U2,LAT,LON] = GRIB('read',FI2,3);
     [V2,LAT,LON] = GRIB('read',FI2,4);
     [~,UABS2] = cart2pol(U2,V2);
      end
      
%% contour3(p)

      if 1
      j=1;
      [~,pCB{i}]=KMLcontour3(LAT,LON,P,'fileName',p{i},... % UABS,
                                 'levels',95000:100:105000,... % 0:3:18,...%
                           'CBcolorTitle','pressure [hPa]',...
                                 'timeIn',FI1.Block(j).Info.Date,...
                                'timeOut',FI1.Block(j).Info.Date+6/24,...
                                'colorbar',i==1,...
                               'zScaleFun',@(z) (z - 1e5).*1e2,...
                             'writeLabels',0);%,...
                                %'colorMap',@(x)colormapgray(0,.8,x)); %,...% @(x)[.7 .7 .7],...'colorSteps',1);
      end

%% contour3(p,U)

      if 0
      j=1;
      [~,uabsCB{i}]=KMLcontour3(LAT,LON,P,UABS,'fileName',uabs{i},... 
                                 'levels',0:3:18,...
                           'CBcolorTitle','speed [m/s]',...
                                 'timeIn',FI.Block(j).Info.Date,...
                                'timeOut',FI.Block(j).Info.Date+6/24,...
                                'colorbar',i==1,...
                               'zScaleFun',@(z) (z - 1e5).*1e2,...
                             'writeLabels',0,...
                                 'extrude',0,...
                                 'zstride',10,...
                               'lineWidth',2);
                                                      
      end
%% quiver

      if 0
      KMLquiver(LAT,LON,V,U,'arrowScale',2e4,...
                                 'timeIn',FI.Block(3).Info.Date,...
                                'timeOut',FI.Block(3).Info.Date+6/24,...
                             'lineWidth',2,...
                             'fillAlpha',1,...
                              'fileName','v.kml',...
                                'stride',31,...
                              'fileName',uv{i})
      end
%% KMLcurvedArrows

      if 0
      dn=1; 
        [X,Y,log]=convertCoordinates(LON,LAT,'CS1.code',4326,'CS2.code',32631); % WGS 84 / UTM zone 31N
        if i==1;x0 = [];y0=[];;end
        [x0,y0]=KMLcurvedArrows(X(1:dn:end,1:dn:end),...
                        Y(1:dn:end,1:dn:end),...
                       {U1(1:dn:end,1:dn:end),U2(1:dn:end,1:dn:end)},...
                       {V1(1:dn:end,1:dn:end),V2(1:dn:end,1:dn:end)},...
                           'timeIn',FI1.Block(3).Info.Date,...
                          'timeOut',FI1.Block(3).Info.Date+6/24,...
                             'time',[],...
                     'coordConvFun',@(x,y,EPSG) convertCoordinates(x,y,EPSG,'CS1.code',32631,'CS2.code',4326),...
                         'n_arrows',400,...
                         'fileName',uvc{i},...
                     'interp_steps',0,...
                            'nhead',2,...
                           'length',3.*3600,...
                         'colorMap',@(m) [1 1 1],...
                       'colorSteps',1,...
                       'x0',x0,'y0',y0);
      end                          
                              
   end  % time
   
   try;KMLmerge_files('filename','p.kml'   ,'sourceFiles',{p{:}   ,pCB{1}{1}})   ;end
   try;KMLmerge_files('filename','uabs.kml','sourceFiles',{uabs{:},uabsCB{1}{1}});end
   try;KMLmerge_files('filename','uvc.kml' ,'sourceFiles',uvc)                   ;end
    
% 2011-05-26 23:00:00
% 2011-05-27 02:00:00 *
% 2011-05-27 05:00:00
% 2011-05-27 08:00:00
% 2011-05-27 11:00:00
% 2011-05-27 14:00:00
% 2011-05-27 17:00:00
% 2011-05-27 20:00:00
% 2011-05-27 23:00:00
% 2011-05-28 02:00:00
% 2011-05-28 05:00:00
% 2011-05-28 08:00:00
% 2011-05-28 11:00:00
% 2011-05-28 14:00:00
% 2011-05-28 17:00:00
% 2011-05-28 20:00:00

% t = nc_cf_time('http://geoservices.knmi.nl/thredds//dodsC/knmi/opendap/HIRLAM_DEMO/HIRLAM72_1306454400.nc')
%
% 2011-05-27 02:00:00
% 2011-05-27 05:00:00
% 2011-05-27 08:00:00 *
% 2011-05-27 11:00:00
% 2011-05-27 14:00:00
% 2011-05-27 17:00:00
% 2011-05-27 20:00:00
% 2011-05-27 23:00:00
% 2011-05-28 02:00:00
% 2011-05-28 05:00:00
% 2011-05-28 08:00:00
% 2011-05-28 11:00:00
% 2011-05-28 14:00:00
% 2011-05-28 17:00:00
% 2011-05-28 20:00:00
% 2011-05-28 23:00:00
% 2011-05-29 02:00:00

% t = nc_cf_time('http://geoservices.knmi.nl/thredds//dodsC/knmi/opendap/HIRLAM_DEMO/HIRLAM72_1306476000.nc')
%
% 2011-05-27 08:00:00
% 2011-05-27 11:00:00
% 2011-05-27 14:00:00
% 2011-05-27 17:00:00
% 2011-05-27 20:00:00
% 2011-05-27 23:00:00
% 2011-05-28 02:00:00
% 2011-05-28 05:00:00
% 2011-05-28 08:00:00
% 2011-05-28 11:00:00
% 2011-05-28 14:00:00
% 2011-05-28 17:00:00
% 2011-05-28 20:00:00
% 2011-05-28 23:00:00
% 2011-05-29 02:00:00
% 2011-05-29 05:00:00
% 2011-05-29 08:00:00
