vis3d/clamp2mesh.m

   1 % function clamp2mesh
   2 % compute the nearest point in the fire mesh to the given coordinates
   3 % use to determine ignition points that lie on the mesh
   4 % run this in the directory with your wrfinput file
   5
   6 % jan mandel, june 2013
   7
   8 d=input('enter domain number: ');
   9 f=sprintf('wrfinput_d%02i',d);
  10 fprintf('reading file %s\n',f);
  11 lon=ncread(f,'FXLONG');
  12 lat=ncread(f,'FXLAT');
  13 min_x=min(lon(:));
  14 max_x=max(lon(:));
  15 min_y=min(lat(:));
  16 max_y=max(lat(:));
  17 [m,n]=size(lon);
  18 [m1,n1]=size(lat);
  19 if(m1 ~= m | n1 ~= n),
  20     error('inconsistent size of FXLONG and FXLAT')
  21 end
  22 fprintf('loaded mesh size %i by %i coordinates %g to %g by %g to %g\n',m,n,min_x,max_x,min_y,max_y);
  23 if(max_x == 0 | max_y == 0)
  24     error('Fire mesh coordinates not set. Please use wrfinput file from a current real.exe or ideal.exe')
  25 end
  26 longlat=input('Enter 1=real run or 0=ideal: ');
  27 if longlat==0, % ideal
  28     fprintf('ideal: ')
  29     unit_fxlon = 1;
  30     unit_fxlat = 1;
  31 else
  32     fprintf('real: ')
  33     lon_ctr=mean(lon(:));
  34     lat_ctr=mean(lat(:));
  35     fprintf('the center of the domain is at coordinates %g %g\n',lon_ctr,lat_ctr)
  36     unit_fxlat=6370*2*pi/360;   % one degree latitude in m
  37     unit_fxlon=cos(lat_ctr*2*pi/360)*unit_fxlat; % one degree longitude in m
  38 end
  39 fprintf('coordinate units are %g %g m\n',unit_fxlat,unit_fxlon)
  40
  41 while 1
  42     x=input('enter the 1st coordinate of ignition point (x or longitude), or enter to exit: ');
  43     if isempty(x),
  44         return
  45     end
  46     y=input('enter the 2nd coordinate of ignition point (y or latitude), or enter to exit: ');
  47     if isempty(y),
  48         return
  49     end
  50     fprintf('entered coordinates %g %g\n',x,y)
  51     if (x<min_x | x>max_x | y<min_y | y>max_y),
  52         error('the point is outside of the domain')
  53     end
  54     % find the nearest point (lon(i,j),lat(i,j)) to (x,y)
  55     d = sqrt((unit_fxlon*(lon-x)).^2 + (unit_fxlat*(lat-y)).^2);
  56     [p,q,minvalue]=find(d==min(d(:)));
  57     for ii=1:length(p)
  58         i=p(ii);
  59         j=q(ii);
  60         fprintf('nearest mesh point %i %i at coordinates %12.10g %12.10g distance %g\n',...
  61             i,j,lon(i,j),lat(i,j),d(i,j))
  62     end
  63
  64 end
  65 % end
  66