Code covered by the BSD License  

Highlights from
WORLDDATAMAP

image thumbnail

WORLDDATAMAP

by

 

28 Apr 2005 (Updated )

Map a region using data from a shapefile or data grid.

Editor's Notes:

This file was a File Exchange Pick of the Week

worlddatamap(varargin)
function hout = worlddatamap(varargin)
%WORLDDATAMAP Map a region using data from a shapefile or data grid
%
%   WORLDDATAMAP with no arguments presents a menu of region names from
%   which to choose, and then displays the boundaries of the selected area
%   using vector data from the default shapefile 'landareas.shp'. The
%   boundaries are symbolized as lines.
%
%   WORLDDATAMAP(REGION) makes a map of the specified region using vector
%   data from the default shapefile, 'landareas.shp'. REGION is a string or
%   a cell array of strings. The projection and limits are suitable to the
%   part of the world specified by REGION. Permissible strings include
%   names of continents and countries as well as 'World', 'North Pole',
%   'South Pole', and 'Pacific'.  The region boundaries are symbolized as
%   lines, provided vector data from the default shapefile exist within the
%   region boundary.
%
%   WORLDDATAMAP(REGIONONLY) adds only the specified region to the map. If
%   any of the names in the REGIONONLY string or cell array of strings end
%   with 'only', only the requested regions are displayed; therfore, only
%   one string needs the 'only' suffix.
%
%   WORLDDATAMAP(LATLIM, LONLIM) makes a map of data within a geographic
%   region specified as latitude and longitude limits in degrees. LATLIM
%   and LONLIM are two-element vectors of the form [southern_limit
%   northern_limit] and [western_limit eastern_limit], respectively.
%   Latitudes range from [-90 90] and longitudes range from [-180 180] for
%   the default shapefile, 'landareas.shp'.
%
%   WORLDDATAMAP(...,VECDISPLAYTYPE) where VECDISPLAYTYPE is a string,
%   controls the rendering of the vector data.   If omitted, VECDISPLAYTYPE
%   defaults to 'line' Permissable values for VECDISPLAYTYPE are:
%
%    Value         Type of render
%    --------      --------------
%    'line'        line  (default)
%    'patch'       patch 
%
%   WORLDATADMAP(Z, REFVEC) derives the map limits from the extent of a
%   regular data grid Z, with a 1-by-3 referencing vector REFVEC. The
%   regular data grid is displayed using MESHM, and the region boundaries
%   are sybmolized as lines, from the default shapefile, 'landareas.shp'.
%
%   WORLDDATAMAP(Z, REFVEC, SURFDISPLAYTYPE) where SURFDISPLAYTYPE is a
%   string, controls the rendering of the surface data.   If omitted,
%   SURFDISPLAYTYPE defaults to 'mesh'. Permissable values for
%   SURFDISPLAYTYPE are:
%
%    Value         Type of render
%    --------      -----------------------------------------
%    'mesh'        surface in z=0 plane using MESHM (default)
%    'mesh3d'      surface, set aspect ratio
%    'dem'         surface, set digital elevation colormap
%    'dem3d'       surface, set digital elevation colormap, aspect ratio
%    'lmesh3d'     surface, apply lighting
%    'ldem3d'      surface, apply lighting, digital elevation colormap
%    'line'        only line boundaries defined by extent of data grid
%    'patch'       only patch boundaries defined by extent of data grid
%
%   WORLDDATAMAP(...,FEATURETYPE) where FEATURETYPE is a string or cell
%   array of strings, controls the display of additional feature data. The
%   features are obtained from the shapefiles, worldcities.shp,
%   worldrivers.shp, or worldlakes.shp. If omitted, FEATURETYPE defaults to
%   'none'. Multiple features are allowed if FEATURETYPE is a cell array of
%   strings. Permissable values for FEATURETYPE are:
%
%    Value         Type of feature
%    --------      --------------
%    'none'        no additional features (default)
%    'cities'      draw all World cities from shapefile
%    'rivers'      draw all World rivers from shapefile
%    'lakes'       draw all World lakes from shapefile
%    'labels'      draw labels (if LabelLat, LabelLon defined in shapefile)
%
%   WORLDDATAMAP(...,DISPLAYTYPEONLY, FEATURETYPE) where DISPLAYTYPEONLY is
%   either VECDISPLAYTYPE ('lineonly' or 'patchonly') or SURFDISPLAYTYPE 
%   ('meshonly', 'mesh3donly', etc) with the suffix 'only', prohibits the
%   features defined by FEATURETYPE ('cities', 'rivers', etc) from being
%   rendered.
%
%   WORLDDATAMAP(...,DISPLAYTYPE, FEATURETYPEONLY) where FEATURETYPEONLY is
%   a FEATURETYPE with the suffix 'only' ('citiesonly', 'riversonly', etc), 
%   and DISPLAYTYPE is either a VECDISPLAYTYPE ('line', 'patch') or
%   SURFDISPLAYTYPE ('mesh', 'dem', etc), symbolizes only the requested
%   feature.
%
%   WORLDDATAMAP(FILENAME, ...) where FILENAME is a string ending in .shp
%   or .dbf, reads the vector data from shapefile FILENAME instead of the
%   default file, landareas.shp.
%
%   WORLDDATAMAP({FILENAME, REGIONATTRIBUTE}, ...) uses the REGIONATTRIBUTE
%   field name rather than the default, 'Name' to search for a specified
%   REGION. If REGIONATTRIBUTE is missing, it defaults to the first match
%   to the following case-insensitive list:
%   'name','cntry_name','country_name','long_name'.
%
%   h = WORLDDATAMAP(...) returns the handle of the map axes.
%
%   WORLDDATAMAP uses TIGHTMAP to set the axis limits tight around the map.
%   If you change the projection, or need more white space around the map
%   frame, use TIGHTMAP again or AXIS AUTO.
%
%   Example 1
%   ---------
%   % World map with coastlines
%   figure
%   worlddatamap('World')
%
%   Example 2
%   ---------
%   % Worldmap with land areas, major lakes and rivers, and cities and
%   % populated places
%   figure
%   worlddatamap('World',{'lakes','rivers','cities'})
%  
%   Example 3
%   ---------
%   % Patch map of Antarctica
%   figure
%   worlddatamap('antarctica','patch')
%
%   Example 4
%   ---------
%   % Map of Africa and India with major cities and populated places
%   figure
%   worlddatamap({'Africa','India'},'cities')
%  
%   Example 5
%   ---------
%   % Map of the geoid over South America and the central Pacific
%   figure
%   worlddatamap([-50 50],[160 -30])
%   load geoid
%   geoshow(geoid, geoidrefvec, 'DisplayType', 'texturemap');
%   load coast
%   geoshow(lat, long)
%
%   Example 6
%   ---------
%   % Map of terrain elevations in Korea
%   figure
%   load korea
%   worlddatamap(map, refvec,'dem')
%  
%   Example 7
%   ---------
%   % Map of the United States of America
%   figure
%   worlddatamap USA
%
%   Remarks:  
%   WORLDDATAMAP requires MATLAB version 7.0.4 (R14SP2) (or higher) and the
%   Mapping Toolbox version 2.1 (or higher).
%
%   See also MESHM, SHAPEREAD, USAMAP, WORLDMAP.

%   Copyright 2005 The MathWorks, Inc.
%   Author: Kelly G. Luetkemeyer

% Define the external data API
api = defineExternalDataApi;

% Verify Mapping Toolbox version for ancillary data
verifyVersion(api);

% Parse and check input arguments
[worldmapArgs, dataset, attribName, displayType, featureType, latlim, lonlim, ...
   region, Z, refvec, regionOnly, listSelect] = parseInputs(varargin{:});

% Verify the dataset
[worldShapeFile, attribName, usingDefaults] = verifyDataSet(api, dataset, attribName);
api.ExternalData(1).Filename = worldShapeFile;

% Change region to lat-lon limit if not using default attribute name.
if ~usingDefaults && ~isempty(region)
  worldmapArgs = getLatLonLimitsFromFile(dataset, attribName, region, worldmapArgs);
end

% Display the base worldmap
if ~listSelect
   h = worldmap(worldmapArgs{:});
else
   % Select the region from the list dialog
   try
      h = worldmap(worldmapArgs{:});
   catch
      % Cancel was clicked.
      close
      return
   end
end
tightmap loose

% Obtain latlim and lonlim from the projection if not specified
if ~isempty('latlim') || ~isempty('lonlim')
   m = gcm;
   latlim = m.maplatlimit;
   lonlim = m.maplonlimit;
end

% Change [0 360] longtitude to [-180 180]
latlim = latlim(:)';
lonlim = lonlim(:)';
if any(lonlim>180)
   lonlim = npi2pi(lonlim);
end
if lonlim == [0 0]
   lonlim = [-180 180];
end

% Display the grid data if specified
if ~isempty(Z) && ~isempty(refvec)
   addGridData(Z, refvec, displayType);
end

% Display the vector data
addWorldVectorData(api, region, regionOnly, latlim, lonlim, ...
   displayType, featureType, worldShapeFile, attribName);

% Return handle if required
if nargout == 1
   hout = h;
end

%--------------------------------------------------------------------------
function api = defineExternalDataApi
% Assign required version (or greater) and data filesnames 
api.requiredVersion   = '2.1';
api.defShapeName      = 'landareas';
api.defAttribName     = 'Name';

api.ExternalData(1).Filename = 'landareas.shp';
api.ExternalData(1).DrawFcn  = '';
api.ExternalData(1).ValidTypes = {};

api.ExternalData(2).Filename = 'worldcities.shp';
api.ExternalData(2).DrawFcn  = @drawCities;
api.ExternalData(2).ValidTypes = {'cities','citiesonly', 'all'};

api.ExternalData(3).Filename = 'worldlakes.shp';
api.ExternalData(3).DrawFcn  = @drawLakes;
api.ExternalData(3).ValidTypes = {'lakes','lakesonly', 'all'};

api.ExternalData(4).Filename = 'worldrivers.shp';
api.ExternalData(4).DrawFcn  = @drawRivers;
api.ExternalData(4).ValidTypes = {'rivers','riversonly', 'all'};

api.ExternalData(5).Filename = 'DEFAULTFILE';
api.ExternalData(5).DrawFcn  = @drawLabels;
api.ExternalData(5).ValidTypes = {'labels','labelsonly', 'all'};

api.mappingShapeFiles = {api.ExternalData(1:4).Filename};
%--------------------------------------------------------------------------
function [worldmapArgs, dataset, attribName, displayType, featureType, latlim, lonlim, ...
  region, Z, refvec, regionOnly, listSelect] = parseInputs(varargin)

% Default values for inputs
dataset = '';
attribName = '';
latlim = [];
lonlim = [];
region = [];
Z      = [];
refvec = [];
regionOnly = false;

if numel(varargin) == 0
   worldmapArgs = varargin;
   displayType = 'line';
   featureType = {'none'};
   listSelect = true;
   return
else
   listSelect = false;
end

% Determine if resolution 'hi', 'lo', has been set.
[worldmapArgs, resolution] = findResolutionArg(varargin);

% Determine if a dataset has been set.
[worldmapArgs, dataset, attribName] = findDataSetArgs(worldmapArgs);

% Determine if type arguments have been set.
[worldmapArgs, displayType, featureType] = findTypeStrings(worldmapArgs);

switch numel(worldmapArgs)
   case 0   % WOLRDDATAMAPA
      listSelect = true;

   case 1   % WORLDDATAMAP(REGION)
      region = worldmapArgs{1};

      %ensure region is a cell array of strings
      if ischar(region)
           region = cellstr(region);
      elseif ~iscell(region)
         error('%s\n\n%s\n\n%s%s%s', ...
            'Function WORLDDATAMAP expected REGION to be one of these types:', ...
            'string or string cell array', ...
            'Instead its type was ', class(region),'.');
      end

      % Parse only from region, if present
      if ~isempty(region)
         [region,regionOnly]=removeOnlyFromRegion(region);
         if regionOnly
            % Remove warning from WORLDMAP for only flag
            %  by re-setting the region name
            worldmapArgs{1} = region;
         end
      else
         region = {'world'};
      end

   case 2   % WORLDDATAMAP(Z, REFVEC)
            % WORLDDATAMAP(LATLIM, LONLIM)

      if numel(worldmapArgs{1}) ~= 2 && numel(worldmapArgs{2}) ~= 2
         % WORLDDATAMAP(Z, REFVEC)
         Z = worldmapArgs{1};
         refvec = worldmapArgs{2};
         if isempty(displayType)
            displayType = 'mesh';
         end
      else
         % WORLDDATAMAP(LATLIM, LONLIM)
         latlim = worldmapArgs{1};
         lonlim = worldmapArgs{2};
      end

   otherwise

      if isempty(displayType) && ischar(worldmapArgs{end})
         typeStrings = getDisplayTypeStrings;
         error('%s\n\n%s\n%s\n','Invalid DisplayType string.', ...
                                'Valid DisplayType strings are:',typeStrings{:});
      end
      error('Incorrect number of input arguments.')
end

% Set default displayType if not specified
if isempty(displayType)
   displayType = 'line';
end

% Add resolution if present
if ~isempty(resolution)
   worldmapArgs = {resolution worldmapArgs{:}};
end


%--------------------------------------------------------------------------
function typeStrings = getDisplayTypeStrings
typeStrings = {'patch','line','patchonly','lineonly','mesh',...
   'meshonly','mesh3d', 'mesh3donly', 'dem','demonly','dem3d','dem3donly',...
   'ldem3d','ldem3donly','lmesh3d','lmesh3donly'};

%--------------------------------------------------------------------------
function typeStrings = getFeatureTypeStrings
typeStrings = { 'cities','citiesonly', ...
   'lakes','lakesonly','rivers','riversonly', ...
   'labels', 'labelsonly', ...
   'all', 'none'};

%--------------------------------------------------------------------------
function [worldmapArgs, resolution] = findResolutionArg(worldmapArgs)
% Find the resolution arg if present.

resolutionArgs = {'hi', 'lo', 'allhi'};
resolution = worldmapArgs{1};

firstArgMatches = ischar(resolution) && any(strcmpi(resolution,resolutionArgs)); 

if firstArgMatches
   worldmapArgs(1) = [];
else
   resolution = '';
end

%--------------------------------------------------------------------------
function [worldmapArgs, dataset, attribName] = findDataSetArgs(worldmapArgs)
% Find the dataset arg if present.

attribName = '';
if iscell(worldmapArgs{1})
   dataset = worldmapArgs{1}{1};
   if numel(worldmapArgs{1}) > 1
      attribName = worldmapArgs{1}{2};
   end
else
   dataset = worldmapArgs{1};
end

if containsShapeExt(dataset);
   worldmapArgs(1) = [];
   if ~isempty(attribName)  && ~ischar(attribName)
       error('%s\n%s%s%s', ...
            'Function WORLDDATAMAP expected REGIONATTRIBUTE to be type char.', ...
            'Instead its type was ', class(attribName),'.');
   end
else
   dataset = '';
   attribName = '';
end

%--------------------------------------------------------------------------
function shapeFilePresent = containsShapeExt(filename)
if ~ischar(filename)
   shapeFilePresent = false;
else
   filename = lower(filename);
   shapeFilePresent = ~isempty( findstr('.shp', filename ) ) || ...
                      ~isempty( findstr('.dbf', filename ) ) ;
end

%--------------------------------------------------------------------------
function verifyVersion(api)
requiredVersion   = api.requiredVersion;
mappingShapeFiles = api.mappingShapeFiles;

for j=1:numel(mappingShapeFiles)
   if ~exist(mappingShapeFiles{j},'file')
      msg1 = sprintf('%s%s\n','Unable to find shapefile: ',mappingShapeFiles{j});
      if ~isdeployed
         mappingVersion = ver('map');
         if isempty(mappingVersion) || isempty(mappingVersion.Version)
            installedVersion = '0.0.0';
            versionStatement = 'The Mapping Toolbox is not installed.';
         else
            installedVersion = mappingVersion.Version;
            versionStatement = ['The installed version of the Mapping Toolbox is ', ...
               installedVersion]; ...
         end

         if (  installedVersion(1) <  requiredVersion(1) ) || ...
            (( installedVersion(1) == requiredVersion(1) ) && ...
             ( installedVersion(3) <  requiredVersion(3) ) )

            msg2 = sprintf('%s%s%s%s%s\n%s\n%s', ...
               'Function ',upper(mfilename), ...
               ' requires Mapping Toolbox version ',requiredVersion, ' or greater.', ...
               versionStatement, ...
               'Please install the latest version of the Mapping Toolbox.');
            error([msg1 msg2]);
         else
            error(msg1);
         end

      else
         msg2 = sprintf('%s%s\n%s%s%s',...
            'The application is missing shapefile:',mappingShapeFiles{j}, ...
            'Compile the application using the -a ', ...
            mappingShapeFiles{j},' option to include the file.');
         error([msg1 msg2]);
      end
   end
end

%--------------------------------------------------------------------------
function [worldShapeFile, attribName, usingDefaults] = verifyDataSet(...
   api, dataset, attribName)

defAttribName     = api.defAttribName;
defShapeFile      = api.defShapeName;

% Assign default shapefile
defaultDataSetFile = [defShapeFile '.shp'];

% Assign default cases if empty
if isempty(dataset) 
   dataSetFile = defaultDataSetFile;
   attribName = defAttribName;
   usingDefaults = true;
else

   dataSetFile = dataset;
   if ~exist(dataSetFile, 'file')
      [path, name] = fileparts(dataSetFile);
      defaultDir = fullfile(fileparts(which(mfilename)), name);
      if ~exist(fullfile(defaultDir,dataSetFile),'file');
         msg = sprintf('%s%s','Failed to locate SHAPEFILE: ',dataSetFile);
         error(msg);
      else
         fprintf('%s%s%s\n','Adding directory ''', defaultDir, ''' to the path.');
         addpath(defaultDir);
      end
   end

   usingDefaults = false;
   if isempty(attribName)
      allowableAttribNames = {'name','cntry_name','country_name','long_name'};
   else
      allowableAttribNames = {attribName};
   end
   info = shapeinfo(dataSetFile);
   for j=1:numel(allowableAttribNames)
      containsName = strcmpi(allowableAttribNames{j}, {info.Attributes.Name});
      if any(containsName)
         attribName = info.Attributes(containsName).Name;
         break;
      end
   end
end
if isempty(attribName)
  msg = sprintf('%s%s%s','REGIONATTRIBUTE for file ''',dataSetFile,''' is empty.');
  error(msg);
end
worldShapeFile = dataSetFile;

%--------------------------------------------------------------------------
function worldmapArgs = getLatLonLimitsFromFile(dataset, attribName, region, worldmapArgs)
% Find the bounding box of the region
s=[];
for k=1:numel(region)
    try
        
    shape = shaperead(dataset, 'UseGeoCoords', true, 'Selector',...
        {@(name) strcmpi(name,region{k}), attribName});
    s = [s, shape];
    catch
        % If any region error's, 
        % then the region is not found in the file,
        % but might be found by worldmap
        return
    end
end
switch numel(s)
    case 0     % s is empty
        return % Region might be found from worldmap
        
    case 1
        bbox = [s.BoundingBox];
        
    otherwise
        bbox = s(1).BoundingBox;
        for num = 2:numel(s)
            bbox = [ min(bbox(1,1)',[s(num).BoundingBox(1,1)]'), ...
            min(bbox(1,2)', [s(num).BoundingBox(1,2)]'); ...
            max(bbox(2,1)', [s(num).BoundingBox(2,1)]'), ...
            max(bbox(2,2)', [s(num).BoundingBox(2,2)]') ];
        end
       
end

lon = bbox(:,1:2:end);
lat = bbox(:,2:2:end);
lonlim = [min(lon(:)) max(lon(:))];
latlim = [min(lat(:)) max(lat(:))];

% Snap map limits to increments of INC, with a 1 degree buffer, except
% for limits that are already exact multiples of INC.
buffer = 1;
if strmatch('usa', dataset)
    inc = 1;  % 1 degree for usa states
else
    inc = 5;  % compatible with worldmap (5 degree border)
end

if mod(latlim(1),inc) ~= 0
    latlim(1) = inc * floor((latlim(1) - buffer)/inc);
end
if mod(lonlim(1),inc) ~= 0
    lonlim(1) = inc * floor((lonlim(1) - buffer)/inc);
end
if mod(latlim(2),inc) ~= 0
    latlim(2) = inc * ceil((latlim(2) + buffer)/inc);
end
if mod(lonlim(2),inc) ~= 0
    lonlim(2) = inc * ceil((lonlim(2) + buffer)/inc);
end

% Ensure that latitude limits remain within [-90 90].
if latlim(1) < -90
    latlim(1) = -90;
end
if latlim(2) > 90
    latlim(2) = 90;
end

if latlim(1) <= -80 && latlim(2) >=80
    % changes projection
    latlim = [-90 90];
end
worldmapArgs(1) = [];
worldmapArgs = {latlim, lonlim, worldmapArgs{:}};

%--------------------------------------------------------------------------
function [worldmapArgs, displayType, featureType] = findTypeStrings(worldmapArgs)
% Find the typeString if present


switch numel(worldmapArgs)
   case 0
      displayType = '';
      featureType = {'none'};
   case 1
     [displayType, featureType, worldmapArgs] = findDisplayFeatureTypeAtEnd(worldmapArgs);

   otherwise
     % ..., DisplayType, FeatureType
     % ..., FeatureType
     % ..., DisplayType
    
     [displayType, featureType, worldmapArgs] = findDisplayFeatureTypeAtEnd(worldmapArgs);
     if isempty(displayType) && isempty(featureType)
       displayType = '';
       featureType = {'none'};
     elseif isempty(displayType) && ischar(worldmapArgs{end}) 
       displayTypeStrings = getDisplayTypeStrings;
       argMatchesDisplayTypeString = any(strcmpi(worldmapArgs{end},displayTypeStrings));
       if argMatchesDisplayTypeString
         displayType = lower(worldmapArgs{end});
         worldmapArgs(end) = [];
       end
    end
end
if ~iscell(featureType)
  featureType = {featureType};
end

%--------------------------------------------------------------------------
function [displayType, featureType, worldmapArgs] = findDisplayFeatureTypeAtEnd(worldmapArgs)
displayType = '';
featureType = '';
argMatchesDisplayTypeString = false;
argMatchesFeatureTypeString = false;

if ischar(worldmapArgs{end}) 
   displayTypeStrings = getDisplayTypeStrings;
   featureTypeStrings = getFeatureTypeStrings;
   argMatchesDisplayTypeString = any(strcmpi(worldmapArgs{end},displayTypeStrings));
   argMatchesFeatureTypeString = any(strcmpi(worldmapArgs{end},featureTypeStrings));
   
elseif iscell(worldmapArgs{end})
   featureTypeStrings = getFeatureTypeStrings;
   type = worldmapArgs{end};
   matchesFeatureTypeString(1:numel(type)) = false;
   for k=1:numel(type)
      matchesFeatureTypeString(k) = any(strcmpi(type{k},featureTypeStrings));
   end
   argMatchesFeatureTypeString = all(matchesFeatureTypeString);
end

if argMatchesDisplayTypeString
   displayType = lower(worldmapArgs{end});
   worldmapArgs(end) = [];
elseif argMatchesFeatureTypeString
   featureType = lower(worldmapArgs{end});
   worldmapArgs(end) = [];
end

%--------------------------------------------------------------------------
function [region,regionOnly] = removeOnlyFromRegion(region)

% detect the 'only case',
regionOnly = false;

%only case with cell array of strings
for i=1:length(region)
   thisregion = region{i};
   if length(thisregion) >= 5 && strcmp(thisregion(end-4:end),' only')
      thisregion = thisregion(1:end-5);
      regionOnly = true;
   elseif length(thisregion) >= 4 && strcmp(thisregion(end-3:end),'only')
      thisregion = thisregion(1:end-4);
      regionOnly = true;
   end
   region{i} = thisregion;
end

%--------------------------------------------------------------------------
function s = loadshape(shapefile, attribName, varargin)
% LOADSHAPE(SHAPEFILE, ATTRIBNAME)
% LOADSHAPE(SHAPEFILE, ATTRIBNAME, BBOX)

args = {shapefile, 'UseGeoCoords', true};
switch numel(varargin)
   case 0
   case 1
      bbox = varargin{1};
      args = {args{:}, 'BoundingBox', bbox};
   otherwise
      bbox = varargin{1};
      args = {args{:}, 'BoundingBox', bbox};
end

if ~isempty(attribName)
   args = {args{:}, 'Attributes', {attribName}};
end

s = shaperead(args{:});

%--------------------------------------------------------------------------
function addGridData(Z, refvec, displayType)

switch displayType
   case {'mesh','dem','meshonly','demonly'}
      meshm(Z,refvec)

   case {'dem3d','mesh3d','ldem3d','lmesh3d',...
         'dem3donly','mesh3donly','ldem3donly','lmesh3donly'}
      meshm(Z,refvec,size(Z),Z)

   otherwise

end


% Control aspect ratio for 3d plots
switch displayType
   case {'dem3d','ldem3d'}
      da = daspect;
      pba = pbaspect;
      da(3) = 7.5*pba(3)/da(3);
      daspect(da);

   case {'mesh3d','lmesh3d'}
      da = daspect;
      pba = pbaspect;
      da(3) = 2*pba(3)/da(3);
      daspect(da);

   otherwise

end

% Add digital elevation colormap
switch displayType
   case {'dem','dem3d','ldem3d','demonly','dem3donly','ldem3donly'}
      demcmap(Z)

   otherwise
end

% Add lighting effects
switch displayType
   case {'lmesh3d','ldem3d'}
      camlight(90,5);
      camlight(0,5);
      lighting phong
      material([0.25 0.8 0])

   otherwise

end

%--------------------------------------------------------------------------
function addWorldVectorData(api, region, regionOnly, latlim, lonlim, ...
   displayType, featureType, worldShapeFile, attribName)


if regionOnly 

   % REGIONONLY is requested with 'only'
   s = readRegionOnly(worldShapeFile, attribName);

   if isempty(s)
     return
   end

   
   if isempty(regexp(featureType{1}, 'only', 'once'))
      switch displayType

        case {'line','lineonly'}
            drawRegionOnlyAsLines(s, region, attribName)

         case {'patch','patchonly'}
            drawRegionOnlyAsPatches(s, region, attribName);

      end
   end

   if ~isequal('none',featureType{1}) && isempty(regexp(displayType, 'only', 'once'))

         for j=1:numel(region)
            indx = find(strcmpi(region{j},{s.(attribName)}));

            if ~isempty(indx) && numel(indx) == 1
               bbox = s(indx).BoundingBox;
            else
               bbox = [ min(lonlim), min(latlim); max(lonlim), max(latlim)];
            end

            drawExternalData(api.ExternalData, featureType, attribName, bbox);
             
         end
   end

else

   % Full region is specified
   [s, bbox] = readRegion(worldShapeFile, latlim, lonlim, attribName);
   if isempty(s)
      return
   end

      
   if isempty(regexp(featureType{1}, 'only', 'once'))
      switch displayType

         case {'lineonly', 'line', 'mesh', 'dem', 'ldem'}
            drawRegionAsLines(s, attribName)


         case {'patch','patchonly'}
            drawRegionAsPatches(s, attribName)

      end
   end
   
  
   if ~isequal('none',featureType{1}) && ...
         isempty(regexp(displayType, 'only', 'once'))

      % Add the external data
       drawExternalData(api.ExternalData, featureType, attribName, bbox);
   end
   
end

%--------------------------------------------------------------------------
function s = readRegionOnly(filename, attribName)
s = loadshape(filename, attribName);
if ~isempty(s) && numel(s(end).Lat) < 3
   % Remove degenerative point 
   s(end) = [];
end

%--------------------------------------------------------------------------
function [s, bbox] = readRegion(filename, latlim, lonlim, attribName)

bbox = [ min(lonlim), min(latlim); max(lonlim), max(latlim)];
if (bbox == [lonlim', latlim'])
   s = loadshape(filename,attribName,bbox);
else
   % May have longitude wrapping, read the whole file
   s = loadshape(filename, attribName);
end
if ~isempty(s) && numel(s(end).Lat) < 3
   % Remove degenerative point 
   s(end) = [];
end

%--------------------------------------------------------------------------
function drawRegionOnlyAsLines(s, region, attribName)

for j=1:numel(region)
   indx = find(strcmpi(region{j},{s.(attribName)}));
   if ~isempty(indx)
     lat = [s(indx).Lat];
     lon = [s(indx).Lon];
     tag = s(indx(1)).(attribName);
   else
     lat = [s.Lat];
     lon = [s.Lon];
     tag = region{j};
   end
   plotm(lat, lon, 'k', 'Tag', tag);
end

%--------------------------------------------------------------------------
function drawRegionAsLines(s, attribName)
for i=1:numel(s)
   plotm(s(i).Lat,s(i).Lon,'k','Tag',s(i).(attribName))
end

%--------------------------------------------------------------------------
function drawRegionOnlyAsPatches(s, region, attribName)
cmap = polcmap(numel(region));
for k=1:numel(region)
   indx = find(strcmpi(region{k},{s.(attribName)}));
   if ~isempty(indx)
     lat = [s(indx).Lat];
     lon = [s(indx).Lon];
     tag = s(indx(1)).(attribName);
     h=patchm(lat, lon, cmap(k,:));
     set(h, 'Tag', tag);
   else
     for j=1:numel(s)
         h = patchesm(s(j).Lat, s(j).Lon, 0.0, ...
                      'Cdata', k, 'FaceColor', 'flat');
         set(h,'Tag',region{k});
     end
   end
end

%--------------------------------------------------------------------------
function drawRegionAsPatches(s, attribName)
alt = 0;
cmap = polcmap(numel(s));
for j=1:numel(s)
    h = patchesm(s(j).Lat, s(j).Lon, alt, ...
                'Cdata', j, 'FaceColor', 'flat');
   set(h, 'Tag', s(j).(attribName));
end
set(gcf,'colormap',cmap);

%--------------------------------------------------------------------------
function drawExternalData(ExternalData, featureType, attribName, bbox)
for j=1:numel(ExternalData)
   for k = 1:numel(featureType)
      if any(strcmpi(featureType{k}, ExternalData(j).ValidTypes))
          if isequal(ExternalData(j).Filename, 'DEFAULTFILE')
             filename = ExternalData(1).Filename;
          else
             filename = ExternalData(j).Filename;
          end
          ExternalData(j).DrawFcn(filename, attribName, bbox);
      end
   end
end

%--------------------------------------------------------------------------
function drawLabels(filename, attribName, bbox)
info = shapeinfo(filename);
if any(strcmpi('LabelLat', {info.Attributes.Name})) && ...
   any(strcmpi('LabelLon', {info.Attributes.Name}))

   m = gcm;
   s = shaperead(filename, 'UseGeoCoords', true, ...
      'Attributes', {'LabelLat','LabelLon', attribName});
  xv = [m.maplonlimit(1), m.maplonlimit(1), m.maplonlimit(2), m.maplonlimit(2)];
  yv = [m.maplatlimit(1), m.maplatlimit(2), m.maplatlimit(2), m.maplatlimit(1)];
   for k=1:numel(s)
       if inpolygon(s(k).LabelLon, s(k).LabelLat, xv, yv)
           h=textm(s(k).LabelLat, s(k).LabelLon, s(k).(attribName), ...
           'fontsize',9, 'fontweight','bold','Tag', s(k).(attribName));
       end
   end
end
  
%--------------------------------------------------------------------------
function drawRivers(shapefile, attribName, bbox)

rivers = loadshape(shapefile, 'Name', bbox);
for i=1:length(rivers)
   linem(rivers(i).Lat, rivers(i).Lon, 'b', 'Tag', rivers(i).Name);
end

%--------------------------------------------------------------------------
function drawLakes(shapefile, attribName, bbox)

lakes = loadshape(shapefile, 'Name', bbox);
for i=1:length(lakes)
   h=patchm(lakes(i).Lat, lakes(i).Lon, 'b');
   set(h,'Tag', lakes(i).Name);
end


%--------------------------------------------------------------------------
function drawCities(shapefile, attribName, bbox)

cities = loadshape(shapefile, 'Name', bbox);
for i=1:length(cities)
   linem(cities(i).Lat,cities(i).Lon,0, '.', ...
         'color' , [0.9000 0.6000 0.1000],  ...
         'markersize', 10, 'Tag', cities(i).Name)
   textm(cities(i).Lat, cities(i).Lon, [' ' cities(i).Name], ...
         'fontsize',9,     'Tag', cities(i).Name)
end

Contact us