gridvecs.m 1.8 KB
Newer Older
Pierre Cazenave's avatar
Pierre Cazenave committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76
function gridvecs(gridfile,x1,x2,y1,y2,ds,thresh,fname)

%clear all; close all;
%gridfile = '~gcowles/Data/NOAA_MITSG_TKE/mtm1b.nc';
%x1 = 836100;
%x2 = 935370;
%y1 = -186230;
%y2 = -125810;
%ds = 2000;
%fname = 'mtm1b_2km_nansound_vecs.dat';

% Find elements closest to a grid of points for cleaning up vector plots 
%
% function [] = gridvecs(fname,x2,y1,y2,ds)
%
% DESCRIPTION:
% Find elements closest to a grid of points for cleaning up vector plots 
%
% INPUT:
%    gridfile:  FVCOM netcdf output file containing x,y   
%    fname:  list of node numbers 
%    x1:  left side of bounding box
%    x2:  right side of bounding box
%    y1:  bottom of bounding box
%    y2:  top of bounding box
%    ds:  grid lengthscale
%    thresh: threshold distance in units of x.  If nearest FVCOM point to the 
%            vector grid is greater than this distance, no point is assigned.
%    
% OUTPUT:
%    fname: file containing indices of node numbers where vectors should be plot
%
% Author(s):  
%    Geoff Cowles (University of Massachusetts Dartmouth)
%
% References
%==============================================================================


% open grid file and read data 
if(~exist(gridfile))
  error('grid file does not exist')
end;

nc = netcdf(gridfile);
x  = nc{'x'}(:);
y  = nc{'y'}(:);
nVerts = numel(x);
fprintf('number of nodes in the mesh %d\n',nVerts);
nc = close(nc);

% create the grid
[X,Y] = meshgrid(x1:ds:x2,y1:ds:y2);
[il,jl] = size(X);
node = zeros(il*jl,1);
cnt   = 0;
for i=1:il
for j=1:jl
   xloc = X(i,j); 
   yloc = Y(i,j);
   dist = sqrt(   (x-xloc).^2 + (y-yloc).^2);  
   [mind,imin] = min(dist);
   if(mind < thresh);  
     cnt = cnt + 1;
     node(cnt) = imin;
   end;
end;
end;
node = node(1:cnt);  

% dump the list
fid = fopen(fname,'w');
for i=1:cnt
  fprintf(fid,'%d\n',node(i));
end;
fclose(fid);