source: trunk/src/evolution_msg.m @ 287

function [times_f,totpix]=evolution_msg(yyyymmddb, yyyymmdde, lonmin,lonmax,latmin,latmax,tseuil)

%EVOLUTION_MSG evolution de la mousson

%+
%
% ===============
% evolution_msg.m
% ===============
%
% .. function:: evolution_msg(yyyymmddb, yyyymmdde, lonmin,lonmax,latmin,latmax,tseuil)
%
% DESCRIPTION
% ===========
%
%   :param yyyymmddb: beginning date yyyymmdd
%   :type yyyymmddb: double
%   :raise  yyyymmddb: required
%
%   :param yyyymmdde: end date yyyymmdd
%   :type yyyymmdde: double
%   :raise yyyymmdde: required
%
%   :param lonmin: longitude min, W < 0
%   :units lonmin: deg
%   :type lonmin: double
%   :raise lonmin: required
%
%   :param lonmax: longitude max, W < 0
%   :units lonmax: deg
%   :type lonmax: double
%   :raise lonmax: required
%
%   :param latmin: latitude min, N > 0
%   :units latmin: deg
%   :type latmin: double
%   :raise latmin: required
%
%   :param latmax: latitude max, N > 0
%   :units latmax: deg
%   :type latmax: double
%   :raise latmax: required
%
%   :param tseuil: seuil de température de brillance
%   :units tseuil: K
%   :type tseuil: double
%   :raise tseuil: required
%
% read all MSG filee :file:`${PROJECT_ID}/MSG/{yyyy}/{mm}/msg-tb108_{yyyy}-{mm}-{dd}_15min.nc` from **yyyymmddb** include to **yyyymmdde** excluded.  
%
% recherche les pixels dont la valeur de température de brillance < tseuil 
% (pixel froids) pour trouver un systeme convectif
%
% dessine l'evolution temporelle du nb de pixels froids
%
% EXAMPLES
% ========
%
% see demo 1
%
% La figure d'évolution s'affiche (tip : si les labels sont tassés, on peut
% en cliquant sur la figure dilater l'axe des X)
%
% SEE ALSO
% ========
%
% :ref:`guide data MSG <data_msg>`
%
% :ref:`varamma_startup.m`
%
% Previous step : :ref:`traite_msg-prod.sh`
%
% Use : :func:`split_yyyymmdd`
%
% TODO
% ====
%
% status de retour
%
% split reading (see :func:`MSGread`) and calculate and plot
%
% filtrage, lissage
%
% check argument
%
% missing value
%
% ascii output
%
% hard coded timestep max per file = 96
%
% EVOLUTIONS
% ==========
%
% $Id$
%
% - pinsard 2011-05-04T09:18:25Z loholt1.ipsl.polytechnique.fr (Linux)
%
%   * replace VARAMMA by PROJECT
%   * make it run on loholt1
%   * add yyyymmddb and yyyymmdde parameters
%
% - fplod 20110228T150040Z aedon.locean-ipsl.upmc.fr (Darwin)
%
%   * example moved to demo
%   * add missing end
%
% - fplod 20110121T131304Z aedon.locean-ipsl.upmc.fr (Darwin)
%
%   * complement on header
%
% - fplod 20100726T160540Z aedon.locean-ipsl.upmc.fr (Darwin)
%
%   * suppress usage of netcdf.open (might be used on climserv 
%     but side effect on zeus running matlab)
%
% - fplod 20100726T124715Z aedon.locean-ipsl.upmc.fr (Darwin)
%
%   *  files must be in ${PROJECT_ID}/MSG/yyyy/mm/
%
% - fplod 20100608T115822Z aedon.locean-ipsl.upmc.fr (Darwin)
%
%   * no more datetick if running octave 3.0.2 on zeus
%
%   * correction of output name detected by matlab
%      ::
%
%      ??? Output argument "t" (and maybe others) not assigned during call to "/.autofs/home/fplod/incas/varamma/varamma_ws/src/evolution_msg.m (evolution_msg)".
% - fplod 20100604T083004Z aedon.locean-ipsl.upmc.fr (Darwin)
%
%   * add usage of varamma_startup : no more chemin parameter, PROJECT_ID value
%     used instead
%
% - fplod 20100525T142706Z aedon.locean-ipsl.upmc.fr (Darwin)
%
%   * revision examples
%
% - lelod 20100521
%
%   * creation
%-
%
global application;
global application_version;
%
global PROJECT_ID;
global PROJECT_OD;
%
global netcdf_open_available
%
h1=1;
% nb de quarts d'heure/fichier
h2=96;

% ++ parameters
lonmin=-15;
lonmax=-10;
latmin=15;
latmax=20;
tseuil=239;

%
arg_info=whos('yyyymmddb');
if ~strcmp(arg_info.class,'double')
   disp(['Incorrect type of arg yyyymmddb ' arg_info.class]);
   whos yyyymmddb
   error(usage);
end
clear arg_info
%
arg_info=whos('yyyymmdde');
if ~strcmp(arg_info.class,'double')
   disp(['Incorrect type of arg yyyymmdde ' arg_info.class]);
   whos yyyymmddb
   error(usage);
end
clear arg_info
%
[yyyyb, mmb, ddb] = split_yyyymmdd(yyyymmddb);
[yyyye, mme, dde] = split_yyyymmdd(yyyymmdde);
datenum_current=datenum(yyyyb,mmb,ddb);
datenum_end=datenum(yyyye,mme,dde);
addtodate_f='day';
addtodate_q=1;
ficlist=[];
% liste des fichiers sur période choisie - toute la periode si possible +todo+
while datenum_current < datenum_end
    disp(['iii : looking for ', datestr(datenum_current,'yyyymmdd')]);
    yyyy_s=datestr(datenum_current,'yyyy');
    mm_s=datestr(datenum_current,'mm');
    dd_s=datestr(datenum_current,'dd');
    filename=['msg-tb108_' yyyy_s '-' mm_s '-' dd_s '_15min.nc'];
    fullfilename=[PROJECT_ID 'MSG/' yyyy_s '/' mm_s '/' filename];
    if (exist(fullfilename)) 
       ficlist=[ficlist fullfilename];
       disp(['iii : data found in ', fullfilename]);
    else
       disp(['iii : no fullfilename ', fullfilename]);
    end
    % next step
    datenum_current=addtodate(datenum_current,addtodate_q, addtodate_f);
    clear yyyy_s
    clear mm_s
    clear dd_s
    clear filename
    clear fullfilename
end
clear datenum_current
clear datenum_end
%
% test if at least one file
nbfic=size(ficlist,1);
if (nbfic == 0)
   error(['eee : no files found during [' num2str(yyyymmddb,'%6.6d') ',' num2str(yyyymmdde,'%6.6d') '[']);
end
%
totpix=[];
times_f=[];
for ijour=1:nbfic
   fullfilename=ficlist(ijour,:);
   % Ouverture du fichier
   if netcdf_open_available
      ncid = netcdf.open(fullfilename, 'NC_NOWRITE');
   else
      ncid = netcdf(fullfilename, 'nowrite');
   end
   disp(['iii : ouverture de ', fullfilename]);
   clear fullfilename;
   if ijour == 1
      % Si 1er fichier lu,
      %  - lecture des longitudes et latitudes
      %  - détermination des indices correspondants à la boite
      % nb : on procéde ainsi pour limiter le nb de données à lire
      % mais cela suppose bien sûr que tous les fichers ont les mêmes
      % caractéristiques (ce qui n'est pas vérifié ici)

      % Récupération des variables.
      if netcdf_open_available
         lon=netcdf.getVar(ncid,1);
         lat=netcdf.getVar(ncid,2);
      else
         lon=file{'LONGITUDE'}(:);
         lat=file{'LATITUDE'}(:);
      end
      indicelat1=find(abs(lat-latmin) == min(abs(lat-latmin)),1);
      indicelat2=find(abs(lat-latmax) == min(abs(lat-latmax)),1);
      clear lat;
      indicelon1=find(abs(lon-lonmin) == min(abs(lon-lonmin)),1);
      indicelon2=find(abs(lon-lonmax) == min(abs(lon-lonmax)),1);
      clear lon;
   end
   if netcdf_open_available
      data=netcdf.getVar(ncid,0);
      % apply scale_factor and add_offset
      data = double(data) * netcdf.getAtt(ncid,0,'scale_factor') + netcdf.getAtt(ncid,0,'add_offset');
      % permute array to have
      % dimension 1 : temps, dimension 2 : lat,  dimension 3 : lon
      data = permute(data,[3 2 1]);
   else
      % Set missing value to NaN
      data(data==file{'DATA'}.missing_value)=NaN;
      % apply scale_factor and add_offset
      data=file{'DATA'}.scale_factor(:)*data+file{'DATA'}.add_offset(:);
   end
   % reduce lat-lon domain
   data=data(:,indicelat1:indicelat2,indicelon1:indicelon2);
   % lecture du temps
   if netcdf_open_available
    hours_since19600101=netcdf.getVar(ncid,3);
   else
    hours_since19600101=file{'TIME'}(:);
   end
   hfin=size(hours_since19600101,1);
   % close NetCDF file
   if netcdf_open_available
      netcdf.close(ncid);
   else
      ncclose(ncid);
   end
   clear ncid
   %
   % conversion du temps lu (HOURS since 1960-01-01 00:00:00) en date
   times_f=[times_f; hours_since19600101/24.+datenum(1960,01,01)];
   clear hours_since19600101;
   %
   tbmsg=data;
   clear data;
   % n pixels froid dans chaque image de chaque jour
   npix=zeros(hfin,1);
   for iheure=h1:hfin
   %   disp(['iii : recherche des pixels froids pour le pas de temps ',num2str(iheure), ' du jour ', num2str(ijour)]);
      % IDL pixfroid=where(tbmsg(iheure) > tseuil,nb);
      pixfroid=find(tbmsg(iheure,:,:) < tseuil);
      nb=size(pixfroid,1);
      clear pixfroid;
 %     nbtot=size(tbmsg(iheure,:,:),2) * size(tbmsg(iheure,:,:),3);
 %     disp(['iii : ', num2str(nb), ' pixels froids parmi ', num2str(nbtot)]);
%++      if nb ~= 0
         % tableau du jour
         npix(iheure)=nb;
%++      end
   end
   clear tbmsg;
   % concatenation
   totpix=[totpix;npix];
end
clear indicelat1;
clear indicelat2;
clear indicelon1;
clear indicelon2;
clear h1;
clear h2;
clear latmin;
clear latmax;
clear lonmin;
clear lonmax;
%
% Create xdata to correspond to the number of points between the start and
% end dates:
nbp=nbfic*96;
times_str=datestr(times_f,6);
%xdata = linspace(times_f(1),times_f(size(times_f,1)),nbp);
%
% sauvegarde des infos en ascii
%PROJECT_OD='./';
%fullfilename=[PROJECT_OD, mfilename, '_', times_str(1,:), '_', times_str(size(times_str,1),:),'.txt'];
%disp(['iii : creation de ', fullfilename]);
%fid=fopen(fullfilename,'w');
% écriture en tête
%fprintf(fid,'# %s\n',fullfilename);
%fprintf(fid,'# temps nb_froid\n');
%fprintf(fid,'!!!!!! ATTENTION CE FICHIER EST INCORRECT !!!!! \n');
%fprintf(fid,'!!!!!! mise au point en cours !!!!! \n');
%fprintf(fid,'%f %5.5d\n',times_f,totpix);
%++ pas ok pb écriture string fprintf(fid,'%s %5.5d\n',times_str(1),totpix);
%fclose(fid);
%clear fullfilename;
%

% plot de nb points froids au cours du temps

figure(1);
clf;
plot(times_f,totpix)

% Set the number of XTicks to the number of points in xdata
%set(gca,'XTick',xdata)

% conversion des données temps en chaine mmdd
%xlabel=datestr(xdata,6);
switch application
  case {'matlab'}
     datetick('x',6);
  case {'octave'}
     switch application_version
        case {'3.0.2'}
           disp(['no datetick available using ', application , ' ', application_version]);
        otherwise
           error(['unknown application version running ', application , ' ' ,application_version]);
     end
  otherwise
     error('unknown application running');
end
%set (gca, 'xticklabel', xlabel);
%get (gca, 'xticklabelmode');
xlabel('date','fontsize',8)
ylabel('nb pixels froids');
title(['evolution du nb de pixels froids dans les donnees MSG entre ', times_str(1,:), ' et ', times_str(size(times_str,1),:)]);

% sauvegarde de la figure
%printer='ps';
%print_printer=['-d', printer];
%fullfilename=[PROJECT_OD, mfilename, '_', times_str(1,:), '_', times_str(size(times_str,1),:),'.', printer];
%clear printer;
%print(print_printer,fullfilename);
%disp(['iii : creation de ', fullfilename]);
%clear print_printer;
%clear fullfilename;
%

end

%!demo
%! clear all
%! varamma_startup
%! more off
%! yyyymmddb=20060801;
%! yyyymmdde=20060802;
%! lonmin=-15;
%! lonmax=-10;
%! latmin=15;
%! latmax=20;
%! tseuil=239;
%! [times_f,totpix]=evolution_msg(yyyymmddb, yyyymmdde, lonmin,lonmax,latmin,latmax,tseuil);