cycling/fixpt.m

   1 function [i,j,new_lat,new_lon]= fixpt(w,pt,varargin)
   2 % compute the nearest point in the fire mesh to the given coordinates
   3 % w - w = read_wrfout_tign(f)  , struct with fxlong,fxlat or red =
   4 %                                                     subset_domain(w)
   5 %  pt   - point to be fixed to fire mesh pt = [lat,lon]
   6 %
   7 % i,j  - indices in the grids fxlong,fxlat  or  xlong,xlt
   8 % new_lat,new_lon  -- fxlat(i,j) fxlong(i,j)
   9 %
  10 % jan mandel, june 2013
  11 % modified by james haley june 2020
  12
  13 lon=w.fxlong;
  14 lat=w.fxlat;
  15
  16 if nargin > 2
  17     lon = w.xlong;
  18     lat = w.xlat;
  19 end
  20 min_x=min(lon(:));
  21 max_x=max(lon(:));
  22 min_y=min(lat(:));
  23 max_y=max(lat(:));
  24 [m,n]=size(lon);
  25 [m1,n1]=size(lat);
  26 if(m1 ~= m | n1 ~= n),
  27     error('inconsistent size of FXLONG and FXLAT')
  28 end
  29 %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);
  30
  31
  32 %%%%%%%% compare with what is in w
  33 %find center opf domain
  34 lon_ctr=mean(lon(:));
  35 lat_ctr=mean(lat(:));
  36 %fprintf('the center of the domain is at coordinates %g %g\n',lon_ctr,lat_ctr)
  37 unit_fxlat=6370*2*pi/360;   % one degree latitude in m
  38 unit_fxlon=cos(lat_ctr*2*pi/360)*unit_fxlat; % one degree longitude in m
  39
  40 %fprintf('coordinate units are %g %g m\n',unit_fxlat,unit_fxlon)
  41
  42 % find point on the fire mesh
  43 x=pt(2);
  44 y=pt(1);
  45 if (x<min_x | x>max_x | y<min_y | y>max_y),
  46     error('the point is outside of the domain')
  47 end
  48 % find the nearest point (lon(i,j),lat(i,j)) to (x,y)
  49 d = sqrt((unit_fxlon*(lon-x)).^2 + (unit_fxlat*(lat-y)).^2);
  50 [p,q,minvalue]=find(d==min(d(:)));
  51 for ii=1:length(p)
  52     i=p(ii);
  53     j=q(ii);
  54     new_lon = lon(i,j);
  55     new_lat = lat(i,j);
  56     %fprintf('nearest mesh point %i %i at coordinates %12.10g %12.10g distance %g\n',...
  57     %    i,j,lon(i,j),lat(i,j),d(i,j))
  58 end
  59
  60
  61
  62 end
  63