Understanding Vector Geodata :: Understanding Map Data (Mapping Toolbox™)

Mapping Toolbox™

Types of Map Data Handled by the Toolbox

Understanding Vector Geodata

On this page…
Points, Lines, Polygons Segments Versus Polygons Mapping Toolbox™ Geographic Data Structures Selecting Data to Read with the shaperead Function

On this page…

Points, Lines, Polygons

Segments Versus Polygons

Mapping Toolbox™ Geographic Data Structures

Selecting Data to Read with the shaperead Function

Points, Lines, Polygons

Vector geospatial data is used to represent linear features such as rivers, coastlines, boundaries, and highways. Vector data can also represent areal features such as water bodies, political units, and enumeration districts. This section familiarizes you with how vector data structures digitally encode geographic entities and how to use this form of data.

In the context of geodata, vector data means "geometric descriptions of geographic objects" rather than its more general mathematical definition, "a quantity specified by a magnitude and a direction." In fact, some vector geodata is specified as points having neither magnitude nor direction. Other geodata—such as postcodes, highway mileposts, or census statistics—only implies an underlying geometry, which vector 2-D coordinate data is required to map or spatially analyze.

In the MATLAB^® workspace, vector data is expressed as pairs of variables that represent the geographic or plane coordinates for a set of points of interest. For example, the following two variables can be mapped as a vector:

lat = [45.6 -23.47 78];
long = [13 -97.45 165];

Note that either row or column vectors can be used, but both variables should have the same shape. For example, lat and long could be defined as columns:

lat = [45.6 -23.47 78]';
long = [13 -97.45 165]';

These values could mean anything. They could represent three locations over which geosynchronous satellites are stationed, and can be communicated by plotting a symbol for each point on a map of the Earth. Alternatively, they might represent a starting point, a midcourse marker, and a finish point of a sailboat race, in which case they can be rendered by plotting two line segments. Or perhaps the values represent the vertices of a triangle bounding a region of interest, and thus constitute a simple polygon.

Note When polygons become graphic objects, they are called patches. In this documentation, the words patch and polygon are often used interchangeably.

Mapping Toolbox™ functions provide for each of these interpretations. For many purposes, the distinction is irrelevant; for others, the choice of a function implies one interpretation over the others. For example, the function plotm displays the data as a line, while fillm displays it as a filled polygon. While you can draw an unfilled polygon with fillm that looks like the output from plotm, the resulting object has a different graphic data type (patch versus line), hence different properties you can set.

A line must contain at least two coordinate elements for each coordinate dimension, and a polygon at least three (note that it is not necessary to duplicate the first point as the last point to define or render a polygon). The toolbox places no limit (beyond available memory) on how large or how complex the shape of a line and polygon can be, other than the restriction that it should not cross itself.

Objects in the real world that vector geodata represents can have many parts, for example, the islands that make up the state of Hawaii. When encoding as vector variables the shapes of such compound entities, you must separate successive entities. To indicate that such a discontinuity exists, the toolbox uses the convention of placing NaNs in identical positions in both vector variables. For example, if a second segment is to be added to the preceding map, the two objects can reside in the same pair of variables:

lat = [45.6 -23.47 78 NaN 43.9 -67.14 90 -89];
lon = [13 -97.45 165 NaN 0 -114.2 -18 0];

Notice that the NaNs must appear in the same locations in both variables. Here is a segment of three points separated from a segment of four points. The NaNs perform two functions: they provide a means of identifying breakpoints in the data, and they serve as pen-up commands when plotting vector maps. The NaNs are used to separate both distinct (but possibly connecting) lines and disconnected patch faces.

Note This convention departs from regular MATLAB graphics, in which NaN-separated polygons cannot be interpreted or displayed as separate patches.

Segments Versus Polygons

Geographic objects represented by vector data might or might not be formatted as polygons. Imagine two variables, latcoast and loncoast, that correspond to a sequence of points that caricature the coast of the island of Great Britain. If this data returns to its starting point, then a polygon describing Great Britain exists. This data might be plotted as a patch or as a line, and it might be logically employed in calculations as either.

Now suppose you want to represent the Anglo-Scottish border, proceeding from the west coast at Solway Firth to the east coast at Berwick-upon-Tweed. This data can only be properly defined as a line, defined by two or more points, which you can represent with two more variables, latborder and lonborder. When plotted together, the two pairs of variables can form a map. The patch of Great Britain plus the line showing the Scottish border might look like two patches or regions, but there is no object that represents England and no object that represents Scotland, either in the workspace or on the map axes.

In order to represent both regions properly, the Great Britain polygon needs to be split at the two points where the border meets it, and a copy of latborder and lonborder concatenated to both lines (placing one in reverse order). The resulting two polygons can be represented separately (e.g., in four variables named latengland, lonengland, latscotland, and lonscotland) or in two variables that define two polygons each, delineated by NaNs (e.g., latuk, lonuk).

The distinction between line and polygon data might not appear to be important, but it can make a difference when you are performing geographic analysis and thematic mapping. For example, polygon data can be treated as line data and displayed with functions such as linem, but line data cannot be handled as polygons unless it is restructured to make all objects close on themselves, as described in Matching Line Segments.

Mapping Toolbox™ Geographic Data Structures

In examples provided in prior chapters, geodata was in the form of individual variables. Mapping Toolbox software also provides an easy means of displaying, extracting, and manipulating collections vector map features that have been organized in a family of specially organized geographic data structures.

Note The data structures discussed in this section are different from the map projection structure (also called an mstruct), which defines a map projection and related display properties for map axes. See the defaultm function reference page for information about mstructs.

The following subsections describe what the geographic data structures contain and how to create and display them. Version 1 of the toolbox used a different kind of geographic data structure (called a display structure), which had a more rigid definition.

Note Version 1 display structures are being phased out of the toolbox and are currently generated only by a few functions. You can use the updategeostruct to convert a display structure containing vector features to a geographic data structure. You must do this in order to export data imported from VMAP level 0 or DCW data sets with shapewrite, for example. For more information, see the displaym reference page, which describes the format and contents of display structures.

Representing Geographic Features

Certain Mapping Toolbox functions read, create, or manipulate vector geodata organized within a geographic data structure. This is a MATLAB structure array that has one element per geographic feature. Each feature is represented by coordinates and any kind and number of attributes. When it holds geographic coordinates (latitude and longitude), it is called a geostruct; when it contains plane (projected x and y) map coordinates, it is called a mapstruct.

Geostructs and mapstructs most frequently originate when vector geodata is imported from a shapefile. They hold only vector features and cannot be used to hold raster data, whether images, regular data grids, or geolocated data grids. The shaperead function returns either a geostruct or a mapstruct that encapsulates some or all of the data stored in a shapefile data set (which stores attributes and coordinates in separate files). The gshhs function also returns a geostruct.

By default, shaperead returns a mapstruct containing X and Y fields. This is appropriate if the data set coordinates are already projected (are in map space). Otherwise, when the coordinate data is unprojected (are in geographic space ), use the parameter-value pair 'UseGeoCoords',true to make shaperead return a geostruct having Lon and Lat fields instead. If you do not know whether coordinates in the shapefile are projected or unprojected, ask your data provider.

To determine what kinds of data a shapefile data set contains, query the data set using the shapeinfo function. shapeinfo returns a 1-by-1 structure containing high-level descriptions; it cannot be used as a geostruct or mapstruct.

Geographic data structures conform to a few simple rules, making them easy to construct programmatically:

Each array element represents an entity of the same geometric class (Point, MultiPoint, Line, or Polygon).
Geostructs: For geographic coordinates, the fields Geometry, Lat, and Lon must be present.
Mapstructs: For plane coordinates, the fields Geometry, Y, and X must be present.
Except for Point geometry, a field called BoundingBox that contains [minX minY; maxX maxY] is provided.
If present, additional fields (for attributes) must contain either strings or scalar numbers.

Coordinate data is stored in fields called X or Lon and Y or Lat, depending on what type of coordinates were read in or placed there. The mandatory coordinate field names tell map display and data export functions that accept geographic data structures whether a map projection has already been applied to coordinates, which can prevent them from generating erroneous data and graphics. Be aware that a mapstruct does not itself identify the map projection or projection parameters to which the X and Y coordinates are referenced. This is also generally true for geodata stored in shapefiles; see the shapeinfo reference page for more information.

Any number of attribute fields can be included in a geostruct or mapstruct, and can contain either a double or a string data type. This type must be consistent for a given attribute across all elements of the structure array. The fields in a geostruct representing imported data can vary according to the type of geometry and the names of attributes that have been read. If a shapefile attribute name cannot be directly used as a field name, shaperead assigns the field an appropriately mangled name, usually by substituting underscores for spaces. shaperead can filter out unwanted attributes; see Selecting Data to Read with the shaperead Function for information on how to set up such filters.

Here is an example of an unfiltered mapstruct returned by shaperead:

S = shaperead('concord_roads.shp')

S = 
609x1 struct array with fields:
    Geometry
    BoundingBox
    X
    Y
    STREETNAME
    RT_NUMBER
    CLASS
    ADMIN_TYPE
    LENGTH

This indicates that the shapefile contains 609 features. In addition to the Geometry, BoundingBox, and coordinate fields (X and Y), there are five attribute fields: STREETNAME, RT_NUMBER, CLASS, ADMIN_TYPE, and LENGTH. Looking at the 10th element,

S(10)

ans = 
       Geometry: 'Line'
    BoundingBox: [2x2 double]
              X: [1x9 double]
              Y: [1x9 double]
     STREETNAME: 'WRIGHT FARM'
      RT_NUMBER: ''
          CLASS: 5
     ADMIN_TYPE: 0
         LENGTH: 79.0347

you can see that this mapstruct contains line features, and that the tenth line has nine vertices. The first two attributes are string-valued (the second happens to be empty), and the final three attributes are numeric. Across the elements of S, X and Y can have various lengths, but STREETNAME and RT_NUMBER must always contain strings, and CLASS, ADMIN_TYPE, and LENGTH must always contain scalar doubles.

Note In mapstructs and geostructs having Line or Polygon geometries, individual features can have multiple parts—separated, disconnected line segments and polygon rings. The parts can include inner rings that run counterclockwise and outline voids. Each disconnected part is separated from the next by a NaN within the X and Y (or Lat and Lon) vectors. You can use the isShapeMultipart function to determine if a feature has NaN-separated parts.

Each NaN-separated multipart line, polygon, or multipoint entity constitutes a single feature and thus has one string or scalar double value per attribute field. It is not possible to assign distinct attributes to the different parts of such a feature; any string or numeric attribute imported with (or subsequently added to) the geostruct or mapstruct applies to all the feature's parts.

How to Construct Geostructs and Mapstructs

Functions such as shaperead or gshhs return geostructs when importing vector geodata. However, you might want to create geostructs or mapstructs yourself in some circumstances. For example, you might import vector geodata that is not stored in a shapefile (for example, from a MAT-file, from an Microsoft^® Excel^® spreadsheet, or by reading in a delimited text file). You also might compute vector geodata and attributes by calling various MATLAB or Mapping Toolbox functions. In both cases, the coordinates and other data are typically vectors or matrices in the workspace. Packaging variables into a geostruct or mapstruct can make mapping and exporting them easier, because geographic data structures provide several advantages over coordinate arrays:

All associated geodata variables are packaged in one container, a structure array.
The structure is self-documenting through its field names.
You can vary map symbology for points, lines, and polygons according to their attribute values by constructing a symbolspec for displaying the geostruct or mapstruct.
A one-to-one correspondence exists between structure elements and geographic features, which extends to the children of hggroups constructed by mapshow and geoshow.

Achieving these benefits is not difficult. Use the following example as a guide to packaging vector geodata you import or create into geographic data structures.

Example — Making Point and Line Geostructs. The following example first creates a point geostruct containing three cities on different continents and plots it with geoshow. Then it creates a line geostruct containing data for great circle navigational tracks connecting these cities. Finally, it plots these lines using a symbolspec.

Begin with a small set of point data, approximate latitudes and longitudes for three cities on three continents:

latparis =  48.87084; lonparis =   2.41306;   % Paris coords
latsant  = -33.36907; lonsant  = -70.82851;   % Santiago
latnyc   =  40.69746; lonnyc   = -73.93008;   % New York City

Build a point geostruct; it needs to have the following required fields:

Geometry (in this case 'Point')
Lat (for points, this is a scalar double)
Lon (for points, this is a scalar double)

% The first field by convention is Geometry (dimensionality).
% As Geometry is the same for all elements, assign it with deal:
[Cities(1:3).Geometry] = deal('Point');

% Add the latitudes and longitudes to the geostruct:
Cities(1).Lat = latparis; Cities(1).Lon = lonparis;
Cities(2).Lat = latsant;  Cities(2).Lon = lonsant;
Cities(3).Lat = latnyc;   Cities(3).Lon = lonnyc;

% Add city names as City fields. You can name optional fields 
% anything you like other than Geometry, Lat, Lon, X, or Y.
Cities(1).Name = 'Paris';
Cities(2).Name = 'Santiago';
Cities(3).Name = 'New York';
% Inspect your completed geostruct and its first member
Cities    

Cities = 
1x3 struct array with fields:
    Geometry
    Lat
    Lon
    Name

Cities(1)

ans = 
    Geometry: 'Point'
         Lat: 48.8708
         Lon: 2.4131
        Name: 'Paris'

Display the geostruct on a Mercator projection of the Earth's land masses stored in the landareas.shp demo shapefile , setting map limits to exclude polar regions:

axesm('mercator','grid','on','MapLatLimit',[-75 75]); tightmap; 
% Map the geostruct with the continent outlines
geoshow('landareas.shp')

% Map the City locations with filled circular markers
geoshow(Cities,'Marker','o',...
    'MarkerFaceColor','c','MarkerEdgeColor','k');

% Display the city names using data in the geostruct field Name.
% Note that you must treat the Name field as a cell array.
textm([Cities(:).Lat],[Cities(:).Lon],...
    {Cities(:).Name},'FontWeight','bold');

Next, build a Line geostruct to package great circle navigational tracks between the three cities:

% Call the new geostruct Tracks and give it a line geometry:
[Tracks(1:3).Geometry] = deal('Line');

% Create a text field identifying kind of track each entry is.
% Here they all will be great circles, identified as 'gc'
% (string signifying great circle arc to certain functions)
trackType = 'gc';
[Tracks.Type] = deal(trackType);

% Give each track an identifying name
Tracks(1).Name = 'Paris-Santiago';
[Tracks(1).Lat Tracks(1).Lon] = ...
        track2(trackType,latparis,lonparis,latsant,lonsant);

Tracks(2).Name = 'Santiago-New York';
[Tracks(2).Lat Tracks(2).Lon] = ...
        track2(trackType,latsant,lonsant,latnyc,lonnyc);

Tracks(3).Name = 'New York-Paris';
[Tracks(3).Lat Tracks(3).Lon] = ...
        track2(trackType,latnyc,lonnyc,latparis,lonparis);

Compute lengths of the great circle tracks:

% The distance function computes distance and azimuth between
% given points, in degrees. Store both in the geostruct.
for j = 1:numel(Tracks)
    [dist az] = ...
        distance(trackType,Tracks(j).Lat(1),...
                           Tracks(j).Lon(1),...
                           Tracks(j).Lat(end),...
                           Tracks(j).Lon(end));
    [Tracks(j).Length] = dist;
    [Tracks(j).Azimuth] = az;
end
% Inspect the first member of the completed geostruct
Tracks(1)

ans = 
    Geometry: 'Line'
        Type: 'gc'
        Name: 'Paris-Santiago'
         Lat: [100x1 double]
         Lon: [100x1 double]
      Length: 104.8274
     Azimuth: 235.8143

Map the three tracks in the line geostruct:

% On cylindrical projections like Mercator, great circle tracks
% are curved except those that follow the Equator or a meridian.

% Graphically differentiate the tracks by creating a symbolspec;
% key line color to track length, using the 'summer' colormap.
% Symbolspecs make it easy to vary color and linetype by
% attribute values. You can also specify default symbologies.

colorRange = makesymbolspec('Line',...
            {'Length',[min([Tracks.Length]) ...
              max([Tracks.Length])],...
             'Color',winter(3)});
geoshow(Tracks,'SymbolSpec',colorRange);

You can save the geostructs you just created as shapefiles by calling shapewrite with a filename of your choice, for example:

shapewrite(Cities,'citylocs');
shapewrite(Tracks,'citytracks');

Making Polygon Geostructs. Creating a geostruct or mapstruct for polygon data is similar to building one for point or line data. However, if your polygons include multiple, NaN-separated parts, recall that they can have only one value per attribute, not one value per part. Each attribute you place in a structure element for such a polygon pertains to all its parts. This means that if you define a group of islands, for example with a single NaN-separated list for each coordinate, all attributes for that element describe the islands as a group, not particular islands. If you want to associate attributes with a particular island, you must provide a distinct structure element for that island.

Be aware that the ordering of polygon vertices matters. When you map polygon data, the direction in which polygons are traversed has significance for how they are rendered by functions such as geoshow, mapshow, and mapview. Proper directionality is particularly important if polygons contain holes. The Mapping Toolbox convention encodes the coordinates of outer rings (e.g., continent and island outlines) in clockwise order; counterclockwise ordering is used for inner rings (e.g., lakes and inland seas). Within the coordinate array, each ring is separated from the one preceding it by a NaN.

When plotted by mapshow or geoshow, clockwise rings are filled. Counterclockwise rings are unfilled; any underlying symbology shows through such holes. To ensure that outer and inner rings are correctly coded according to the above convention, you can invoke the following functions:

ispolycw — True if vertices of polygonal contour are clockwise ordered
poly2cw — Convert polygonal contour to clockwise ordering
poly2ccw — Convert polygonal contour to counterclockwise ordering
poly2fv — Convert polygonal region to face-vertex form for use with patch in order to properly render polygons containing holes

Three of these functions check or change the ordering of vertices that define a polygon, and the fourth one converts polygons with holes to a completely different representation.

For more information about working with polygon geostructs, see the Mapping Toolbox "Converting Coastline Data (GSHHS) to Shapefile Format" demo mapexgshhs.

Mapping Toolbox™ Version 1 Display Structures

Prior to Version 2, when geostructs and mapstructs were introduced, a different data structure was employed when importing geodata from certain external formats to encapsulate it for map display functions. These display structures accommodated both raster and vector map data and other kinds of objects, but lacked the generality of current geostructs and mapstructs for representing vector features and are being phased out of the toolbox. However, you can convert display structures that contain vector geodata to geostruct form using updategeostruct. For more information about Version 1 display structures and their usage, see Version 1 Display Structures in the reference page for displaym. Additional information is located in reference pages for updategeostruct, extractm, and mlayers.

Selecting Data to Read with the shaperead Function

The shaperead function provides you with a powerful method, called a selector, to select only the data fields and items you want to import from shapefiles.

A selector is a cell array with two or more elements. The first element is a handle to a predicate function (a function with a single output argument of type logical). Each remaining element is a string indicating the name of an attribute.

For a given feature, shaperead supplies the values of the attributes listed to the predicate function to help determine whether to include the feature in its output. The feature is excluded if the predicate returns false. The converse is not necessarily true: a feature for which the predicate returns true may be excluded for other reasons when the selector is used in combination with the bounding box or record number options.

The following examples are arranged in order of increasing sophistication. Although they use MATLAB function handles, anonymous functions, and nested functions, you need not be familiar with these features in order to master the use of selectors for shaperead.

Example 1: Predicate Function in Separate File

Define the predicate function in a separate file. (Prior to Release 14, this was the only option available.) Create a file named roadfilter.m, with the following contents:

 function result = roadfilter(roadclass,roadlength)
 mininumClass = 4;
 minimumLength = 200;
 result = (roadclass  >= mininumClass) && ...
          (roadlength >= minimumLength);
 end

You can then call shaperead like this:

roadselector = {@roadfilter, 'CLASS', 'LENGTH'}

roadselector = 
    @roadfilter    'CLASS'    'LENGTH'

s = shaperead('concord_roads', 'Selector', roadselector)

s = 
115x1 struct array with fields:
    Geometry
    BoundingBox
    X
    Y
    STREETNAME
    RT_NUMBER
    CLASS
    ADMIN_TYPE
    LENGTH

or, in a slightly more compact fashion, like this:

s = shaperead('concord_roads',...
              'Selector', {@roadfilter, 'CLASS', 'LENGTH'})

s = 
115x1 struct array with fields:
    Geometry
    BoundingBox
    X
    Y
    STREETNAME
    RT_NUMBER
    CLASS
    ADMIN_TYPE
    LENGTH

Prior to Version 7 of the Mapping Toolbox software, putting the selector in a file or subfunction of its own was the only way to work with a selector.

Note that if the call to shaperead took place within a function, then roadfilter could be defined in a subfunction thereof rather than in an M-file of its own.

Example 2: Predicate as Function Handle

As a simple variation on the previous example, you could assign a function handle, roadfilterfcn, and use it in the selector:

roadfilterfcn = @roadfilter
s = shaperead('concord_roads',...
              'Selector', {roadfilterfcn, 'CLASS', 'LENGTH'})
roadfilterfcn = 
@roadfilter
s = 
115x1 struct array with fields:
    Geometry
    BoundingBox
    X
    Y
    STREETNAME
    RT_NUMBER
    CLASS
    ADMIN_TYPE
    LENGTH

Example 3: Predicate as Anonymous Function

Having to define predicate functions in M-files of their own, or even as subfunctions, may sometimes be awkward. Anonymous functions allow the predicate function to be defined right where it is needed. For example:

roadfilterfcn = ...
    @(roadclass, roadlength) (roadclass >= 4) && ...
    (roadlength >= 200)

roadfilterfcn = 
    @(roadclass, roadlength) (roadclass >= 4) ...
               && (roadlength >= 200)

s = shaperead('concord_roads','Selector', ...
              {roadfilterfcn, 'CLASS', 'LENGTH'})

s = 
115x1 struct array with fields:
    Geometry
    BoundingBox
    X
    Y
    STREETNAME
    RT_NUMBER
    CLASS
    ADMIN_TYPE
    LENGTH

Example 4: Predicate (Anonymous Function) Defined Within Cell Array

There is actually no need to introduce a function handle variable when defining the predicate as an anonymous function. Instead, you can place the whole expression within the selector cell array itself, resulting in somewhat more compact code. This pattern is used in many examples throughout Mapping Toolbox documentation and M-file help.

s = shaperead('concord_roads', 'Selector', ...
    {@(roadclass, roadlength)...
    (roadclass >= 4) && (roadlength >= 200),...
    'CLASS', 'LENGTH'})

s = 
115x1 struct array with fields:
    Geometry
    BoundingBox
    X
    Y
    STREETNAME
    RT_NUMBER
    CLASS
    ADMIN_TYPE
    LENGTH

Example 5: Parameterizing the Selector; Predicate as Nested Function

In the previous patterns, the predicate involves two hard-coded parameters (called minimumClass and minimumLength in roadfilter.m), as well as the roadclass and roadlength input variables. If you use any of these patterns in a program, you need to decide on minimum cut-off values for roadclass and roadlength at the time you write the program. But suppose that you wanted to wait and decide on parameters like minimumClass and minimumLength at run time?

Fortunately, nested functions provide the additional power that you need to do this; they allow you utilize workspace variables in as parameters, rather than requiring that the parameters be hard-coded as constants within the predicate function. In the following example, the workspace variables minimumClass and minimumLength could have been assigned through a variety of computations whose results were unknown until run-time, yet their values can be made available within the predicate as long as it is defined as a nested function. In this example the nested function is wrapped in an M-file called constructroadselector.m, which returns a complete selector: a handle to the predicate (named nestedroadfilter) and the two attibute names:

 function roadselector = ...
     constructroadselector(minimumClass, minimumLength)
 roadselector = {@nestedroadfilter, 'CLASS', 'LENGTH'};
     function result = nestedroadfilter(roadclass, roadlength)
         result = (roadclass  >= minimumClass) && ...
                  (roadlength >= minimumLength);
     end
 end

The following four lines show how to use constructroadselector:

minimumClass = 4;     % Could be run-time dependent
minimumLength = 200;  % Could be run-time dependent

roadselector = constructroadselector(...
    minimumClass, minimumLength);

s = shaperead('concord_roads', 'Selector', roadselector)

s = 
115x1 struct array with fields:
    Geometry
    BoundingBox
    X
    Y
    STREETNAME
    RT_NUMBER
    CLASS
    ADMIN_TYPE
    LENGTH

Types of Map Data Handled by the Toolbox

Understanding Raster Geodata