cycling/fire_gradients.m

   1 function [dx,dy,nan_msk] = fire_gradients(lon,lat,tign,unit)
   2 %computes the gradients of the fire arrival time cone
   3 % inputs : tign - fire arrival time matrix, lon = fxlon, lat = fxlat
   4 % output:  ux,uy   - vector components of gradients in the x and y directions
   5 % unit = 1 ==> return unit vectors
   6 E = wgs84Ellipsoid;
   7 % [n,m] = size(lon);
   8 % hx = distance(lat(1,round(m/2)),lon(1,1),lat(1,round(m/2)),lon(end,end),E)/m;
   9 % hy = distance(lat(1,1),lon(round(n/2),1),lat(end,end),lon(round(n/2),1),E)/n;
  10
  11
  12 [aspect,slope,dy,dx] = gradientm(lat,lon,tign,E);
  13 %set all gradients at top of cone to NaN
  14 t_top = tign > max(tign(:))-0.1;
  15 dx(t_top) = NaN;
  16 dy(t_top) = NaN;
  17
  18 %[dx,dy] = gradient(tign);
  19 %using Sobel edge detector
  20 % [dx,dy] = imgradientxy(tign);
  21 % dx = dx/8;
  22 % dy = dy/8;
  23 %make unit vectors
  24 if unit == 1
  25 x = dx./sqrt(dx.^2+dy.^2);
  26 y = dy./sqrt(dx.^2+dy.^2);
  27 dx = x;
  28 dy = y;
  29 end
  30
  31 % figure,contour(fxlong,fxlat,tign)
  32 % hold on
  33 % quiver(fxlong,fxlat,ux,uy)
  34 % hold off
  35
  36 %[gy,gx] = imgradientxy(tign);
  37 mx = isnan(dx);
  38 my = isnan(dy);
  39 nan_msk = logical(mx.*my);
  40 %plot the vectors
  41 fprintf('there were %d NaN in the gradient \n',sum(nan_msk(:)));
  42 % figure,quiver(lon(~nan_msk),lat(~nan_msk),dx(~nan_msk),dy(~nan_msk))
  43 % title('Gradient')
  44 % figure,scatter(lon(nan_msk),lat(nan_msk));
  45 % title(['Location of ',num2str(sum(nan_msk(:))),' NaN in gradient'])
  46 end % function