source: trunk/src/precipPlot.m @ 327

%+
%
% .. _precipPlot.m:
%
% ============
% precipPlot.m
% ============
%
% DESCRIPTION
% ===========
%
% plot precipitation data according to files originaly transmitted to
% Soukèyé Cisse in :file:`Données climatiques_Soukèye.zip`
%
% This zip file contains Excel files. Some have been converted in CSV.
%
% +todo+ :ref:`pre_precip.sh`
%
% This file in mainly a serie of exercices for Soukèyé Cisse to learn octave.
%
% EXAMPLES
% ========
%
% ::
%
%   $ octave
%   octave> more off
%   octave> varamma_startup
%   octave> precipPlot
%
% TODO
% ====
%
% - trouver une solution pour visualiser sur zeus et/ou sur le poste de travail lodyn301
%   les fichiers dans $VARAMMA_OD
%
%        * sur zeus pas de xv ni de gv !!
%        * sur lodyn301 le /usr/temp/soclod différent de /usr/temp/soclod de zeus !!
%
% - embrouille dans les fichiers/données de moyennes mensuelles : pas de figure pour l'année 1995
%   alors qu'il y a des données
%
% - calcul et tracé des cumuls par decades
% - avec subplot mettre toutes figures sur un meme graphe
%
% - last day of mont is not always 31
%   tip eomday function
%
% EVOLUTIONS
% ==========
%
% $Id$
%
% $URL$
%
% - fplod 20100720T165423Z aedon.locean-ipsl.upmc.fr (Darwin)
%
%   * sauvegarde des données des moyennes mensuelles dans
%     $VARAMMA_OD/precipPlot_station_meanmonthly.txt
%     sous forme ASCII (header + 3 colonnes de données)
%   * make it ok running matlab : replace " by '
%
% - soclod+fplod 20100720T16000Z
%
%   * sauvegarder des figures des moyennes mensuelles dans
%     $VARAMMA_OD/precipPlot_station_meanmonthly_yyyy.ps
%
%     implique l'utilisation de varamma_startup
%
%  * ajout de l'exemple dans l'en tête
%
% - soclod 20100720T150000Z
%
%   * changer l'axe des abscisses avec startdate et enddate dans le tracé des moyenne mensuelles
%
% - soclod 20100719TT150000Z
%
%   * calcul de moyennes mensuelles et tracé des moyennes mensuelles des mois de mousson
%
% - soclod 20100716TT140000Z
%
%   * change ian-imonth loops to line loop because of implicit missing values
%
% - soclod 20100716TT130000Z
%
%   * année sur 2 digits sur le graphique precip_journaliere toutes les
%     années
%
% - fplod 20100716T080507Z aedon.locean-ipsl.upmc.fr (Darwin)
%
%   * tick one day on 10 on figure with several month
%
% - fplod 20100716T073903Z aedon.locean-ipsl.upmc.fr (Darwin)
%
%   * bug fix : replace xdecad by xdata, add figure(ifigure) to move
%     one plot
%   * make it ok running matlab : replace " by '
%
% - soclod 20100715
%
%   * dessin precip_journaliere avec label x yyyymmdd
%   * desin precip_journaliere avec subplot
%
% - fplod 20100715T075514Z aedon.locean-ipsl.upmc.fr (Darwin)
%
%   * add a figure using datetime label
%
% - 20100708 soclod+fplod
%
%   * ajout precip journaliere
%
% - fplod 20100709T072024Z aedon.locean-ipsl.upmc.fr (Darwin)
%   * add loop on figure
%   * make it work using matlab 7.4
%
% .. note::
%
%    sous octave load('petitecote.xls.csv') charge le tableau petitecote_xls
%    sous matlab load('petitecote.xls.csv') charge le tableau petitecote
%
% precipitation journaliere de juin a octobre
% a partir du fichier Barkedji_pluviom_SP-2008.xls
%
% pour faire la conversion Ascii de Barkedji_pluviom_SP-2008.xls en :
%  - ne pas avoir de , entre les cellules
%  - ne pas avoir de ligne vide
%  - ne pas avoir de ligne qui coimmencent par de "
%
%
% ::
%
%   $ xls2csv -c" " -q1 Barkedji_pluviom_SP-2008.xls | \
%   sed -e "/^ *$/d" -e "/^\"/d"  -e "/^jour/d" | tr -d '\f' \
%   > $PROJECT/src/Barkedji_pluviom_SP-2008.xls.csv
%
%  pour faire la conversion Ascii de petite\ côte.xls :
%  ++ pas fini
%  xls2csv -c" " -q1 petite\ côte.xls| sed -e "s/^ 19/\" \" 19/" -e "s/^  janvier/\"station\" \"année\" janvier/" > /usr/home/fplod/incas/varamma/varamma_ws/src/petitecote.xls.csv
%
%-
%
ifigure=1;
%
%  Load data from PDXprecip.dat and plot it with symbols

load PDXprecip.dat;         %  read data into PDXprecip matrix
month = PDXprecip(:,1);     %  copy first column of PDXprecip into month
precip = PDXprecip(:,2);    %  and second column into precip

plot(month,precip,'o');     %  plot precip vs. month with circles

xlabel('month of the year');              %  add axis labels and plot title
ylabel('mean precipitation (mm)');
title('Mean monthly precipitation at Portland International Airport');

% initialize tableau petitecote
petitecote=load('petitecote.xls.csv');
% initialize le tableau month
month=1:12;
% start  boucle of  year
for ii=1:2;
   % prepare affichage figure after (not delete the figure before)
   ifigure=ifigure + 1;
   figure(ifigure)
   % initialize tableau content precipitation during the courant month
   precip=petitecote(ii,2:13);
   % initialize year
   iyear=petitecote(ii,1);
   % trace courant figure
   plot(month,precip,'o');
   % define title with courant month and courant year
   xlabel(['month of the year ' num2str(iyear)]);
   % define au title with precipitation in (mm)
   ylabel('precipitation (mm)');
   % give title of figure
   title('monthly precipitation from file petitecote.xls.csv');
% fin de la boucle on year
end
clear ii;
clear month;
clear precip;
clear petitecote

% same plot than above using datetime label
petitecote=load('petitecote.xls.csv');
month=1:12;
for ii=1:2
   ifigure=ifigure + 1;
   iyear=petitecote(ii,1);
   startdate=datenum(iyear,1,1);
   enddate=datenum(iyear,12,31);
   xdata=linspace(startdate,enddate,size(month,2));
   precip=petitecote(ii,2:13);
   plot(xdata,precip,'o');
   set(gca,'XTick',xdata)
   % nb : because of lack of datetick in octave 3.0.2 on zeus
   % we use xticklabel instead of datetick
   xticklabel=datestr(xdata,'yyyy/mm');
   set(gca,'XTickLabel',xticklabel)
   xlabel(['month of the year ' num2str(iyear)]);
   ylabel('precipitation (mm)');
   title('monthly precipitation from file petitecote.xls.csv with datetime');
end
clear ii;
clear month;
clear precip;

load('precip_journaliere')
% plot 1 figure / month
day=1:31;
for imonth=6:10
   ifigure=ifigure + 1;
   figure(ifigure);
   precip=precip_journaliere(:,imonth-4);
   plot(day, precip, 'o');
   xlabel(['day of the month ' num2str(imonth)]);
   ylabel('precipjournaliere (mm)');
   title('precipitation journaliere');
end
load('precip_journaliere')
day=1:31;
year=2008;
for imonth=6:10
    ifigure=ifigure+1;
    figure(ifigure);
    startdate=datenum(year,imonth,1);
    enddate=datenum(year,imonth,31);
    xdata=linspace(startdate,enddate,size(day,2));
    precip=precip_journaliere(:,imonth-4);
    plot(xdata,precip,'o');
    set(gca,'xtick',xdata);
    xticklabel=datestr(xdata,'dd');
    set(gca,'xticklabel',xticklabel);
    xlabel(['day of the month ' num2str(imonth) ' ' num2str(year)]);
    ylabel('precipjournaliere (mm)');
    title('precipitation journaliere');
end

load('precip_journaliere')
day=1:31;
year=2008;
ifigure=ifigure+1;
figure(ifigure);
nl=2;
nc=3;
for imonth=6:10
    % compute the day number of the 1st date on current month plot
    startdate=datenum(year,imonth,1);
    % compute the day number of the last date on current month plot
    % +todo+ last day of mont is not always 31
    enddate=datenum(year,imonth,31);
    % set x values where we want a tick : one tick every ten days
    xdata=startdate:10:enddate;
    % set text associated to each x values where we want a tick: format dd
    xticklabel=datestr(xdata,'dd');
    % fill precip with current month daily precipitation
    precip=precip_journaliere(:,imonth-4);
    % set x to day numbers
    x=day+startdate-1;
    % define position of plot on the figure
    subplot(nl,nc,imonth-5);
    plot(x,precip,'o');
    % control number of ticks and associated text
    set(gca,'xtick',xdata);
    set(gca,'xticklabel',xticklabel);
    xlabel(['day of the month ' num2str(imonth) ' ' num2str(year)]);
    ylabel('precipjournaliere (mm)');
    title('precipitation journaliere');
end

ifigure=ifigure+1;
figure(ifigure)
precip_journaliere=load('velingra_q_1961-2000.xsl.csv');
% initialize le tableau of the years
an=1961:2000;
% initialize the number d'années
nban=size(an,2);
% initialize tableau x sachant que x contienne des numéros de jour
% ayant comme 1ére valeur le 1er jan 1961 et
% comme derniére valeur le 31 dec 2000
% initialize array month with value of month in one year
month=1:12;
x=1+datenum(1961,1,1):nban*size(month,2)*31+datenum(1961,1,1);
% xdata is an array of nban day numbers yyyy/01/01 where we want ticks
xdata=datenum(an,1,1);
% xticklabel is an array of nban strings "yyyy/mm/dd" we want as associated
% to ticks
xticklabel=datestr(xdata,'yy');
% initialize le tableau content precipitation de tous jours, tous les
% mois et de toutes les années
% ++ NaN au lieu de zero
precip=zeros(1,size(an,2)*size(month,2)*31);
for iligne=1:size(precip_journaliere,1)
    % get year in 1st column
    ian=precip_journaliere(iligne,1);
    % get month in 2d column
    imonth=precip_journaliere(iligne,2);
    % d = position of 1st current month and current year in precip array
    d=datenum(ian,imonth,1)-datenum(1961,1,1)+1;
    f=d+31-1;
    precip(d:f)=precip_journaliere(iligne,3:33);
end

% tracer the courant figure
plot(x,precip,'-');
% control number of ticks and associated text
set(gca,'xtick',xdata);
set(gca,'xticklabel',xticklabel);
% definir au title with number of day for courant month and courant year
xlabel('1st day of the month each year');
% definir au title with precipitation in (mm)
ylabel('precipjournaliere (mm)');
% give title of figure
title('annual velingra precipitation');
%
ifigure=ifigure + 1;
figure(ifigure);
% one figure from June 1st to October 30th
load precip_journaliere
month=6:10;
x=1:size(month,2)*31;
plot(x(1:31),precip_journaliere(:,2),'o', x(31+1:31+31), precip_journaliere(:,3),'-',x(31+31+1:31*3), precip_journaliere(:,4),'+',x(31*3+1:31*4), precip_journaliere(:,5),'x',x(31*4+1:31*5),precip_journaliere(:,6),'*');
xlabel('day number (from June 1st to October 30th)');
ylabel('precipjournaliere (mm)');
title('precipitation journaliere');


ifigure=1;
%nl=3;
%nc=3;
station='velingra';
precip_journaliere=load('velingra_q_1961-2000.xsl.csv');
% initialise les annee
domaine_an=1961:2000;
nban=size(domaine_an,2);
% initialise le tableau des mois pour chaque annee
domaine_month=1:12;
nbmonth=size(domaine_month,2);
% initialise le tableau contenant les sommes et les moyennes
% mensuelles de toutes les annee
summensuelle=zeros(nban,nbmonth);
meanmensuelle=zeros(nban,nbmonth);
% boucle sur les annee
for ian=1:nban;
   % boucle sur les mois
   for imonth=domaine_month;
      % calcul de la ligne contenant les precipitations du mois
      % courant/annee courante
      iligne=imonth+((ian-1)*12);
      if (iligne <= size(precip_journaliere,1));
         % calcul de la somme mensuelle des precipitations
         summensuelle(ian,imonth)=sum(precip_journaliere(iligne,3:eomday(ian+1961-1,imonth)+2));
         % calcul de la moyenne mensuelle des precipitations
         meanmensuelle(ian,imonth)=summensuelle(ian,imonth)/eomday(ian+1961-1,imonth);
      else
         summensuelle(ian,imonth)=NaN;
         meanmensuelle(ian,imonth)=NaN;
      end
   end
end
%
% trace des moyennes mensuelles
domaine_an=1961:2000;
nban=size(domaine_an,2);
domaine_month=6:10;
xdata=zeros(1,size(domaine_month,2));
% boucle sur les années à dessiner
for ian=1:nban;
    startdate=datenum(ian+(1961-1),6,1);
    enddate=datenum(ian+1961-1,10,31);
    % set x values were we want tick:one tick for month
    for imonth=domaine_month
        xdata(imonth-5)=datenum(ian+1961-1,imonth,1);
    end
    xticklabel=datestr(xdata,'yyyy/mm/dd');
    % ++todo+++ subplot(nl,nc,(ian+1)-1961);
    % dessine les precipitations moyennes mensuelles de chaque mois de
    % l'annee courante
    if (find(meanmensuelle(ian,domaine_month) > 0 ))
      ifigure=ifigure+1;
      figure(ifigure);
      plot(xdata, meanmensuelle(ian,domaine_month),'o');
      % on impose les bornes x et y des graphiques pour avoir la même echelle
      % sur toutes les figures
      set(gca,'xtick',xdata);
      set(gca,'xticklabel',xticklabel);
      xlabel(['month of the year ' num2str(ian+1961-1)]);
      ylabel('mean monthly of precipitation (mm)');
      title('mean  precipitation velingra');
      % sauvegarde de la figure de l'année courante
      % dans un fichier POSTSCRIPT ${PROJECT_OD}/precipPlot_velingra_meanmonthly_yyyy.ps
      filename=[PROJECT_OD mfilename '_' station '_meanmonthly_' num2str(ian+1961-1) '.ps'];
      print(filename,'-dps');
    else
     disp(['no data in year ' num2str(ian+1961-1)])
    end
end
close all
% save monthly means in an ASCII file
fullfilename=[PROJECT_OD mfilename '_'  station '_meanmonthly.txt'];
disp(['iii : creation de ', fullfilename]);
fid=fopen(fullfilename,'w');
% write header
fprintf(fid,'# %s\n',fullfilename);
fprintf(fid,'# year month mean\n');
for iyear=1961:2000
    for imonth=6:10
       fprintf(fid,'%4.4d %2.2d %5.1f\n', iyear, imonth, meanmensuelle(iyear-1961+1,imonth-5));
    end
end
fclose(fid);
clear fullfilename;