function [out ohead]=sir_reproject(a,ahead,rhead,opt)
%
% function [out ohead]=sir_reproject(a,ahead,rhead <,option>)
%
% reprojects image a with sir header ahead to the same map projection
% and size as rhead (caution: tends to be slow since the lat/lon of each
% pixel in both map projections have to be computed)
%
% Remapping algorithm summary: 
%  Option 1 (averaging + nearest-neighbor) [default]
%    1. For each input pixel, compute lat/lon of pixel center using the SIR
%       pixel -> lat/lon routines
%    2. For each input pixel lat/lon, determine the output pixel in the output
%       image using SIR lat/lon -> pixel routines.  If multiple input pixels
%       hit the same output pixel, they are averaged. (nodata values not used)
%    3. Scan output image.  For each empty output pixel, compute the
%       corresponding pixel center lat/lon
%    4. Compute the input pixel location from the lat/lon and if sufficently
%       close, copy the value.
%  Option 2 (nearest-neighbor)
%    1. For each pixel in the output image, compute lat/lon of pixel center 
%       using the SIR pixel -> lat/lon routines
%    2. Compute the input pixel location from the lat/lon and copy the value.

% written by DGL at BYU 05 Oct 2022

option=1;
if nargin>3, option=opt; end;

% create output file header
ohead=ahead;
ohead(1:9)=rhead(1:9); % copy projection info
ohead(17)=rhead(17); % copy projection type
ohead(127:128)=rhead(127:128); % copy projection storage info
ohead(169)=rhead(169);
ohead(240)=rhead(240);
ohead(255)=rhead(255);
%anodata=sirheadvalue('anodata',ahead);
anodata=rhead(49);
if anodata==0, anodata=-33;end;

out=zeros([rhead(2),rhead(1)])+nan;
if option==2 % Nearest neighbor pixel copying.  Fill output array with 
             % nearest input array pixel
  % for each each output pixel, find nearest input pixel
  ind1=1:rhead(1)*rhead(2);
  [i,j]=ind2sub([rhead(2),rhead(1)],ind1);
  [lon,lat]=pix2latlon(i,j,rhead);
  [x,y]=latlon2pix(lon,lat,ahead);
  xc=round(x); yc=round(y);
  ind=find(xc>0 & xc<=ahead(1) & yc>0 & yc<=ahead(2));
  ind2=sub2ind([ahead(2),ahead(1)],xc(ind),yc(ind));
  out(ind1(ind))=a(ind2);

else % % averaging + nearest-neighbor (slower)   
  % copy input to output, averaging as necessary while handling nodata pixels
  
  ind1=1:ahead(1)*ahead(2);
  [i,j]=ind2sub([ahead(2),ahead(1)],ind1);
  [lon,lat]=pix2latlon(i,j,ahead);
  [x,y]=latlon2pix(lon,lat,rhead);
  xc=round(x); yc=round(y);
  ind=find(xc>0 & xc<=rhead(1) & yc>0 & yc<=rhead(2));
  ind2=sub2ind([rhead(2),rhead(1)],xc(ind),yc(ind));
  cnt=zeros(size(out));
  for n=1:length(ind)
    if a(ind1(ind(n)))>anodata+0.001 
      if isnan(out(ind2(n)))
	out(ind2(n))=a(ind1(ind(n)));
	cnt(ind2(n))=1;
      else
	out(ind2(n))=out(ind2(n))+a(ind1(ind(n)));
	cnt(ind2(n))=cnt(ind2(n))+1;
      end
    end
  end
  ind=find(cnt>1);
  out(ind)=out(ind)./cnt(ind);
   
  % fill any remaining nan pixels with nearest neighbor 
  ind1=find(isnan(out));
  [i,j]=ind2sub([rhead(2),rhead(1)],ind1);
  [lon,lat]=pix2latlon(i,j,rhead);
  [x,y]=latlon2pix(lon,lat,ahead);
  xc=round(x); yc=round(y);
  ind=find(xc>0 & xc<=ahead(1) & yc>0 & yc<=ahead(2));
  ind2=sub2ind([ahead(2),ahead(1)],xc(ind),yc(ind));
  out(ind1(ind))=a(ind2);
    
end
out(isnan(out))=anodata;