WebSVN – lagranto.ecmwf – /trunk/figure/matlab/plot_airparcels.m

% ------------------------------------------------------------------------
% Preparations
% ------------------------------------------------------------------------

% Add some matlab paths
addpath('/usr/local/matlabtools/mexnc/');
addpath('/usr/local/matlabtools/snctools/');
addpath('/usr/local/matlabtools/cdf_io/');
addpath('/home/sprenger/lagranto/matlab/');

% Set the datadirectories
cdfdir   = '/lhome/sprenger/lagranto/test/cdf/';
tradir   = '/home/sprenger/lagranto/matlab/';

% Read a trajectory file
tra = read_trajectory( [ tradir '/trajectory2.lsl' ]);

% Set the tabel for input files
filelist.n = 4;
filelist.date(1) = cellstr('19891020_00'); filelist.time(1) =  0;
filelist.date(2) = cellstr('19891020_06'); filelist.time(2) =  6;
filelist.date(3) = cellstr('19891020_12'); filelist.time(3) = 12;
filelist.date(4) = cellstr('19891020_18'); filelist.time(4) = 18; 

% Get the grid parameters from first data file
inp                     = ncget([ cdfdir 'P' char(filelist.date(1)) ],'U','PS');
grid.nx                 = size(inp.U.data,3);
grid.ny                 = size(inp.U.data,2);
grid.nz                 = size(inp.U.data,1);
grid.xmin               = inp.cstfile.lonmin;
grid.ymin               = inp.cstfile.latmin;
grid.dx                 = inp.cstfile.dellon;
grid.dy                 = inp.cstfile.dellat;
grid.xmax               = grid.xmin + double(grid.nx-1) * grid.dx;
grid.ymax               = grid.ymin + double(grid.ny-1) * grid.dy;
grid.mdv                = -999.
grid.lon1               = grid.xmin + double(0:grid.nx-1) * grid.dx;
grid.lat1               = grid.ymin + double(0:grid.ny-1) * grid.dy;
[ grid.lon2 grid.lat2 ] = meshgrid(grid.lon1,grid.lat1);
grid.aklay              = inp.cstfile.aklay;
grid.aklev              = inp.cstfile.aklev;
grid.bklay              = inp.cstfile.bklay;
grid.bklev              = inp.cstfile.bklev;

% ------------------------------------------------------------------------
% Loop over al trajectory times
% ------------------------------------------------------------------------

% Set flags for load manager
load0 = -1;
load1 = -1;

% No trajectories are selected
select = zeros(tra.ntra,1);

% Init the textlines
for j=1:tra.nfield
    textline(j) = cellstr(' ');
end
first = 1;

% Loop over all trajectory times
i     = 1;
stat  = 1;
while ( stat ~= 0 )

    % Get the time 
    time    = tra.time(i);
    
    % Decide which data files are needed
    need0 = find( time <= filelist.time , 1 );
    need1 = need0 + 1;
    if ( need1 > filelist.n ) 
       need1 = need0;
    end
       
    % Load files if not yet done
    if ( need0 == load1 )
        inp0  = inp1;
        load0 = need0;
    end
    if ( need1 == load0 )
        inp1  = inp0;
        load1 = need1;
    end
    if ( need0 ~= load0 ) 
        inps        = ncget([ cdfdir 'S' char(filelist.date(need0)) ],'TH','PV','RH');
        inpp        = ncget([ cdfdir 'P' char(filelist.date(need0)) ],'Q','U','V','T','PS');
        inp0        = inpp;
        inp0.TH     = inps.TH;
        inp0.PV     = inps.PV;
        inp0.RH     = inps.RH;
        inp0.P      = inp.PS;
        inp0.P.data = kron(inp0.PS.data,grid.bklay) + kron(ones(grid.ny,grid.nx),grid.aklay);
        inp0.P.data = reshape(inp0.P.data,grid.nz,grid.ny,grid.nx);
        load0       = need0;
    end
    if ( need1 ~= load1 ) 
        inps        = ncget([ cdfdir 'S' char(filelist.date(need1)) ],'TH','PV','RH');
        inpp        = ncget([ cdfdir 'P' char(filelist.date(need1)) ],'Q','U','V','T','PS');
        inp1        = inpp;
        inp1.TH     = inps.TH;
        inp1.PV     = inps.PV;
        inp1.RH     = inps.RH;
        inp1.P      = inp.PS;
        inp1.P.data = kron(inp1.PS.data,grid.bklay) + kron(ones(grid.ny,grid.nx),grid.aklay);
        inp1.P.data = reshape(inp1.P.data,grid.nz,grid.ny,grid.nx);
        load1       = need1;
    end

    % Get the mean pressure of the trajectory sample
    if ( any(select) )
        mask = ( find(select) - 1) * tra.ntime + i;
    else
        mask = ( (1:tra.ntra) -1 ) * tra.ntime + i;
    end
    hori.level = mean ( tra.p ( mask ) );
    
    % Interpolate fields to a pressure surface
    [ xpos  ypos zpos ] = meshgrid( grid.lon1, grid.lat1, hori.level );
    ind                 = get_index(inp0.P.data,[ zpos(:) ypos(:) xpos(:) ],grid);
    hori0.U             = reshape( int_index(inp0.U.data,ind,grid.mdv), 1,grid.ny,grid.nx);
    hori0.V             = reshape( int_index(inp0.V.data,ind,grid.mdv), 1,grid.ny,grid.nx);
    hori0.VEL           = sqrt(hori0.U .^ 2 + hori0.V .^ 2);

    [ xpos  ypos zpos ] = meshgrid( grid.lon1, grid.lat1, hori.level );
    ind                 = get_index(inp1.P.data,[ zpos(:) ypos(:) xpos(:) ],grid);
    hori1.U             = reshape( int_index(inp1.U.data,ind,grid.mdv), 1,grid.ny,grid.nx);
    hori1.V             = reshape( int_index(inp1.V.data,ind,grid.mdv), 1,grid.ny,grid.nx);
    hori1.VEL           = sqrt(hori1.U .^ 2 + hori1.V .^ 2);

    % Do a time interpolation (and get rid of 1-dimensions) 
    if ( need1 ~= need0 )
        frac     = ( time - filelist.time(need0) ) / ( filelist.time(need1) - filelist.time(need0) );
        hori.U   = squeeze( (1-frac) * hori0.U   + frac * hori1.U   );
        hori.V   = squeeze( (1-frac) * hori0.V   + frac * hori1.V   );
        hori.VEL = squeeze( (1-frac) * hori0.VEL + frac * hori1.VEL );
    else
        hori.U   = squeeze( hori0.U   );
        hori.V   = squeeze( hori0.V   );
        hori.VEL = squeeze( hori0.VEL );
    end
    
    % Open a new figure ans set the geographical projection
    figure(1);
    clf;
    load coast
    %h=axesm('MapProjection','stereo','origin',[ 90 0 ]);
    h=axesm( 'eqdcylin','MapLatLimit',[ 20 80 ], 'MapLonLimit', [ -100 40 ] );
    gridm;
    h=plotm(lat,long,'Color','k','LineWidth',1.5);
    %axis([-0.6 0.67 -1.2 -0.1]);
    axis([ -1.0    1.0    0.4000    1.4000 ]);
    daspect([ 1 0.75 1] );
    title ( [ 'P = ' num2str(round(hori.level)) ' hPa / T = ' num2str(time) ' h'],...
           'FontSize',22)

    % Plot velocity and overlay wind vectors
    c_map = col_scale([0:2.5:50 ],hori.VEL);
    [C,h] = contourfm(grid.lat1,grid.lon1,c_map.data,c_map.xtick);
    colormap('default');
    col_load( [ pwd '/col_table.txt' ],c_map.xtick,44); 
    caxis(c_map.caxis)
    q=colorbar('vert');
    set(q,'ytick',c_map.xtick,'YTickLabel',c_map.label,'FontSize',15); 
    [ c_map.x c_map.y ] = meshgrid(grid.lon1,grid.lat1);
    c_map.u             = squeeze(hori.U);
    c_map.v             = squeeze(hori.V);
    maskx               = 1:2:size(grid.lon2,2);        
    masky               = 1:2:size(grid.lon2,1);
    hold on
    q=quiverm(grid.lat2(masky,maskx),grid.lon2(masky,maskx),...
              hori.V   (masky,maskx),hori.U   (masky,maskx),'k');
    set(q,'Linewidth',1,'Color','k')
    hold off  
    
    % Handle user commands
    stat = 2;
    while ( stat == 2 )
        
        % Show the positions of the trajectory points
        hold on
        for j=1:tra.ntra
            if select(j) 
                mask = ( tra.time == time ) & ( tra.label == j );
                plotm(tra.lat(mask),tra.lon(mask),'o','markersize',10,'color','r',...
                      'MarkerEdgeColor','k','MarkerFaceColor','r');
            else
                mask = ( tra.time == time ) & ( tra.label == j );
                plotm(tra.lat(mask),tra.lon(mask),'o','markersize',10,'color','k',...
                      'MarkerEdgeColor','k','MarkerFaceColor','k');
            end
        end
        hold off
      
        % Get the cursor position
        [ x y ] = ginput(1)
        plotdom = axis;

       % Step backward
       if ( x < plotdom(1) ) 
           i = i - 1;
           if ( i < 1)
               i = tra.ntime;
           end
           stat = 1;
           
       % Step forward
       elseif ( x > plotdom(2) ) 
           i = i + 1;
           if ( i > tra.ntime )
               i = 1;
           end
           stat = 1;
           
       % Exit
       elseif ( y > plotdom(4) )
           stat = 0;
        
       % Write trajectory info
       elseif ( (y < plotdom(3)) & any(select) )
           for j=1:tra.nfield
              mask = ( find(select) - 1) * tra.ntime + i;
              val  = tra.(char(tra.field(j)));
              val  = mean( val(mask) );
              if ( first == 0 ) 
                 text(plotdom(1),plotdom(3)-0.04*j,char(textline(j)),...
                      'FontSize',15,'Color','g' );
              end
              textline(j) = cellstr([ char(tra.field(j)) ' = ' num2str(val)]);
              text(plotdom(1),plotdom(3)-0.04*j,char(textline(j)),...
                   'FontSize',15,'Color','b' );
              first = 0;
           end
        
       % Select trajectory points
       else
           [ lat0, lon0 ] = minvtran(x,y);
           mask           = (tra.time == time );
           lon1           = tra.lon(mask);
           lat1           = tra.lat(mask);
           dist           = distance(lat0,lon0,lat1,lon1);
           label          = find( min(dist) == dist );
           if  ( dist(label) < 5 ) 
              if ( select(label) == 0 ) 
                 select(label) = 1;
              elseif ( select(label) == 1 )
                 select(label) = 0;
              end
           else
               stat = 1;
           end
       end
       
    end

% End loop 
end
Subversion Repositories lagranto.ecmwf

(root)/trunk/figure/matlab/plot_airparcels.m – Rev 3
