function hout = worldmap(varargin)

%WORLDMAP maps a country or region using the WORLDLO atlas data
%  
%   WORLDMAP maps a region or country using the WORLDLO atlas data. The region 
%   is selected interactively from a list. The map is created in the current
%   axes. 
%  
%   WORLDMAP REGION or WORLDMAP('REGION') makes a map of the specified 
%   region.  Recognized region strings are 'Africa', 'Antarctica', 'Asia', 
%   'Europe',.  'North America', 'North Pole', 'Pacific', 'South America', 
%   'South Pole', 'World' and the names of the countries in the WORLDLO atlas 
%   data.  REGION may also be a padded string matrix or a cell array of strings 
%   containing multiple country names.
%   
%   WORLDMAP REGIONONLY adds only the specified country to the map.  For 
%   example, WORLDMAP italyonly.  If any of the country names in REGION string 
%   matrix or cell array of strings end with 'only', only the requested 
%   countries are displayed.
%  
%   WORLDMAP('REGION','TYPE') controls the atlas data displayed. Type 'line' or
%   'patch' creates a map with atlas data of those types. If the size of the map
%   permits, point and text annotation of countries and cities is included'. 
%   Type 'lineonly' or 'patchonly' suppresses the point and text annotation.
%   Type 'none' suppresses all atlas data. If omitted, type 'line' is assumed.
%  
%   WORLDMAP(latlim,lonlim) and WORLDMAP(latlim,lonlim,'TYPE') use the supplied
%   geographic limits to define the map. The geographic limits are two-element 
%   vectors of the form [start end], in angular units of degrees. 
%  
%   WORLDMAP(map,maplegend) and WORLDMAP(map,maplegend,'TYPE') use the supplied
%   regular matrix map to define the geographic limits. The matrix map is displayed
%   using MESHM unless type is 'none', 'lineonly', 'patch' or 'patchonly'. Use type
%   'meshonly' to display only the matrix map.
%  
%   h = WORLDMAP(...) returns the handle of the map axes.
%  
%   WORLDMAP uses TIGHTMAP set the axis limits tight around the map. If you change 
%   the projection, or just want more white space around the map frame, use TIGHTMAP
%   again or AXIS AUTO.
%  
%   See also USAMAP, GRIDM, MLABEL, PLABEL, FRAMEM, WORLDLO

%  Copyright 1996-1998 by Systems Planning and Analysis, Inc. and The MathWorks, Inc.
%  $Revision: 1.5 $ $Date: 1998/10/01 19:51:13 $ 
%  Written by:  W. Stumpf


% Parse and check input arguments

if nargin == 0; % EASYMAP
	region = [];
	type = 'line';
	latlim = [];
	lonlim = [];
elseif nargin ==1 % EASYMAP REGION
	region = varargin{1};
	type = 'line';
	latlim = [];
	lonlim = [];
elseif nargin ==2 & ( (isstr(varargin{1}) | iscell(varargin{1})) & isstr(varargin{2})) % EASYMAP(region,type)
	region = varargin{1};
	type = varargin{2};
	latlim = [];
	lonlim = [];
elseif nargin ==2 & ( ~isstr(varargin{1}) & ~isstr(varargin{2})) % EASYMAP(latlim,lonlim),  EASYMAP(map,maplegend)
	region = [];
	latlim = varargin{1};
	lonlim = varargin{2};
	type = 'line';
	if all( size(lonlim) == [1 3]) & size(latlim,1) > 1
		map = varargin{1};
		maplegend = varargin{2};
		[latlim,lonlim] = limitm(map,maplegend);
		lonlim = lonlim + [-1 0]*(10*epsm('deg')); % avoids trimming
		type = 'mesh';
	end	
elseif nargin ==3 & ( ~isstr(varargin{1}) & ~isstr(varargin{2})) % EASYMAP(latlim,lonlim,type),  EASYMAP(map,maplegend,type)
	region = [];
	latlim = varargin{1};
	lonlim = varargin{2};
	type = varargin{3};		
	if all( size(lonlim) == [1 3])
		map = varargin{1};
		maplegend = varargin{2};
		[latlim,lonlim] = limitm(map,maplegend);
		lonlim = lonlim + [-1 0]*(10*epsm('deg')); % avoids trimming
	end
else
	error('Incorrect number or type of arguments')
end

% ensure row vectors

latlim = latlim(:)';
lonlim = lonlim(:)';

if ~isempty(latlim);
	if ~isequal(size(latlim),[1,2]) & ~isequal(size(lonlim),[1,2])
		error('Latlim and lonlim must be two element vectors')
	end
end

if ~isempty(type)
	typeindx = strmatch(type,{'patch','line','patchonly','lineonly','none','mesh','meshonly'},'exact');
	if isempty(typeindx); error(sprintf('Map type is not ''patch'', ''line'', ''patchonly'', ''lineonly'', \n''mesh'', ''meshonly''  or ''none''. '));end
end




onlyflag = 0;

% If no geographic limits or region supplied, interactive 
% selection of region

if nargin == 0;

	% Load the worldlo database. 
	load worldlo
	
	countrylist = {...
		'Africa',...
		'Antarctica',...
		'Asia',...
		'Europe',...
		'North America',...
		'North Pole',...
		'Pacific',...
		'South America',...
		'South Pole',...
		'World',...
		POpatch.tag};
	
	indx = listdlg('ListString',countrylist,...
                   'SelectionMode','single',...
					    'Name','Select a region');
						
	if isempty(indx); return; end
	
	region = countrylist{indx};
end

% If no geographic limits supplied, use region to look it up

if isempty(latlim)

	if ~iscell(region)
		switch lower(region)
		case {'af','africa'}
			latlim = [-40 40];
			lonlim = [-20 54];
		case {'an','antarctica'}
			latlim = [-90 -65];
			lonlim = [-180 180];
		case {'pa','pacific'}
			latlim = [-50 30];
			lonlim = [110 210];
		case {'sa','south america'}
			latlim = [-60 15];
			lonlim = [-90 -30];
		case {'e','europe'}
			latlim = [30 60];
			lonlim = [-15 45];
		case {'a','asia'}
			latlim = [-10 85];
			lonlim = [25 195];
		case {'na','north america'}
			latlim = [14 85];
			lonlim = [-170 -49];
		case {'np','north pole'}
			latlim = [84 90];
			lonlim = [-180 180];
		case {'sp','south pole'}
			latlim = [-60 -90];
			lonlim = [-180 180];
		case {'wo','world'}
			latlim = [-90 90];
			lonlim = [-180 180];
		end
	end
	
	if isempty(latlim) % country name from worldlo's POpatch

		
% detect the 'only case', 
		
		% single string 
		if isstr(region) & size(region,1) ==1 
		
			if length(region) >= 5 & strcmp(region(end-4:end),' only') 
				region = region(1:end-5);
				onlyflag=1;
			elseif length(region) >= 4 & strcmp(region(end-3:end),'only')
				region = region(1:end-4);
				onlyflag=1;
			end
		
		% string matrix - convert to cell array of strings, which need not be padded
		elseif isstr(region) & size(region,1) > 1
			for i=1:size(region,1)
				cellregion{i} = region(i,:);
			end
			region = deblank(cellregion);
		end
	
		%only case with cell array of strings
		if iscell(region)
			for i=1:length(region)
				thisregion = region{i};
				if length(thisregion) >= 5 & strcmp(thisregion(end-4:end),' only') 
					thisregion = thisregion(1:end-5);
					onlyflag=1;
				elseif length(thisregion) >= 4 & strcmp(thisregion(end-3:end),'only')
					thisregion = thisregion(1:end-4);
					onlyflag=1;
				end
				region{i} = thisregion;
			end
		end
	
% special abbreviations

		if ~iscell(region) & (strcmp(lower(region),'us') | strcmp(lower(region),'usa'))
			region = 'united states';
		end

% load data if needed

		if ~exist('POpatch','var'); load worldlo; end

% extract lat-long for countries. Error out if region unrecognized

		try
			[lat,lon,indx] = extractm(POpatch,region);
		catch
			linebreak = sprintf('\n');
			error(['Valid region strings are country names or' 	linebreak ...
					'   ''an'' or ''antarctica'','				linebreak ...
					'   ''a''  or ''asia'','					linebreak ...
					'   ''af'' or ''africa'','					linebreak ...
					'   ''pa'' or ''pacific'','					linebreak ...
					'   ''e''  or ''europe'','					linebreak ...
					'   ''na'' or ''north america'','			linebreak ...
					'   ''sa'' or ''south america'','			linebreak ...
					'   ''np'' or ''north pole'','			linebreak ...
					'   ''sp'' or ''south pole'','			linebreak ...
					'   ''af'' or ''africa''.' ])
		end
	
% make nice latitude and longitude limits from data

		inc = 5;
		buffer = 1;
		latlim = [ inc*floor(min(lat-buffer)/inc) inc*ceil(max(lat+buffer)/inc) ]; % Snap to 5 degree increments, with a 1 degree buffer
		lonlim = [ inc*floor(min(lon-buffer)/inc) inc*ceil(max(lon+buffer)/inc) ];

		if ~iscell(region) & ~isempty(strmatch('russia',lower(region))) % special case crossing the dateline
			lonlim = [ 5*floor(min(zero22pi(lon))/5) 5*ceil(max(zero22pi(lon))/5) ];
		end
		if ~iscell(region) & ~isempty(strmatch('united states',lower(region))) % special case crossing the dateline
			lonlim = [ 5*floor(min(zero22pi(lon))/5) 5*ceil(max(zero22pi(lon))/5) ]-360;
		end
	end

end

if lonlim(1) > lonlim(2)
	lonlim(2) = lonlim(2)+360; % ensure ascending lonlim
end

% Now have latlim, lonlim, in the range -180 + 200 (in cases computed from region)

% Compute a nice increment for labels and grid.
% Pick a ticksize that gives at least 3 grid lines

maxdiff = max(abs([diff(latlim) diff(lonlim)]));
mindiff = min(abs([diff(latlim) diff(lonlim)]));
ticcandidates = [.1 .5 1 2 2.5 5 10 15 20 30 45 ] ;
[ignored,indx] = min( abs( 3 - (mindiff ./ ticcandidates) ));

ticksize = ticcandidates(indx);
roundat = 0;
if mod(ticksize,1) ~= 0; roundat = -1; end

% Select a projection based on latlim,lonlim.

if isequal(latlim,[-90 90]) & isequal(lonlim,[-180 180]) % entire globe
	projection = 'robinson';
elseif max(abs(latlim)) < 20% straddles equator, but doesn't extend into extreme latitudes
	projection = 'mercator';
elseif abs(diff(latlim)) <= 90 & abs(sum(latlim)) > 20  & max(abs(latlim)) < 90 % doesn't extend to the pole, not stradling equator
	projection = 'eqdconic';
elseif abs(diff(latlim)) < 85  & max(abs(latlim)) < 90 % doesn't extend to the pole, not stradling equator
	projection = 'sinusoid';
elseif max(latlim) == 90 & min(latlim) >= 84
	projection = 'upsnorth';
elseif min(latlim) == -90 & max(latlim) <= -80
	projection = 'upssouth';
elseif max(abs(latlim)) == 90 & abs(diff(lonlim)) < 180
	projection = 'polycon';
elseif max(abs(latlim)) == 90 
	projection = 'eqdazim';
else
	projection = 'miller';
end


% Set up the axes with the selected projection and parameters.
% More than one path because of error messages setting parallels
% when projections don't have any.

if strcmp(projection,'upsnorth') | strcmp(projection,'upssouth')

	axesm('MapProjection','ups','Zone',projection(4:end),...
        	'Frame','on',...
    		'Grid','on',...
    		'LabelRotation','on',...
    		'MeridianLabel','on',...
    		'ParallelLabel','on',...
    		'MLabelParallel',0 ...
            )

else
    
    mstruct = defaultm(projection);
    
    if strcmp(projection,'eqdazim')
       axesm(...
    		'MapProjection',projection,...
    		'FLatLimit',[-Inf abs(diff(latlim))],...
    		'Origin',[90*sign(latlim(2)) mean(lonlim) 0],...
    		'MLabelLoc',ticksize,...
    		'MLineLoc',ticksize,...
    		'PLabelLoc',ticksize,...
    		'PLineLoc',ticksize,...
    		'MLabelRound',roundat,...
    		'PLabelRound',roundat,...
    		'GColor',[.75 .75 .75],... %[0.36471      0.72941      0.79216],... % fog
    		'GLineStyle',':',...
    		'Frame','on',...
    		'Grid','on',...
    		'LabelRotation','on',...
    		'MeridianLabel','on',...
    		'ParallelLabel','on',...
    		'MLabelParallel',0 ...
    		)
    
    elseif mstruct.nparallels > 0
    	axesm(...
    		'MapProjection',projection,...
    		'MapLatLimit',latlim,...
    		'MapLonLimit',lonlim,...
    		'MapParallels',[],...
    		'MLabelLoc',ticksize,...
    		'MLineLoc',ticksize,...
    		'PLabelLoc',ticksize,...
    		'PLineLoc',ticksize,...
    		'MLabelRound',roundat,...
    		'PLabelRound',roundat,...
    		'MLabelPar',0,...
    		'GColor',[.75 .75 .75],... %[0.36471      0.72941      0.79216],... % fog
    		'GLineStyle',':',...
    		'Frame','on',...
    		'Grid','on',...
    		'LabelRotation','on',...
    		'MeridianLabel','on',...
    		'ParallelLabel','on'...
    		)
    else
    	axesm(...
    		'MapProjection',projection,...
    		'MapLatLimit',latlim,...
    		'MapLonLimit',lonlim,...
    		'MLabelLoc',ticksize,...
    		'MLineLoc',ticksize,...
    		'PLabelLoc',ticksize,...
    		'PLineLoc',ticksize,...
    		'MLabelRound',roundat,...
    		'PLabelRound',roundat,...
    		'MLabelPar',0,...
    		'GColor',[.75 .75 .75],... %[0.36471      0.72941      0.79216],... % fog
    		'GLineStyle',':',...
    		'Frame','on',...
    		'Grid','on',...
    		'LabelRotation','on',...
    		'MeridianLabel','on',...
    		'ParallelLabel','on'...
    		)
    end
end

% Display basemap information

if onlyflag % just the region requested ( | iscell(region))
	switch lower(type)
	case {'line','lineonly'}
	
		if ~exist('POpatch','var'); load worldlo; end
	
		[lat,lon,indx] = extractm(POpatch,region);
		for i=1:length(indx)
			plotm(POpatch(indx(i)).lat,POpatch(indx(i)).long,'k','Tag',POpatch(indx(i)).tag)
		end
		
	case {'patch','patchonly'}
	
		if ~exist('POpatch','var'); load worldlo; end
	
		h = displaym(POpatch,region);
		zdatam(h,-1)

	end
else
	
	switch lower(type)
	case {'line','lineonly','mesh'}
		
		if exist('map','var') 
			meshm(map,maplegend)
		end
	

	
		if ~exist('POpatch','var'); load worldlo; end
	
		displaym(POline)
		displaym(DNline)
		
		for i=1:length(DNpatch) % convert large lake patches to lines
			DNpatch(i).type = 'line';
			DNpatch(i).otherproperty = {'color','k'};
		end
		displaym(DNpatch);
	
	case {'patch','patchonly'}
	
		if ~exist('POpatch','var'); load worldlo; end
	
		h = displaym(POpatch);
		zdatam(h,-1)
		h = displaym(DNpatch);
		set(h,'Facecolor','w')
		zdatam(h,-0.5)
		
		polcmap
	case 'meshonly'

		if exist('map','var')
			meshm(map,maplegend)
		end
	
	end

% Add populated places, if space permits

	switch lower(type)
	case {'line','patch','mesh'}
		
		if exist('map','var')
			meshm(map,maplegend)
		end
	
		if ~exist('POpatch','var'); load worldlo; end
	
		indx = find(([PPtext.lat]  > latlim(1)) & ([PPtext.lat]  < latlim(2)) & ...
					([PPtext.long] > lonlim(1)) & ([PPtext.long] < lonlim(2)) );

		if length(indx) < 15
						
			displaym(PPpoint(1)) % City point markers
				
			for i=indx
				PPtext(i).string = [ ' ' PPtext(i).string ' '];
			end
				
			h = displaym(PPtext(indx)); % City text
			set(h,'fontsize',9)
			rotatetext(h)
		end	
		
		indx = find(([POtext.lat]  > latlim(1)) & ([POtext.lat]  < latlim(2)) & ...
					([POtext.long] > lonlim(1)) & ([POtext.long] < lonlim(2)) );
		
		if length(indx) < 15
			h = displaym(POtext(indx));
			set(h,'fontsize',9)
			rotatetext(h)
		end
	end
end

tightmap



if nargout ==1;
	hout = gca;
end
