Maskmod click.m

De WAVEWATCH-III Wiki

% EXAMPLE SCRIPT FOR MODIFYING MASK FILE FOR MULTI-GRID WW3 % (Note : To run this script you have to run 'create_grid_regional.m' % and 'create_grid_global.m' first because this example modifies the mask % file generated from that script)

% 0. Define paths

 bin_dir = '/export/home/ardhuin/TOOLS/MATLAB/GRIDGEN/bin';             % matlab scripts location
 in_dir = '/export/home/ardhuin/TOOLS/MATLAB/GRIDGEN/examples/data';  % reference data location
 out_dir = '/export/home/ardhuin/TOOLS/MATLAB/GRIDGEN/examples/data';   % output directory (for grid files)
  addpath(bin_dir,'-end')
  

% 1. Read particulars from grid to be modified

 fname = 'norgas_2m';                                             % file name prefix
 fname = 'SEMREV_6s';
 [lon,lat] = read_ww3meta([in_dir,'/',fname,'.meta']);             % read parameter file
 Nx = length(lon);
 Ny = length(lat);
 m = read_mask([in_dir,'/',fname,'.maskorig_ascii'],Nx,Ny);        % read mask file
 m2 =m; %read_mask([in_dir,'/',fname,'.mask'],Nx,Ny);        % read mask file
 bot = read_bot([in_dir,'/',fname,'.depth_ascii'],Nx,Ny);
 [Sx1,Sy1]  = read_obstr([in_dir,'/',fname,'.obstr_lev1'],Nx,Ny);

% 2. Read particulars from base grid (grid with which data has to be exchanged)

% fnameb = 'duck_1m';

% fnameb = 'NOOS_sub6';

% fnameb = 'glob_30m';
 fnameb = 'norgas_2m';
%fnameb = 'atn_10m';                                             % file name prefix
 [lonb,latb] = read_ww3meta([in_dir,'/',fnameb,'.meta']);
 Nxb = length(lonb);
 Nyb = length(latb);
 mb = read_mask([in_dir,'/',fnameb,'.maskorig_ascii'],Nxb,Nyb);
% mb = read_mask([in_dir,'/',fnameb,'.mask'],Nxb,Nyb);
 [Sx,Sy] = read_obstr([in_dir,'/',fnameb,'.obstr_lev1'],Nxb,Nyb);

% 3. Define the polygon that describes the computational area figure(1); clf; pcolor(lon,lat,m); shading flat; colorbar; title('Original Mask for grid','fontsize',14); set(gca,'fontsize',14);

figure(12); clf; pcolor(lonb,latb,mb); shading flat; colorbar; title('Mask for base grid','fontsize',14); set(gca,'fontsize',14);


figure(3);

 clf;
 %pcolor(lon,lat,m2-m_new);
 pcolor(lon,lat,bot./1000);%+(Sx1+Sy1).*1);
 shading flat;
 colorbar;
 caxis([-150,10]);
 title(['Final Land-Sea Mask '],'fontsize',14);
 set(gca,'fontsize',14);
 hold on;
 contour(lonb,latb,mb);
 
%interactive generation of the polygon
 morepoly ='Y';
 np=0;
 pxall=[];
 pyall=[];
 npall=[];
 npoly=0;
 while (morepoly == 'Y')
    [px,py]=ginput;
    npoly=npoly+1;
    %close the polygon
    np0=size(px,1);
    px(np0)=px(1);
    py(np0)=py(1);
    hold on;
    plot(px,py);
    npall=[npall np];   
    pxall=[pxall' px']';
    pyall=[pyall' py']';
    
    m_tmp = modify_mask(m,lon,lat,px,py,mb,lonb,latb,0);         % save the mask info for every polygon in 
                                                              % a different variable
    if (npoly==1) 
        m_new = m_tmp;
    else
        
       loc = find(m_tmp~=3);                                        % determine the active cells
       m_new(loc) = m_tmp(loc);                                     % update the final mask for only those active
       clear loc;                                                   % cells
    end
 
    morepoly = input('Add another polygon? Y/N [N]: ', 's');
    if isempty(morepoly)
    morepoly = 'N';
    end

end

 %pause
 %px=[2  51 80  80 109  50  2   2]/60.-76.;
 %py=[32 32 61 151 180 180 180 32]/60.+34.5;
 px=pxall;
 py=pyall;
 np=npall;
 save polygon_mask px py np
 
 %save mask_mod_gascogne_141107.mat lon lat px py m
 %load mask_mod_gascogne.mat px py
 hold on; plot(px,py)
 %fid = fopen('alaska_10minzone.ascii','r');
 %[a1,count] = fscanf(fid,'%f');
 %px = a1(1:2:count);
 %py = a1(2:2:count);

% 4. Compute the new mask

 %m_new = modify_mask(m,lon,lat,px,py,mb,lonb,latb);           % create the final mask
 m_memo=m_new;
 I=find(m==0);
 m_new(I)=0;

lat2=repmat(lat,Nx,1)'; lon2=repmat(lon',1,Ny)'; %I=find(lat2 > 51 & lat2 < 53 & lon2 > -3 & lon2 < -1); %%I=find(lat2 > 51.2 & lat2 < 51.6 & lon2 > -2.99 & lon2 < -3.034); %m_new(I)=1-m_new(I);

% 5. You can have multiple polygons if needed

% px1 = [178 184 184 178 178]; % py1 = [70.5 70.5 72 72 70.5];

% m_tmp = modify_mask(m,lon,lat,px1,py1,mb,lonb,latb);         % save the mask info for every polygon in 
                                                              % a different variable
 %loc = find(m_tmp~=3);                                        % determine the active cells
 %m_new(loc) = m_tmp(loc);                                     % update the final mask for only those active
 %clear loc;                                                   % cells

% 6. Write out new mask file

 write_ww3file([out_dir,'/',fname,'.mask'],m_new);

% 7. Vizualization (this step can be commented out if resources are limited)

figure(13); clf; pcolor(lon,lat,m_new); shading flat; colorbar; title('Final Mask for grid','fontsize',14); set(gca,'fontsize',14); %

I=find(m_new==2); nb=size(I); lat2=repmat(lat,Nx,1); lon2=repmat(lon',1,Ny); latI=lat2(I); lonI=lon2(I); zerI=zeros(nb,1); oneI=zerI+1;

% Information for coarse grid: [ lonI latI zerI zerI oneI]

% Information for fine grid: dx=(lon(Nx)-lon(1))/(Nx-1) IX=1+(lonI-lon(1))/dx; dy=(lat(Ny)-lat(1))/(Ny-1) IY=1+(latI-lat(1))/dy;

[IX IY ] % END OF SCRIPT