function [xyh_out]=gttrack(gtfile,xax,yax,xyh_in);
%[out]=gttrack(gt_file,xax,yax,in);
%magic program to do a "gridtrack" deal on a really nasty geotiff file
[nrin,ncin]=size(xyh_in);
xcol=1;
ycol=2;
xyh_out=[xyh_in,(zeros(nrin,1))];
gtcol=ncin+1;
xmax=max(xax);
xmin=min(xax);
ymax=max(yax);
ymin=min(yax);
nx=length(xax);
ny=length(yax);
xstep=mode(abs(xax(1:nx-1)-xax(2:nx)));
ystep=mode(abs(yax(1:ny-1)-yax(2:ny)));
%disp([xstep,ystep]);
xmax2=xmax+xstep/2;
xmin2=xmin-xstep/2;
ymax2=ymax+ystep/2;
ymin2=ymin-ystep/2;
%disp(ystep)
%disp(size(yax))
for n=1:nrin
    x=xyh_in(n,xcol);
    y=xyh_in(n,ycol);
    %disp([x])
    %disp(y)
   % disp(ymax);
   % disp(ymin2);
    if x<=(xmax2) && x>=(xmin2)
            %disp([x,y])
            fdx=find(abs(xax-x)<=(xstep/2));
    else
        if x>xmax2
            fdx=find(xax==xmax);
        else 
            fdx=find(xax==xmin);
        end
    end    
    if y<=(ymax2)  && y>=(ymin2)
        fdy=find(abs(yax-y)<=(ystep/2));
    else
        if y>ymax2
            fdy=find(yax==ymax);
        else 
            fdy=find(yax==ymin);
        end
        %disp('y')
        %fdy0=find(abs(yax-y)<=(ystep/2));
        %fdy=emptysub(fdy0,yax,y);
    end
    xyh_out(n,gtcol)=gtfile(fdy(1),fdx(1));
end
%zerodex=find(xyh_out(gtcol)==0);
%nzdex=length(zerodex);
%if nzdex>=1
%    for n2=1:nzdex
%        n3=zerodex(n2);
        