function plotBYUhr(speed, dir, lat, lon, land)
%
% plotBYUhr(speed, dir, lat, lon, land)
% sample matlab function to plot BYU high resolution wind fields
% BYU Aug 2003
%
% only the speed and dir arguments are required. additional
% arguments result in additional plots being generated
%
% use loadBYUhr to read arrays

% myfigure can be replaced with figure

% plot images in along-track/cross-track

myfigure(1)
colormap('gray');
ind=find(dir>180);
dir(ind)=dir(ind)-360;
imagesc(dir,[-180 180]);

% create a circular grayscale
cm=gray;
cm2=cm;
[ncm k]=size(cm);
for k=1:ncm/2
  cm2(k,:)=1-cm(2*k-1,:);
  cm2(ncm-k+1,:)=1-cm(2*k-1,:);
end
colormap(cm2);
colorbar;
xlabel('Alongtrack pixel (2.5 km/pixel)');
ylabel('Crosstrack pixel (2.5 km/pixel)');
title('Wind direction relative to north (deg)');
%axis off


myfigure(2)
imagesc(speed);
%colormap('gray');
colorbar;
xlabel('Alongtrack pixel (2.5 km/pixel)');
ylabel('Crosstrack pixel (2.5 km/pixel)');
title('Wind speed (m/s)');
%axis off

if nargin > 2

  myfigure(3)
  imagesc(lat);
  colorbar;
  xlabel('Alongtrack pixel (2.5 km/pixel)');
  ylabel('Crosstrack pixel (2.5 km/pixel)');
  title('Latitude (deg)');

  myfigure(4)
  imagesc(lon);
  colorbar;
  xlabel('Alongtrack pixel (2.5 km/pixel)');
  ylabel('Crosstrack pixel (2.5 km/pixel)');
  title('Longitude (deg)');

  % wind speed/dir arrows
  
  %myfigure(5)
  %quiver(lon,lat,speed.*cos(dir*pi/180),speed.*sin(dir*pi/180));

end

% generate an image of wind speed and direction on at lat/lon
% projection.  Using this simple mapping approach, some output
% pixels are not filled in.  To make the image look subjectively
% better, a median fill is added if the land mask array is passed in.

if nargin > 2  % only do if lat, lon are provided

  latmin=min(min(lat));
  latmax=max(max(lat));
  dellat=latmax-latmin;

  lonmin=min(min(lon));
  lonmax=max(max(lon));
  dellon=lonmax-lonmin;
  if dellon > 180  % switch from [-180 +180] to [0 360]
    ind=find(lon<0)
    lon(ind)=lon(ind)+360;
    lonmin=min(min(lon));
    lonmax=max(max(lon));
    dellon=lonmax-lonmin;
  end

  sc=40;  % pixels/deg
  M=floor(dellon*sc);
  N=floor(dellat*sc);
  latax=([1:N]-1)/sc+latmin;
  lonax=([1:M]-1)/sc+lonmin;

  indx=floor((lon-lonmin)*sc)+1;
  ind=find(indx>M);
  indx(ind)=M;
  indy=N-floor((lat-latmin)*sc)+1;
  ind=find(indy>N);
  indy(ind)=N;

  ind=sub2ind([N M],indy,indx);
  
  arr=zeros([N,M]);
  arr(ind)=speed;
  if nargin > 4  % median fill projection "holes" to look better
    arr=medianfill(arr);
    arr(ind)=arr(ind).*land;
  end
  myfigure(6)
  imagesc(lonax,latax,arr)
  colorbar;
  title('Wind Speed (m/s)');
  xlabel('Longitude (deg)');
  ylabel('Latitude (deg)');
%  axis off

  arr=zeros([N,M]);
  arr(ind)=dir;
  if nargin > 4  % median fill projection "holes" to look better
    arr=medianfill(arr);
    arr(ind)=arr(ind).*land;
  end
  myfigure(7)
  imagesc(lonax,latax,arr)
  colorbar;
  title('Wind Direction (deg)');
  xlabel('Longitude (deg)');
  ylabel('Latitude (deg)');
%  axis off

end


function b=medianfill(a)
%
% a=medianfill(a)
%
% fills zero-valued "holes" in the array a with median
% of nearby values, assuming at least 3 nearby values are non-zero
%
b=a;
[m n]=size(a);
ind=find(a==0);
[x y]=ind2sub([m n], ind);
ind1=find(x>1 & x<m & y>1 & y < n);
ind=ind(ind1);
N=length(ind);
if N < 1
  return
end
indm=[-m-1 -m -m+1 -1 1 m-1 m m+1];
for j=1:N
  ind2=ind(j)+indm;
  c=a(ind2);
  i=find(c ~=0);
  if length(i) > 2
    b(ind(j))=median(c(i));
  end
end










