% AVHRR_png.m
% Trace les images AVHRR du catalogue:
% \\Pcrougier\GOLTS\Osis\CMS_GOLTS\png_gray
% Permet de choisir l'image à tracer.
% Fabrique une palette de couleur centrée sur zmin, zmax :
%   -zmin : valeur min des pixels de l'image # 0.
%   -zmax : valeur max des pixels de l'image # 255.
% Si cette palette ne convient pas, possibilité de réajuster zmin et zmax /
% aux pixels ayant la plus forte occurence : méthode histogramme.
% NB: la vraie palette (avec nuages en blanc et terre en gris) n'est pas
% affichée en entier (cf pg AVHRR_test.m).
% Plot les isobathes 100, 200, 500, 1000, 1500, 2000 m

%adaptation pour etudiants ayant telecharge eux-memes les images
% format png a partir de osis (AP 09)

close all;
clear all;
ici=pwd;

%GOL
ww=500;  % width image
hh=260      % height image
ll= 850;    % left
tt=140;     % top


% désarchivage de l'image satellite sélectionnée
% if isunix
%     addpath '/homes/grougier/golts/matlab';
%     path_png='/export/home/GOLTS/Osis/CMS_GOLTS/png_gray';
%     [f_png,p_png]=uigetfile([path_png,'/*.png']);
% else
%     %path_png='\\Pcrougier\GOLTS\Osis\CMS_GOLTS\png_gray';
%     path_png='g:\golts\AVHRR\png_gray';
    [f_png,p_png]=uigetfile(['*.png']);
% end

%(rappel: user_entry =  chemin du r€pertoire
%             d = dir(*.extension) (contient donc les noms des fichiers)
%             suivi du r€pertoire où effectuer la d€compression (sans 
%oublier de laisser un espace))

% extraction de la date de l image
[p_raw,f_raw,ext]=fileparts([p_png,f_png]);
[pxx,name,ext,version]=fileparts(f_raw);
titre_image=name;
tan=titre_image(1:4);
tmois=titre_image(5:6);
tjour=titre_image(7:8); 
theure=titre_image(9:10);
tminute=titre_image(11:12);

%lecture pour le canal 4
imc4=imread([p_png,f_png]);
imc4=double(imc4);

% TRACET DE L'IMAGE figure 1
figure(1);
hold on;
axis([2.5 7.5 41 43.6]);
%asp=[0.3 0.3 1];
asp=[1 cos(43*pi/180) 1];
set(gca,'DataAspectRatioMode','manual','DataAspectRatio',asp);

%positionnement de l'image en lat,lon
xim=(-6+(ll-1)*0.01)+[1:ww]*0.01;
%xim=2.49+[1:ww]*0.01;
yim=(45-tt*0.01-hh*0.01)+[hh:-1:1]*0.01;
%yim=41.+[hh:-1:1]*0.01;
him=image(xim,yim,imc4);
[xc,yc]=meshgrid(xim,yim);
[C,h]=contour(xc,yc,imc4,[50 50],'k'); % 85~13.5°
% choisir la couleur de l'isocontour
set(h,'color',[0.1 0.1 0.6]);

% colormap centrée sur les min et max pour les images infra rouge
ind=find(imc4~=0 & imc4~=255);
zmax=max(imc4(ind));
zmin=min(imc4(ind));
jet_ir=colormap(jet(zmax-zmin+1));
%construction de la barre de couleur
map_fin(1,:)=[1 1 1]; % blanc pour les nuages
for i=2:zmin
    map_fin(i,:)=jet_ir(1,:);
end
map_fin(zmin+1:zmax+1,:)=jet_ir;
for i=zmax+2:254
    map_fin(i,:)=jet_ir(end,:);
end
map_fin(255,:)=[0.5 0.5 0.5]; % gris pour la terre
c=colormap(map_fin);
co=colorbar;
% pour avoir une légende de barre de couleur en température (et non pas en count) strictement monotone
set(co,'ylim',[zmin,zmax],'ytick',[zmin:10:zmax],'yticklabel',[(zmin/10+5):1:(zmax/10+5)])  

%set(gca,'ydir','normal','ylim',[41.5 43.6],'ytick',[41.5:0.5:43.5],'yticklabel',[41.5:0.5:43.5],...
 %   'xlim',[3 7])
set(gca,'layer','top','fontsize',12)
title(['image AVHRR du ',num2str(tjour),'/',num2str(tmois),'/',num2str(tan),' - ',num2str(theure),'h',num2str(tminute),''])

%----fin figure 1------

choixg=menu('voulez-vous modifier la colorbar ?','oui','non');
if choixg==1
    choix=menu('Que voulez-vous faire ?','ajuster avec l''histogramme ','fixer les limites');
    
    % HISTOGRAMME figure 2
    if choix==1
        figure(2)
        range=255;
        [m,n]=size(imc4);
        [cnts,bins]=histc(reshape(imc4,[1 m*n]),[0:range-1]);
        bar([1:range],cnts,'histo');  
        set(gca,'xlim',[0 255],'xtick',[0 50 100 120 150 200 255])
        title('Histogramme');
        
        % Pointage d'un minimum et d'un maximum
        disp(' Pointez le minimum')
        [zmin,rien]=ginput(1);
        zmin=ceil(zmin);
        disp(' Pointez le maximum')
        [zmax,rien]=ginput(1);
        zmax=ceil(zmax);
        close(2)
    else % fixe les limites par commande au clavier
        Tmin=input('Quelle Tmin>=5 ? valeur entiere : ');
        Tmax=input('Quelle Tmax ? valeur entiere :');
        
%AJOUT du ROUND on 28/08/08 car hic sur k (colorbar monotone)
% les limitres rentrées sont des temperatures à transformer en count
        zmin=round((Tmin-5)*10);
        zmax=round((Tmax-5)*10);
    end
    
    display([zmin zmax])
        
    % RETRACET DE L'IMAGE figure 3
    figure(3)
    hold on
    axis([2.5 7.5 41 43.6]);
    %axis([3 7 41.5 43.6]);
    %axis([4.8 6.1 42.2 43.45]);
    %asp=[0.3 0.3 1];
    asp=[1 cos(43*pi/180) 1];
    set(gca,'DataAspectRatioMode','manual','DataAspectRatio',asp);
    
    xim=2.49+[1:500]*0.01;
    yim=41.+[260:-1:1]*0.01;
    him=image(xim,yim,imc4);
    [xc,yc]=meshgrid(xim,yim);
    %[C,h]=contour(xc,yc,imc4,[80 80],'k'); % 85~13.5°
    %% choisir la couleur de l'isocontour
    %set(h,'color',[0.1 0.1 0.6])
    
    % colormap pour les images infra rouge
    jet_ir=colormap(jet(zmax-zmin+1));
    %construction de la barre de couleur
    clear map_fin
    map_fin(1,:)=[1 1 1];
    for i=2:zmin
        map_fin(i,:)=jet_ir(1,:);
    end
    map_fin(zmin+1:zmax+1,:)=jet_ir;
    for i=zmax+2:254
        map_fin(i,:)=jet_ir(end,:);
    end
    map_fin(255,:)=[0.5 0.5 0.5];
    c=colormap(map_fin);
    
    co=colorbar;
    %% Position de la colorbar
    Pc=get(co,'Position');  
    
% calcul des temperatures
    Tmin=zmin/10+5;
    Tmax=zmax/10+5;
    
    % pour avoir une légende de barre de couleur strictement monotone
    set(co,'ylim',[zmin,zmax],'ytick',[zmin:10:zmax],'yticklabel',[Tmin:1:Tmax])

    box on
%    set(gca,'ydir','normal','ylim',[41.5 43.6],'ytick',[41.5:0.5:43.5],'yticklabel',[41.5:0.5:43.5],...
%        'xlim',[3 7])
    set(gca,'layer','top','fontsize',12)
   title(['image AVHRR du ',num2str(tjour),'/',num2str(tmois),'/',num2str(tan),' - ',num2str(theure),'h',num2str(tminute),''])

   close(1); %fermeture de l'image non contrastee
end   



% ajout du point SOFI
%Sofi=[5.13 43.07];
%plot(Sofi(1,1),Sofi(1,2),'xr');

% % AJOUT DE LA BATHYMETRIE
%     %isobathes à partir de la bathymétrie réaliste
%     cd d:\bathy;
    load line100.txt;
    load line200.txt;
    load line500.txt;
    load line1000.txt;
    load line1500.txt;
    load line2000.txt;
    %load line2500.txt;
    
    h1=plot(line100(:,2),line100(:,1),'k');
    set(h1,'color',[0.6 0.6 0.6]);
    h2=plot(line200(:,2),line200(:,1),'k');
    %set(h2,'color',[0.2 0.2 0.2],'linewidth',2)  % ligne des 200 m en noir
    set(h2,'color',[0.6 0.6 0.6]);
    h3=plot(line500(:,2),line500(:,1),'k');
    set(h3,'color',[0.6 0.6 0.6]);
    h4=plot(line1000(:,2),line1000(:,1),'k');
    set(h4,'color',[0.6 0.6 0.6]);
    h5=plot(line1500(:,2),line1500(:,1),'k');
    set(h5,'color',[0.6 0.6 0.6]);
    h6=plot(line2000(:,2),line2000(:,1),'k');
    set(h6,'color',[0.6 0.6 0.6]);

xlabel('Longitude','fontsize',12);
ylabel('Latitude','fontsize',12);
cd (ici);

