% copyright 2018, David G. Long at BYU
%
% This matlab script illustrates reading and plotting of the contents of
% some of the key data arrays in the BYU-standard L2H UHR wind product
% files for ASCAT.  These files are in netcdf form, and can be read by 
% built in Matlab utilties.  Here reading and display for most arrays is
% illustrated with the idea that the user will keep and use only those of
% particular interest for their application.
%

clear

%% define input file name
fname='/auto/temp/brown/ascat_winds/output/hawaii_winds2011/Data/360/msfa-a-E2T11-360-360-26903.l2h.nc';
fprintf('reading and plotting file %s\n',fname);

%% begin reading .nc file
if 0 % optionally read SIR file header info
  % this information can be useful in map transformations
  % (not used in this script)
  sir_geom_nsx = double(ncreadatt(fname, '/', 'sir_geom_nsx'));
  sir_geom_nsy = double(ncreadatt(fname, '/', 'sir_geom_nsy'));
  sir_geom_iopt = double(ncreadatt(fname, '/', 'sir_geom_iopt'));
  sir_geom_xdeg = double(ncreadatt(fname, '/', 'sir_geom_xdeg'));
  sir_geom_ydeg = double(ncreadatt(fname, '/', 'sir_geom_ydeg'));
  sir_geom_ascale = double(ncreadatt(fname, '/', 'sir_geom_ascale'));
  sir_geom_bscale = double(ncreadatt(fname, '/', 'sir_geom_bscale'));
  sir_geom_a0 = double(ncreadatt(fname, '/', 'sir_geom_a0'));
  sir_geom_b0 = double(ncreadatt(fname, '/', 'sir_geom_b0'));
end

% read and plot selected arrays for the UHR wvc pixels
% note: all pixels will have lat/lon locations, even if empty
% use "namb>0" to determine which wvc pixels have wind estimates
namb=double(ncread(fname,'number_ambiguities'));
ambselect=double(ncread(fname,'ambiguity_select'));
lat=double(ncread(fname,'lat'));
lon=double(ncread(fname,'lon'));
wspds=double(ncread(fname,'wspeeds'));
wdirs=double(ncread(fname,'wdirs'));
%mles=double(ncread(fname,'mles'));

% determine selected wind vector
N=size(wspds);
wspd=zeros(size(N(2),N(3)));
wdir=zeros(size(N(2),N(3)));
%mle=zeros(size(N(2),N(3)));
for x=1:N(2)
  for y=1:N(3)
    wspd(x,y)=wspds(max(1,ambselect(x,y)),x,y);
    wdir(x,y)=wdirs(max(1,ambselect(x,y)),x,y);
    %mle(x,y)=mles(max(1,ambselect(x,y)),x,y);
  end
end
% save space
clear wspds wdirs

% plot selected UHR winds
myfigure(1);
imagesc(rot90(wspd));; colorbar;
title('wind speed (m/s)');
axis off
myfigure(2);
imagesc(rot90(wdir),[-180 180]);; colorbar;
axis off
title('wind direction (m/s)');

if 0  % Optionally create 
  myfigure(3);
  colormap('gray');
  imagesc(rot90(wdir),[-180 180]);

  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;
  axis off
  title('wind direction grayscale (m/s)');
end
  
if 1  % optionally read and plot land mask
  land=double(ncread(fname,'land'));
  myfigure(4)
  imagesc(rot90(land));
  title('Land mask (blue=ocean, yellow=land)');
  axis off;
  clear land
end

if 1  % optionally plot lat/lon. 
  % Note 90 deg rotation so image "looks" right when displayed
  myfigure(5);
  imagesc(rot90(lat)); colorbar;
  title('Latitude');
  axis off;
  myfigure(6);
  imagesc(rot90(lon)); colorbar;
  title('Longitude');
  axis off;
end


if 1 % optionally read and plot low res (L2B) winds (from L2B file)
  lspd=ncread(fname,'L2B_speed');
  myfigure(7);
  imagesc(rot90(lspd)); colorbar;
  axis off;
  title('L2B model speed (m/s)');
  ldir=ncread(fname,'L2B_dir');
  myfigure(8);
  imagesc(rot90(ldir),[-180 180]); colorbar;
  axis off;
  title('L2B model dir (deg)');
  clear lspd ldir
end

if 1 % optionally read and plot model winds (originally from L2B file)
  mspd=ncread(fname,'L2B_Mspeed');
  myfigure(9);
  imagesc(rot90(mspd)); colorbar;
  axis off;
  title('Model speed (m/s)');
  mdir=ncread(fname,'L2B_Mdir');
  myfigure(10);
  imagesc(rot90(mdir),[-180 180]); colorbar;
  axis off;
  title('Model dir (deg)');
  clear mspd mdir
end

if 1 % optionally plot number of ambiguities and ambiguity selection
  myfigure(11)
  imagesc(rot90(namb),[0 4]); colorbar;
  title('Number of UHR ambiguities');
  axis off;
  myfigure(12)
  imagesc(rot90(ambselect),[0 4]); colorbar;
  title('UHR ambiguity selected');
  axis off;
end

if 1  % optionally read and plot distance to nearest land value
  ldis=double(ncread(fname,'ldis'));
  myfigure(13)
  imagesc(rot90(ldis)); colorbar;
  title('Distance to land (km) (clipped to <256)');
  axis off;
  clear ldis
end

if 1  % optionally read and plot pixel time
  time=double(ncread(fname,'time'));
  time0=min(time(time>1));
  time(time<time0)=time0-1;
  myfigure(14)
  imagesc(rot90(time-time0)); colorbar;
  title(sprintf('Measurement time (seconds from %d\n seconds from 2000-01-01 00:00:00)',time0));
  axis off;
  clear time
end

if 1 % optionally plot sigma-0, inc, and azimuth angles
  sig=ncread(fname,'sig');
  myfigure(15);
  subplot(2,2,1)
  imagesc(rot90(squeeze(sig(1,:,:))),[-30,-5]); colorbar;
  axis off; axis image;
  title('fore sigma-0 (dB)');
  subplot(2,2,2)
  imagesc(rot90(squeeze(sig(2,:,:))),[-30,-5]); colorbar;
  axis off; axis image;
  title('mid sigma-0 (dB)');
  subplot(2,2,3)
  imagesc(rot90(squeeze(sig(3,:,:))),[-30,-5]); colorbar;
  axis off; axis image;
  title('aft sigma-0 (dB)');
  clear sig

  azi=ncread(fname,'azi');
  myfigure(16);
  subplot(2,2,1)
  imagesc(rot90(squeeze(azi(1,:,:))),[-180,180]); colorbar;
  axis off; axis image;
  title('fore azimuth angle (deg)');
  subplot(2,2,2)
  imagesc(rot90(squeeze(azi(2,:,:))),[-180,180]); colorbar;
  axis off; axis image;
  title('mid azimuth angle (deg)');
  subplot(2,2,3)
  imagesc(rot90(squeeze(azi(3,:,:))),[-180,180]); colorbar;
  axis off; axis image;
  title('aft azimuth angle (deg)');
  clear azi

  inc=ncread(fname,'inc');
  myfigure(17);
  subplot(2,2,1)
  imagesc(rot90(squeeze(inc(1,:,:))),[0,65]); colorbar;
  axis off; axis image;
  title('fore incidence angle (deg)');
  subplot(2,2,2)
  imagesc(rot90(squeeze(inc(2,:,:))),[0,65]); colorbar;
  axis off; axis image;
  title('mid incidence angle (deg)');
  subplot(2,2,3)
  imagesc(rot90(squeeze(inc(3,:,:))),[0,65]); colorbar;
  axis off; axis image;
  title('aft incidence angle (deg)');
  clear inc
end

