function consolplot(basen,varargin)
%
% function consolplot(location,varargin)
%
% Script to plot consolidated databased iceberg positions in fancy manner
% with various options to enhance capability
%
% optional arguments:
%  'after'  : add "after" tag to title, use Fig 1
%  'before' : add "before" tag to title, use Fig 3
%  '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/

afterflag=0;
beforeflag=0;
pauseflag=0;
noplotflag=0;
printflag=0;
quietflag=0;
singleflag=0;
zoomflag=0;
prefix='';
minstart=0;
fig=1;  % figure number for iceberg track plot
if nargin > 1
  for k=1:length(varargin)
    if strcmp(varargin(k),'after')
      afterflag=1;
      prefix='after';
    end
    if strcmp(varargin(k),'before')
      beforeflag=1;
      prefix='before';
      fig=3;
    end
    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';  % list of colors
colorsc='cybkgkbybkgk'; % contrasting colors

% 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]);

    if ~noplotflag
      % plot landmap image and add coastline overlay
      myfigure(fig); % myfigure() can be replaced with figure()
      clf;
      colormap(gray)
      imagesc(1-map1,[0 1]); % show land map
      hold on;
      % plot coastline
      ind=find(apen==3);
      n=length(ind); 
      for k=1:n-1
	lon=alon(ind(k):ind(k+1)-1);
	lat=alat(ind(k):ind(k+1)-1);
	[x,y]=latlon2pix(lon,lat,head);
	plot(x,y,'g');
      end
    end
    
    % date and time of year
    year=floor(csv.date/1000);
    jd=mod(csv.date,1000);
    time=year+(jd-1)/365;

    % extract positions for each sensor
    ncol=length(csv.Properties.VariableNames);
    atmp=char(csv.Properties.VariableNames);
    ncnt=0;
    legh=[];
    legn=[];
    for icol=1:ncol
      sen=char(csv.Properties.VariableNames(icol));
      if sen(end)=='1'  % location
	ncnt=ncnt+1;
	iflag=1;
	switch sen
	  case 'nic_1'
	    lon=csv.nic_2;
	    lat=csv.nic_1;
	    
	    % check for and extract interpolation column
	    idx=find(all(ismember(atmp(:,1:length('nic_3')),'nic_3')'));
	    if length(idx)>0
	      interpflag=csv.nic_3;
	    else
	      interpflag=zeros(size(year));
	    end
	    
	    % check for and extract size columns
	    jdx=find(all(ismember(atmp(:,1:length('size_1')),'size_1')'));
	    if length(jdx)>0
	      isize=csv.size_1.*csv.size_2*(1.852).^2;
	    else
	      isize=zeros(size(year));
	    end
	  
	    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('NIC size (sq km)');
	      title(sprintf('size %s (%d)',file_name,length(csv.date)))
	    end
	  
	  case 'qscat_1'
	    lon=csv.qscat_2;
	    lat=csv.qscat_1;
	    isize=zeros(size(year));
	    idx=find(all(ismember(atmp(:,1:length('qscat_3')),'qscat_3')'));
	    if length(idx)>0
	      interpflag=csv.qscat_3;
	    else
	      interpflag=zeros(size(year));
	    end

	  case 'ers_1'
	    lon=csv.ers_2;
	    lat=csv.ers_1;
	    isize=zeros(size(year));
	    idx=find(all(ismember(atmp(:,1:length('ers_3')),'ers_3')'));
	    if length(idx)>0
	      interpflag=csv.ers_3;
	    else
	      interpflag=zeros(size(year));
	    end

	  case 'ascat_1'
	    lon=csv.ascat_2;
	    lat=csv.ascat_1;
	    isize=zeros(size(year));
	    interpflag=zeros(size(year));
	    idx=find(all(ismember(atmp(:,1:length('ascat_3')),'ascat_3')'));
	    if length(idx)>0
	      interpflag=csv.ascat_3;
	    else
	      interpflag=zeros(size(year));
	    end

	  case 'nscat_1'
	    lon=csv.nscat_2;
	    lat=csv.nscat_1;
	    isize=zeros(size(year));
	    idx=find(all(ismember(atmp(:,1:length('nscat_3')),'nscat_3')'));
	    if length(idx)>0
	      interpflag=csv.nscat_3;
	    else
	      interpflag=zeros(size(year));
	    end

	  case 'oscat_1'
	    lon=csv.oscat_2;
	    lat=csv.oscat_1;
	    isize=zeros(size(year));
	    idx=find(all(ismember(atmp(:,1:length('oscat_3')),'oscat_3')'));
	    if length(idx)>0
	      interpflag=csv.oscat_3;
	    else
	      interpflag=zeros(size(year));
	    end

	  case 'seawinds_1'
	    lon=csv.seawinds_2;
	    lat=csv.seawinds_1;
	    isize=zeros(size(year));
	    idx=find(all(ismember(atmp(:,1:length('seawinds_3')),'seawinds_3')'));
	    if length(idx)>0
	      interpflag=csv.seawinds_3;
	    else
	      interpflag=zeros(size(year));
	    end

	  case 'sass_1'
	    lon=csv.sass_2;
	    lat=csv.sass_1;
	    isize=zeros(size(year));
	    idx=find(all(ismember(atmp(:,1:length('sass_3')),'sass_3')'));
	    if length(idx)>0
	      interpflag=csv.sass_3;
	    else
	      interpflag=zeros(size(year));
	    end

	  case 'size_1'
	    iflag=0;
	  otherwise
	    iflag=0;
	    fprintf('Unrecognized sensor: %s\n',sen);
	end
      
	if iflag % only do if known sensor

	  [x,y]=latlon2pix(lon,lat,head);
	
	  if ~noplotflag % plot iceberg positions 
	    myfigure(fig) % myfigure() can be replaced with figure()
	    % filter out 0,0 locations from plot
	    index=find(lat==0 & lon==0);
	    x(index)=nan;
	    y(index)=nan;
	    legh(ncnt)=plot(x,y,colors(ncnt));
	    plot(x,y,[colors(ncnt) '*']);
	    legn{ncnt}=sen(1:end-2);
	    if length(interpflag)>0 plot(x(find(interpflag==1)),y(find(interpflag==1)),[colors(ncnt) 'o']);
	    end
	  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(ncnt),'*']);
	    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))));
	      switch(ncnt)
		case 1
		  set(T,'Color',[1,0,0]);
		case 2
		  set(T,'Color',[0,1,0]);
		case 3
		  set(T,'Color',[1,1,0]);
		case 4
		  set(T,'Color',[1,0,1]);
		case 5
		  set(T,'Color',[0,1,1]);
	      end
	    end
	  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);
      
	  icnt=length(find(interpflag==1));
	  pcnt=length(find(lat<0));
	  isize=isize(isize>0);
	  fprintf('%d %s %8s %d %4d-%3.3d to %4d-%3.3d (%4d %4d of %4d) size: %0.f-%.0f\n',ifcnt,file_name,sen(1:end-2),length(jnd),ymin1,dmin1,ymax1,dmax1,icnt,pcnt,length(csv.date),min(isize),max(isize));
      
	else
	  continue;
	end
      end
    end

    if ~noplotflag % add legend and title 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

      if length(legh)>0
	legend(legh,legn);
      end
      title(sprintf('%s %s %4d-%3.3d to %4d-%3.3d (%d)',prefix,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
