femwind/topo.m

   1 % a quick script to prepare terrain height for fewmind
   2 % from elevation in a geotiff
   3
   4 filename='GranitePeak-DEM.tif';
   5 sr_x=10
   6 sr_y=10
   7 disp(['reading ',filename])
   8 t = Tiff(filename,'r');
   9 info = geotiffinfo(filename)
  10 dx = info.PixelScale(1)
  11 dy = info.PixelScale(2)
  12 d = read(t);
  13 ds = size(d);
  14 disp(['grid size ',num2str(size(d)),' min ',num2str(min(d(:))),' max ',num2str(max(d(:)))])
  15 isza = floor(ds(1)/sr_x)
  16 iszf = isza * sr_x
  17 jsza = floor(ds(2)/sr_y)
  18 jszf = jsza * sr_y
  19 zsf = d(1:iszf,1:jszf);
  20 [xf,yf] = ndgrid(1:iszf,1:jszf);
  21 [xa,ya] = ndgrid(([1:isza]-0.5)*sr_x,([1:jsza]-0.5)*sr_y);
  22
  23 ht1 = interp2(1:jszf,1:iszf,zsf,ya(:),xa(:));
  24 if isnan(sum(ht1)),error('NaN in interpolant'),end
  25 ht = reshape(ht1,size(xa));
  26
  27 figure(1), mesh(xf,yf,zsf)
  28 %figure(2),
  29 hold on
  30 mesh(xa,ya,ht)
  31 hold off
  32 drawnow
  33
  34
  35 fzsf='input_zsf';
  36 disp(['writing ',fzsf])
  37 write_array_2d(fzsf,zsf)
  38
  39 fht='input_ht';
  40 disp(['writing ',fht])
  41 write_array_2d(fht,ht)
  42