function statsplot(basen,varargin)
%
% function stats(location,varargin)
%
% Script to plot stats database iceberg positions in fancy manner
% with various options to enhance capability
%
% optional arguments:
%  'noplot' : do not generate (disable) plot if present
%  'pause'  : pause after each iceberg if present and plotting is enabled
%  'print'  : generate .png file if present and plotting is enabled
%  'quiet'  : do not include name of file being read if present
%  'single' : plot only specified file if present, otherwise do all files
%              in directory
%  'zoom'   : zoom into area of iceberg track if present and single enabled
%

% Written by D.G. Long, Nov 2017 www.scp.byu.edu
% Uses SIR file projection utilities loadsir, latlon2pix, myfigure
% available from the Scatterometer Climate Record Pathfinder www.scp.byu.edu
% ftp://ftp.scp.byu.edu/pub/software/matlab/

pauseflag=0;
noplotflag=0;
printflag=0;
quietflag=0;
singleflag=0;
zoomflag=0;
minstart=0;
fig=2;  % figure number for iceberg track plot
if nargin > 1
  for k=1:length(varargin)
    if strcmp(varargin(k),'pause')
      pauseflag=1;
    end
    if strcmp(varargin(k),'noplot')
      noplotflag=1;
    end
    if strcmp(varargin(k),'print')
      printflag=1;
    end
    if strcmp(varargin(k),'quiet')
      quietflag=1;
    end
    if strcmp(varargin(k),'single')
      singleflag=1;
    end
    if strcmp(varargin(k),'zoom')
      zoomflag=1;
    end
    if isnumeric(cell2mat(varargin(k)))
      minstart=cell2mat(varargin(k));
    end
  end
end

persistent map1; % used so that map file needs only be loaded once
persistent head;
if isempty(map1) | isempty(head) 
  % load Antarctic land mask .sir file
  % This file is available from the scatterometer climate record pathfinder
  % (SCP) www.scp.byu.edu  The direct URL is
  % ftp://ftp.scp.byu.edu/data/qscat/info/landmasks/qus-Ant.sir.lmask.gz
  fname='/auto/archive/cia/qus-Ant.sir.lmask';
  [map1 head]=loadsir(fname);
end

if ~noplotflag
  persistent ant; % used so that map file needs only be loaded once
  if isempty(ant)
    % load Antarctica coastline background
    % This file is available from the scatterometer climate record pathfinder
    % (SCP) www.scp.byu.edu  The direct URL is
    % ftp://ftp.scp.byu.edu/outgoing/data/coastlines/filt_antarctica_gshhs.dat.gz
    ant=load('/auto/archive/cia/filt_antarctica_gshhs.dat');
  end
  alon=ant(:,2);
  alat=ant(:,1);
  apen=ant(:,3);
  rscale=0.3;
end

% set colors for display
colors='brgcykbrgcyk';
colorsc='cybkgkbybkgk';

% determine what file(s) to load
if singleflag
  files=1;
  [basen,file_name,ext]=fileparts(basen);
else
  files = dir([basen '/*.csv']);
end

% load each iceberg file
ifcnt=0;
for file = files'
  %disp(file.name);
  ifcnt=ifcnt+1; % file counter
  if ifcnt>minstart  % option to start later in list
  
   xmin=1e10; ymin=1e10; xmax=-1e10; ymax=-1e10;
  
   if ~singleflag
    [dummy,file_name,ext]=fileparts(file.name);
   end
   if ~quietflag
    disp(['Loading ' basen '/' file_name ext]);
   end
   csv = readtable([basen '/' file_name ext]);

   % date and time of year
   year=floor(csv.date/1000);
   jd=mod(csv.date,1000);
   time=year+(jd-1)/365;

   lat=csv.lat;
   lon=csv.lon;
   isize=csv.size;
   %disp=csv.disp;
   %flags=csv.flags;
   mask=csv.mask;
   angles=csv.vel_angle;
  
   [x,y]=latlon2pix(lon,lat,head);
  
   if ~noplotflag
    % plot iceberg size vs time
    myfigure(5); % myfigure() can be replaced with figure()
    plot(time,isize,'b');
    hold on;plot(time,isize,'*'); hold off
    xlabel('year'); ylabel('size (sq km)');
    title(sprintf('size %s (%d)',file_name,length(csv.date)))
    % plot iceberg mask vs time
    myfigure(6); % myfigure() can be replaced with figure()
    plot(time,mask,'b');
    hold on;plot(time,mask,'*'); hold off
    xlabel('year'); ylabel('mask flag');
    title(sprintf('mask flag %s (%d)',file_name,length(csv.date)))
    % plot iceberg angle vs time
    myfigure(7); % myfigure() can be replaced with figure()
    plot(time,angles,'b');
    hold on;plot(time,angles,'*'); hold off
    xlabel('year'); ylabel('rotation angle (sq km)');
    title(sprintf('angle %s (%d)',file_name,length(csv.date)))
    drawnow
   end

   if ~noplotflag
    % plot landmap image and add coastline overlay
    myfigure(fig); % myfigure() can be replaced with figure()
    clf;
    colormap(gray)
    % plot land map
    imagesc(1-map1,[0 1]);
    hold on;
    % plot coastline
    ind=find(apen==3);
    n=length(ind); 
    for k=1:n-1
      llon=alon(ind(k):ind(k+1)-1);
      llat=alat(ind(k):ind(k+1)-1);
      [xx,yy]=latlon2pix(llon,llat,head);
      plot(xx,yy,'g');
    end

    % plot iceberg positions 
    % filter out 0,0 locations from plot
    index=find(lat==0 & lon==0);
    x(index)=nan;
    y(index)=nan;
    plot(x,y,colors(1));
    plot(x,y,[colors(1) '*']);
   end

   % clip positions to size of image
   knd=find(floor(x)>0 & floor(x)<=head(1) & floor(y)>0 & floor(y)<=head(2));
   x=x(knd);
   y=y(knd);
  
   xmin=min([min(x(x~=0 & y~=0)),xmin]);
   xmax=max([max(x(x~=0 & y~=0)),xmax]);
   ymin=min([min(y(x~=0 & y~=0)),ymin]);
   ymax=max([max(y(x~=0 & y~=0)),ymax]);
	
   % check to see if an iceberg position is contained within land mask
   % note: land mask can be inaccurate near iceshelves
   ind=sub2ind(size(map1),floor(y),floor(x));
   jnd=find(map1(ind)>0.5);
   if ~noplotflag % plot iceberg over land positions 
    if length(jnd)>0
      plot(x(jnd),y(jnd),[colorsc(1),'*']);
    end
    
    % add dates
    Txoff=0.02; Tyoff=0.02; % displacement of date string from position
    for n=1:length(x)
      T=text(x(n)+Txoff,y(n)+Tyoff,sprintf('%4d-%3.3d',year(knd(n)),jd(knd(n))));
      set(T,'Color',[1,0,0]);
    end
  
    ymin1=min(year(~isnan(lat)));
    dmin1=floor((min(time(~isnan(lat)))-ymin1)*365+0.0001);
    ymax1=max(year(~isnan(lat)));
    dmax1=floor((max(time(~isnan(lat)))-ymax1)*365+0.0001);
    
    % add to the position plot, zoom
    myfigure(fig) % myfigure() can be replaced with figure()
    if zoomflag && singleflag && abs(xmin)<1e5 && abs(xmin)<1e5
      del=max(floor(xmax-xmin+10),floor(ymax-ymin+10))/2;
      xmin=floor((xmax+xmin)/2-del);
      xmax=floor((xmax+xmin)/2+del);
      ymin=floor((ymax+ymin)/2-del);
      ymax=floor((ymax+ymin)/2+del);
      axis(floor([xmin xmax ymin ymax]));
    end    
    
    title(sprintf('%s %4d-%3.3d to %4d-%3.3d (%d)',file_name,year(1),jd(1),year(end),jd(end),length(csv.date)));
    axis square;
    axis off
  
    if printflag
      print([file_name '.png'],'-dpng');
    end
    
    if pauseflag
      fprintf('pausing');
      pause;
      fprintf('... continuing\n');
    end
   end
  end
end

