function [latArray,lonArray]=...
		plotSingleImage(imageArray, latArray, lonArray, figNum, ...
                         keywords, latLonLimits, resolutionDifference)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% function plotSingleImage(imageArray, latArray, lonArray, figNum, ...
%%                         keywords, latLonLimits, resolutionDifference)
%%
%% plotSingleImage uses matlab's mapping toolbox to project imageArray
%% of data onto the earth using latArray and lonArray as lat/lon. coastlines 
%% are shown globally, and state boundaries are shown within the US.  
%%
%% Suzanne Wetzel Seemann
%% CIMSS/SSEC
%% swetzel@ssec.wisc.edu
%% October 2001 
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%% NOTE: if you want to specify resolutionDifference but you want to 
%%       use the default latLonLimits, you may have 6 input arguments
%%       with the last being resolutionDifference and no latLonLimits.
%%       the same is NOT true for "keywords" -- a cellstring must be 
%%       provided as the fifth argument if the 6th and/or 7th are needed,
%%       even if no keyword options are desired.  this can be any dummy
%%       cellstring, e.g. {'none'}, {'dummy'}, {'etc'}  
%%
%% required input arguments:
%%      imageArray      array of data to be projected onto a map
%%      latArray        array of latitudes corresponding to imageArray
%%      lonArray        array of longitudes corresponding to imageArray
%%      figNum          figure number for plot
%%
%% optional input arguments:
%%      keywords        cellstring of keywords to specify various options
%%                      (a cellstring is {'keyword1', 'keyword2', etc.} 
%%                      see list of keywords and their meanings below.
%%                      if no keywords are needed, but input arguments
%%                      6 or 7 are required, use {'none'} or other dummy 
%%                      argument for keywords 
%%
%%      latLonLimits    optional limits of lat/lon for the map.  should be
%%                      in degrees and of the format: 
%%                            [southLat northLat westLon eastLon]
%%                      if this argument is not specified, the max/min 
%%                      latitudes and longitudes of the data are used
%%
%%      resolutionDifference  where the lat/lon arrays are at a lower
%%                            resolution than the data arrays (such as
%%                            modis output from shaima's modis_level1b_read),
%%                            resolutionDifference must be specified:  
%%       (size of lat/lon arrays)*resolutionDifference = (size of data arrays) 
%%                            at present, this assumes the same factor for
%%                            both lat and lon and in all dimensions.  
%%
%%
%% LIST OF KEYWORDS CURRENTLY SUPPORTED:
%%
%%
%%  note on keywords: 1) if your specific application requires a projection that
%%                    is not listed here, a color that is not currently used,
%%                    a map axes option that is not here, or anything else
%%                    different, let me know.  the current keywords include the
%%                    common cases but are easy to tailor to specific cases. 
%%
%%                    2) keywords are case insensitive
%%    
%%                    3) use cellstring format.  
%%                       e.g., {'gridLines','colorbar','mercator'} 
%%
%%  MISCELLANEOUS KEYWORDS
%%
%%         'gridLines'  includes lat/long gridlines on map (note this sometimes
%%                      alters the display so the map occupies less real estate
%%                      on the figure window, so try it on and off). default no
%%                      gridlines.  gridline color is the same as coastline 
%%                      color.  if no coastline color is specified both coasts 
%%                      and grids are black.
%%
%%         'compassLabels' lat/lon line labels are displayed as N,S,E,W instead
%%                         of the default signed labels (+ for N/E, - for S/W)
%% 
%%         'noMap'  supresses plotting of the geographic boundaries
%%         'overwrite'  allows for plotting multiple granules on the same
%%                      figure.  when this keyword is selected, the figure is 
%%                      not cleared and to save time the coastlines, etc. are 
%%                      not redrawn.  default no overwrite. 
%%          
%%         'invertColorMap'  colormap is currently set to "gray" with the 
%%                           default being high values white and low values 
%%                           black. 'invertColorMap' keyword will invert the 
%%                           colormap to be more realistic for IR images. 
%%
%%         'discreteColorMap20' sets colormap to have 20 discrete colors 
%%         'discreteColorMap15' sets colormap to have 15 discrete colors 
%%         'discreteColorMap10' sets colormap to have 10 discrete colors 
%%         'discreteColorMap5' sets colormap to have 5 discrete colors 
%%
%%         'jetColorMap'    sets colormap to "jet" instead of the default "gray"
%%         'winterColorMap' sets colormap to "winter" instead of the default "gray"
%%         'boneColorMap' sets colormap to "bone" instead of the default "gray"
%%         'mcidasColorMap' uses the mcidas colormap from goes sounder products
%%         'cloudMaskColorMap'  
%%
%%         'colorbar'   plots the colorbar on the figure to identify the values
%%                      that correspond to a given color
%%
%%         'frame'      draws a frame around the map
%%         'hidem'      makes background all grey
%%
%%         'centerProjOrigin'  centers the projection origin around the
%%                             map lon.  default is projection origin
%%                             centered around lat = 0, lon = 0, z = 0
%%                             (needed where the image crosses +/- 180)
%%         'highResMap' uses a high resolution coastline map
%%
%%  COLOR FOR COASTLINE and STATELINE: DEFAULT is BLACK 
%%         'magentaMap'  
%%         'blueMap'
%%         'yellowMap'
%%         'whiteMap'
%%         'redMap'
%%    
%%  PROJECTION: DEFAULT is 'eqaconic', Equal Area Conic (Albers)
%%         'mercator'  Mercator Cylindrical
%%         'lambcyln'  Lambert Cylindrical
%%         'lambert'   Lambert Conformal Conic 
%%         'mollweid'  Mollweid
%%         'robinson'  Robinson
%%         'miller'    Miller Cylindrical
%%         'gortho'    Gall Orthographic
%%         'gstereo'   Gall Stereographic
%%         'giso'      Gall Isographic
%%         'eqdcylin'  Equidistant Cylindrical
%%         'eqacylin'  Equal Area Cylindrical
%%         'eqdconic'  Equidistant Conic  
%%         'aitoff'    Aitoff
%%	   'eckert3'   Eckert III
%%         'eckert4'   Eckert IV
%%         'loximuth'  Loximuthal
%%         'bsam'      Bolshoi Sovietskii Atlas Mira 
%%         'cassini'   Cassini Cylindrical Projection
%%
%% notes: 1) to get a list of other projections supported by the 
%%       matlab mapping toolbox that could be included here
%%       (but are not currently), type 'getm mapprojection'.  the
%%       currrent code is not set up for azimuthal projections, as
%%       they use polar coordinates instead of lat/lon.
%%       
%%       2) to see a global example of most of these projections,
%%       use the matlab gui 'viewmaps' 
%%
%% sample run statements:
%%   plotSingleImage(imageArray, latArray, lonArray, 1);
%%
%%   plotSingleImage(imageArray, latArray, lonArray, 1, ...
%%       {'invertColorMap', 'magentaMap', 'mercator'});
%%
%%   plotSingleImage(imageArray, latArray, lonArray, 1, ...
%%       {'saveData','gridLines','whiteMap','mcidasColorMap'}, 5);
%%
%%   plotSingleImage(imageArray, latArray, lonArray, 1, ...
%%       {'gridLines','colorbar','invertColorMap','yellowMap', ...
%%        'robinson'}, [25 40 -90 -60]);
%%
%%   plotSingleImage(imageArray, latArray, lonArray, 1, ... 
%%       {'colorbar','invertColorMap','magentaMap','miller'}, ...
%%       [25 40 -90 -60], 5); 
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Error checking and processing of input
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if isa(latArray,'double') == 0 | isa(lonArray,'double') == 0
   latArray = double(latArray);
   lonArray = double(lonArray);
end
 
if isa(imageArray,'double') == 0 
   imageArray = double(imageArray);
end

nargs = nargin;
if nargs < 4 
  error(['Function plotSingleImage requires at least 4 input arguments.']);
else
  %% error check figNum
  if (~isnumeric(figNum) | (length(figNum) ~= 1))
      error(['figNum argument must be a single numeric']);
  end
  %if (figNum - round(figNum) ~= 0)
  %    error(['figNum must be an integer']);
  %end
end

if nargs > 4 
  if ~iscellstr(keywords)
     error(['Fifth argument must be a cellstring of keywords: ' ...
            ' {''keyword1'',''keyword2'', etc...}']);
  else
     keywords = lower(keywords);
  end
else
  keywords = {' '};
end

resdif = 'n';
latlon = 'n';

if nargs == 6 
    if ~isnumeric(latLonLimits)
       error(['sixth argument must be numeric']);
    end

    %% decide whether latLonLimits or resolutionDifference is first
    if length(latLonLimits) == 1 
       resolutionDifference = latLonLimits;
       latLonLimits = [];
       resdif = 'y';
  
    elseif length(latLonLimits) == 4
       latlon = 'y';

       %% check that imageArray, latArray, lonArray are the same size because
       %% there is no resolutionDifference argument
       if (size(latArray) ~= size(lonArray)) | ...
          (size(latArray) ~= size(imageArray))
          error(['latArray and lonArray must be the same size as ' ...
                 'imageArray if no resolutionDifference argument is ' ...
                 'specified']);
       end
    else
       error(['Sixth input argument must either be a single integer for ' ...
              'resolutionDifference or a 4 element array for latLonLimits.']); 
    end
end
 
%% if 7 input arguments, set both resdif and latlon to 'y', indicating
%% both resolutionDifference and latLonLimits were specified 

if nargs == 7 
  resdif = 'y';
  latlon = 'y';
end 

if nargs > 7
  error(['plotSingleImage accepts a maximum of 7 input arguments']);
end 

%% if latLonLimits specified, check latLonLimit input

if strcmp(latlon,'y')
  if max(abs(latLonLimits(1:2))) > 90  | max(abs(latLonLimits(3:4))) > 180
     error(['Latitude cannot exceed +/- 90 deg and Longitude cannot ' ...
            'exceed +/- 180 deg.']); 
  end
  if (latLonLimits(1) > latLonLimits(2)) | (latLonLimits(3) > latLonLimits(4))
     error(['latLonLimits must be of the form: [southLat northLat westLon ' ...
            'eastLon]']);
  end

else
  %% create latLonLimits if they're not given as input
  %% check for lat,lon == -999 first and don't include this 
  %% value as minimum lat/lon 
   minLon = min(lonArray(:));
   minLat = min(latArray(:));

   badIndsLon = find(lonArray == -999);
   badIndsLat = find(latArray == -999);

   if ~isempty(badIndsLon)
        temp = lonArray;
        temp(badIndsLon) = 999;
        minLon = min(temp(:));
   end
   if ~isempty(badIndsLat)
        temp = latArray;
        temp(badIndsLat) = 999;
        minLat= min(temp(:));
   end

   latLonLimits = [minLat max(latArray(:)) minLon max(lonArray(:))];
end

%% if resolutionDifference specified, interpolate lat/lon data by the 
%% given factor 

if strcmp(resdif,'y') 
  if (resolutionDifference - round(resolutionDifference) ~= 0)
      error(['resolutionDifference must be an integer']);
  end
   
   %% interpolate latitude and longitude to increase size by a factor 
   %% of resolutionDifference

   disp([' ']);
   disp(['Interpolating latitude and longitude to increase size of arrays ' ...
         'by a factor of ' num2str(resolutionDifference) '...']);
 
% shaima's code (didn't work)
%   latArray = GetInterpArray(latArray, resolutionDifference, ...
%				      size(imageArray));
%   lonArray = GetInterpArray(lonArray, resolutionDifference, ...
%				      size(imageArray));

% my quick 'n dirty approach:
    [X,Y] = meshgrid([1:5:5*size(latArray,2)],[1:5:5*size(latArray,1)]);
    XI = [1:5*size(latArray,2)];
    YI = [1:5*size(latArray,1)]';
    latArray = interp2(X,Y,latArray,XI,YI);
    lonArray = interp2(X,Y,lonArray,XI,YI);
end

figure(figNum);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% set projection based on keywords
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
numProj = 0;

if ismember('mercator',keywords)
   projection = 'mercator';
   numProj = numProj + 1;   
end
if ismember('lambcyln',keywords)
   projection = 'lambcyln';
   numProj = numProj + 1;   
end
if ismember('lambert',keywords)
   projection = 'lambert';
   numProj = numProj + 1;   
end
if ismember('mollweid',keywords)
   projection = 'mollweid';
   numProj = numProj + 1;   
end
if ismember('robinson',keywords)
   projection = 'robinson';
   numProj = numProj + 1;   
end
if ismember('miller',keywords)
   projection = 'miller';
   numProj = numProj + 1;   
end
if ismember('gortho',keywords)
   projection = 'gortho';
   numProj = numProj + 1;   
end
if ismember('gstereo',keywords)
   projection = 'gstereo';
   numProj = numProj + 1;   
end
if ismember('giso',keywords)
   projection = 'giso';
   numProj = numProj + 1;   
end
if ismember('eqdcylin',keywords)
   projection = 'eqdcylin';
   numProj = numProj + 1;   
end
if ismember('eqacylin',keywords)
   projection = 'eqacylin';
   numProj = numProj + 1;   
end
if ismember('eqdconic',keywords)
   projection = 'eqdconic';
   numProj = numProj + 1;   
end
if ismember('aitoff',keywords)
   projection = 'aitoff';
   numProj = numProj + 1;   
end
if ismember('eckert3',keywords)
   projection = 'eckert3';
   numProj = numProj + 1;   
end
if ismember('eckert4',keywords)
   projection = 'eckert4';
   numProj = numProj + 1;   
end
if ismember('loximuth',keywords)
   projection = 'aitoff';
   numProj = numProj + 1;   
end
if ismember('bsam',keywords)
   projection = 'aitoff';
   numProj = numProj + 1;   
end
if ismember('cassini',keywords)
   projection = 'cassini';
   numProj = numProj + 1;   
end

if numProj > 1
  error(['only one projection can be specified in keywords']);
elseif numProj < 1
  projection = 'eqaconic';
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% define axes if no overwrite
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if ismember('magentamap',keywords)
  mapColor = 'm';
elseif ismember('bluemap',keywords)
  mapColor = 'b';
elseif ismember('redmap',keywords)
  mapColor = 'r';
elseif ismember('yellowmap',keywords)
  mapColor = 'y';
elseif ismember('whitemap',keywords)
  mapColor = 'w';
else
  mapColor = 'k';
end

if ~ismember('overwrite',keywords)
  %cla; clf;
  if ismember('gridlines',keywords)  
      numberOfGridLines = 4;

      deltaLat = latLonLimits(2) - latLonLimits(1);
      deltaLon = latLonLimits(4) - latLonLimits(3);
     
      lineSeparation = max(ceil(deltaLon/numberOfGridLines),ceil(deltaLat/numberOfGridLines));
      if lineSeparation > 5 
          lineSeparation = roundn(lineSeparation,1);
      end
      if lineSeparation > 30
          lineSeparation = 30;
      end
      axesm('mapprojection', projection, 'maplatlimit', latLonLimits(1:2), ...
            'maplonlimit', latLonLimits(3:4), 'grid', 'on', 'gcolor', ...
            mapColor, 'meridianlabel', 'on', 'parallellabel', 'on', ...
            'mlinelocation',lineSeparation,'plinelocation',lineSeparation,...
            'mlabellocation',lineSeparation,'plabellocation',lineSeparation);
  else 
      axesm('mapprojection', projection, 'maplatlimit', latLonLimits(1:2), ...
            'maplonlimit', latLonLimits(3:4), 'grid', 'off');
  end
end

if ismember('compasslabels',keywords)
  setm(gca,'labelformat','compass');
else
  setm(gca,'labelformat','signed');
end

if ismember('frame',keywords)
  setm(gca, 'frame', 'on');
end

if ismember('centerprojorigin',keywords)
  midLon = mean([latLonLimits(3:4)]);  
  setm(gca,'origin',[0 midLon 0]);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% plot data 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% note: time testing showed that "quicklook" actually
% doubled the time to view for a 596x1249 array.  
% reduceNbyN needs to be optimized
%if ismember('quicklook',keywords)
   %% call shaima's code to reduce size
%   avgnum = 5; 
%   latArray = ReduceNbyN(latArray, avgnum);
%   lonArray = ReduceNbyN(lonArray, avgnum);
%   imageArray = ReduceNbyN(imageArray, avgnum);
%end

%size(latArray)
%size(imageArray)
surfm(latArray,lonArray,imageArray);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% colormap  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
%% set number of colors
if ismember('discretecolormap20',keywords)
   numColors = 20;
elseif ismember('discretecolormap15',keywords)
   numColors = 15;
elseif ismember('discretecolormap10',keywords)
   numColors = 10;
elseif ismember('discretecolormap5',keywords)
   numColors = 5;
else
   numColors = 264;
end

%% set color scheme
if ismember('jetcolormap',keywords)
   colormap(jet(numColors));
elseif ismember('wintercolormap',keywords)
   colormap(winter(numColors));
elseif ismember('bonecolormap',keywords)
   colormap(bone(numColors));
elseif ismember('mcidascolormap',keywords)
   cmap=eu('rest','GA20P.ET');
   colormap(cmap);
elseif ismember('cloudmaskcolormap',keywords)
   cmap = [1.0 1.0 1.0; 1.0 1.0 0; 0 0.8 0.4; 0 0 1.0];
   colormap(cmap); 
   caxis([0 3]);
else
   colormap(gray(numColors))
end

if ismember('invertcolormap',keywords)
   cc = colormap; colormap(flipud(cc));
end
 
% attempt to make discrete colormap bounds to have standard endpoints:
cax = caxis;
deltaCax = (cax(2)-cax(1))/numColors;
if deltaCax >= .01 & deltaCax <= 1000 & numColors ~= 264
  intervalChoices=[0.01 0.02 0.05 0.1 0.2 0.5 1 2 5 10 20 50 100 200 500 1000];
  iii = intervalChoices(find(intervalChoices - deltaCax > 0) - 1); 
  caxMin = iii(1);
  caxMax = caxMin + caxMin*(numColors);
  caxis([caxMin caxMax]);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% draw maps if no overwrite and if no noMap
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if ~ismember('overwrite',keywords) & ~ismember('nomap',keywords)
	if ismember('highresmap',keywords)
           % need to do this only once
 %          gshhs('gshhs_h.b','createindex');
           s = gshhs('gshhs_h.b',latLonLimits(1:2), latLonLimits(3:4));
	   % this line to have filled continents
	   %displaym(s); tightmap;
	   % this line to have continents outlined
	   [mapLat, mapLon] = extractm(s); 
	   h = plotm(mapLat,mapLon,1,mapColor); tightmap; 
	   set(h,'linewidth',0.5);
        else
           load coast; plotm(lat, long, 1, mapColor);
           load usalo; plotm(statelat,statelon,mapColor); plotm(uslat,uslon,mapColor);
           plotm(gtlakelat,gtlakelon, 1, mapColor); tightmap;
        end
%xlabel('<-- W                                                          E -->');
%ylabel('<-- S                                                          N -->');

  %fillm(gtlakelat,gtlakelon,1,'facecolor','cyan','edgecolor','none');
end 

if ismember('colorbar',keywords)
  colorbar;
end

if ismember('hidem',keywords)
  hidem(gca);
end

if ismember('cloudmaskcolormap',keywords)
  a = get(gcf,'Children');
  set(a(2),'ytick',[0.45 1.18 1.8 1.9 2.6 2.7]);
  set(a(2),'yticklabel',{'Cloudy','Uncertain','Clear','Probably','Clear','Confident'});
  set(a(2),'ticklength',[0 0]);
end

