gridcbass2

%
 cellsize in degrees
 gridsize determines output grid  -gridsize:cellsize:gridsize

0001 function [XI YI ZI hits] = gridcbass2(x,y,data,cellsize,gridsize,sigma,convsize)
0002 %%
0003 % cellsize in degrees
0004 % gridsize determines output grid  -gridsize:cellsize:gridsize
0005 
0006 %Setup a grid
0007 
0008 ti=-1*gridsize:cellsize:gridsize;
0009 gridpoints = 1+(2*gridsize);
0010 
0011 [XI,YI]=meshgrid(ti,ti);
0012 
0013 ZI = zeros(size(XI));
0014 hits = zeros(size(XI));
0015 ZIsize = size(ZI,1);
0016 %%
0017 %Now bin up the data
0018 % First find all the data that is within the grid
0019 ingrid =  (abs(x)<=gridsize) & (abs(y)<=gridsize);
0020 
0021 xbin = (floor(x(ingrid)./cellsize) + ceil(size(ti,2)/2));
0022 ybin = (floor(y(ingrid)./cellsize) + ceil(size(ti,2)/2));
0023 datain = data(ingrid);
0024 
0025 x_cell = (x(ingrid)./cellsize)+ ceil(size(ti,2)/2);
0026 y_cell = (y(ingrid)./cellsize)+ ceil(size(ti,2)/2);
0027 
0028 for k = 1:length(datain);
0029     i=xbin(k);
0030     j=ybin(k);
0031 % for p=-convsize:1:convsize; % units of cells ie is an integer e.g. convisze=1 gives a 3x3 grid
0032 %     for q = -convsize:1:convsize;
0033 %       dist = sqrt((x_cell(k) - (i+p))^2 + (y_cell(k) - (j+q))^2 ); % in cells
0034 %        dist = dist*cellsize;
0035 %        if ((i+p>0) & (i+p<=ZIsize)  & (j+q>0) & (j+q<=ZIsize))
0036           ZI(i,j) = ZI(i,j) + datain(k);
0037           hits(i,j) = hits(i,j) + 1;
0038            %ZI(i+p,j+q) = ZI(i+p,j+q) + datain(k)*interp_func(dist,sigma);
0039         
0040 %      end
0041 % end
0042 end
0043 ZI = ZI./hits;
0044 blah = hits>0;
0045 ZI(~blah)=0;
0046 surface(XI,YI,ZI)
0047 end
0048 %%
0049 function f = interp_func(dist,sigma)
0050 %dist is in real units (degrees)
0051 %sigma - default for raster scans shoudl be 0.2degrees (spacing of rasters
0052 %in elevation)
0053 
0054 f = exp(-dist^2/(2*sigma^2));
0055 end       
0056