Accelerating the pace of engineering and science

# Documentation

## Understanding Vector Geodata

### Points, Lines, and Polygons

In the context of geodata, vector data means points, lines, and polygons that represent geographic objects. Vector geospatial data is used to represent point features, such as cities and landmarks; linear features, such as rivers and highways; and areal features, such as bodies of water and voting districts.

#### Displaying a Point

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, a single point, the location of the Eiffel Tower can be mapped as a vector.

1. Display a map of France.

h = worldmap('France');
geoshow (landareas, 'FaceColor', [1 1 .5]);
2. Save the location of the Eiffel Tower in vector form.

TowerLat = 48.85;
TowerLon = 2.28;

3. Place a red dot on the map to indicate the tower and label it.

geoshow(TowerLat, TowerLon, 'Marker','.','MarkerEdgeColor','red')
textm(TowerLat,TowerLon + 0.5, 'Eiffel Tower');

#### Displaying a Line

This simple example demonstrates how to use the function linem to display vector data for three short lines branching from one common endpoint.

axesm sinusoid; framem;
linem([15; 0; -45; -25],[-100; 0; 70; 110],'r-')
linem([15; -30; -60; -65],[-100; -20; 100; 150],'b-')
linem([15; 20; 40; 20],[-100; -20; 40; 50], 'g-')

Rivers are examples of lines. Enter the following code to add rivers to a map of France:

h = worldmap('France');
geoshow (landareas, 'FaceColor', [1 1 .5]);
geoshow(rivers, 'Color', 'blue')

Find out more about the rivers structure array:

rivers
rivers =

128x1 struct array with fields:
Geometry
BoundingBox
Lon
Lat
Name

rivers contains 128 world rivers. Type the following at the command line to view the structure for the eighth river:

rivers(8)
ans =

Geometry: 'Line'
BoundingBox: [2x2 double]
Lon: [129.6929 128.9659 128.7473 NaN]
Lat: [63.3965 63.4980 63.5220 NaN]
Name: 'Lena'

The rivers are stored as shapes of type 'Line'. Data for the eighth river, Lena, is stored in Lat and Lon vectors. Note that each vector ends with a NaN.

#### Displaying a Polygon

Many common map objects, such as state boundaries, islands, and continents, are polygons. Some polygon objects in the real world 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. A NaN separator serves as a "pen-up" command, a command to stop drawing one polygon and start drawing another. The example below demonstrates how NaN separators divide the data for simple polygons.

1. Copy and paste the following vector variables at the command line:

x = [40  55  33  10  0  5  10  40  NaN  10  25  30  25  10...
10  NaN  90  80  65  80  90  NaN];
y = [50  20  0  0  15  25  55  50  NaN  20  10  10  20  30...
20  NaN  10  0  20  25  10  NaN];

These vector variables appear as rows in the following table, along with a row listing the indices.

Notice that the NaNs appear in the same locations in both x and y vectors. Columns 9, 16, and 22 of the table have NaNs. These mark the division between separate polygons. Also, notice that the x and y data for vertices 1 and 8 are the same. This is the point where segments join together to close the polygon.

2. Use the mapshow function to display the polygon.

mapshow(x,y,'DisplayType','polygon')

3. In this example, the vector variables contain data that displays as a compact, multipart polygon with a hole. Compare the data from the table to the illustration. Note that the vertices in the illustration have been labelled to correspond to the indices. (Your output will not contain these labels.)

Individual contours in x and y are assumed to be external contours if their vertices are arranged in clockwise order; otherwise they are assumed to be internal contours. You can see that the "hole" has vertices that appear in counter-clockwise order.

Now consider an example of polygons on a map of the United States. Enter the following code to display a map of the U.S. (excluding Alaska and Hawaii):

figure; ax = usamap('conus');
set(ax, 'Visible', 'off')
names = {states.Name};
indexConus = 1:numel(states);
stateColor = [0.5 1 0.5];
geoshow(ax, states(indexConus),  'FaceColor', stateColor)
setm(ax, 'Frame', 'off', 'Grid', 'off',...
'ParallelLabel', 'off', 'MeridianLabel', 'off')

Examine the structure for one of the states:

states(4)
ans =

Geometry: 'Polygon'
BoundingBox: [2x2 double]
Lon: [1x183 double]
Lat: [1x183 double]
Name: 'Arkansas'
LabelLat: 34.8350
LabelLon: -91.8861

You can see that Arkansas is of shape type 'Polygon'. View the last entry in the Lat vector:

states(4).Lat(end)
ans =

NaN

The NaN serves as a separator between polygons.

Compare the first and next-to-last entries:

states(4).Lat(1)
states(4).Lat(182)
ans =

33.0200

ans =

33.0200

The first and next-to-last entries are the same to close the polygon.

### 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 of vector map features organized in geographic data structures.

A geographic data structure is a MATLAB structure array that has one element per geographic feature. Each feature is represented by coordinates and attributes. A geographic data structure that holds geographic coordinates (latitude and longitude) is called a geostruct, and one that holds map coordinates (projected x and y) is called a mapstruct. Geographic data structures hold only vector features and cannot be used to hold raster data (regular or geolocated data grids or images).

#### Shapefiles

Geographic data structures most frequently originate when vector geodata is imported from a shapefile. The Environmental Systems Research Institute designed the shapefile format for vector geodata. Shapefiles encode coordinates for points, multipoints, lines, or polygons, along with non-geometrical attributes.

A shapefile stores attributes and coordinates in separate files; it consists of a main file, an index file, and an xBASE file. All three files have the same base name and are distinguished by the extensions .shp, .shx, and .dbf, respectively. (For example, given the base name 'concord_roads' the shapefile filenames would be 'concord_roads.shp', 'concord_roads.shx', and 'concord_roads.dbf').

#### The Contents of Geographic Data Structures

The shaperead function reads vector features and attributes from a shapefile and returns a geographic data structure array. The shaperead function determines the names of the attribute fields at run-time from the shapefile xBASE table or from optional, user-specified parameters. If a shapefile attribute name cannot be directly used as a field name, shaperead assigns the field an appropriately modified name, usually by substituting underscores for spaces.

Fields in a Geographic Data Structure

Geometry

String

One of the following shape types: 'Point', 'MultiPoint', 'Line', or 'Polygon'.

For a 'PolyLine', the value of the Geometry field is simply 'Line'.

BoundingBox

2-by-2 numerical array

Specifies the minimum and maximum feature coordinate values in each dimension in the following form:

$\left[\begin{array}{cc}\mathrm{min}\left(X\right)& \mathrm{min}\left(Y\right)\\ \mathrm{max}\left(X\right)& \mathrm{max}\left(Y\right)\end{array}\right]$

Omitted for shape type 'Point'.

X, Y, Lon, or Lat

1-by-N array of class double

Coordinate vector.

Attr

String or scalar number

Attribute name, type, and value.

Optional. There are usually multiple attributes.

The shaperead function does not support any 3-D or "measured" shape types: 'PointZ', 'PointM', 'MultipointZ', 'MultipointM', 'PolyLineZ', 'PolyLineM', 'PolygonZ', 'PolylineM', or 'Multipatch'. Also, although 'Null Shape' features can be present in a 'Point', 'Multipoint', 'PolyLine', or 'Polygon' shapefile, they are ignored.

PolyLine and Polygon Shapes.  In geographic data structures with Line or Polygon geometries, individual features can have multiple parts—disconnected line segments and polygon rings. The parts can include counterclockwise inner rings that outline "holes." For an illustration of this, see Displaying a Polygon. 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 multipoint or NaN-separated multipart line or polygon 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 in combination.

Mapstructs and Geostructs.  By default, shaperead returns a mapstruct containing X and Y fields. This is appropriate if the data set coordinates are already projected (in a map coordinate system). Otherwise, if the data set coordinates are unprojected (in a geographic coordinate system), use the parameter-value pair 'UseGeoCoords',true to make shaperead return a geostruct having Lon and Lat fields.

Coordinate Types. If you do not know whether a shapefile uses geographic coordinates or map coordinates, here are some things you can try:

• Use shapeinfo to obtain the BoundingBox. By looking at the ranges of coordinates, you may be able to tell what kind of coordinates you have.

• Examine the optional .prj file, if one has been provided. The .prj file is written in well-known text, a text mark-up language. If your .prj file contains the term PROJCS, you have map coordinates. If your .prj file contains the term GEOGCS, but not the term PROJCS, you have geographic coordinates.

The geoshow function displays geographic features stored in geostructs, and the mapshow function displays geographic features stored in mapstructs. If you try to display a mapstruct with geoshow, the function issues a warning and calls mapshow. If you try to display a geostruct with mapshow, the function projects the coordinates with a Plate Carree projection and issues a warning.

#### Examining a Geographic Data Structure

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

The output appears as follows:

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

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.

Look at the 10th element:

S(10)

The output appears as follows:

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

This mapstruct contains 'Line' features. The tenth line has nine vertices. The first two attributes are string-valued. The second happens to be empty. 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.

In this example, shaperead returns an unfiltered mapstruct. If you want to filter out some attributes, see Selecting Data to Read with the shaperead Function for more information.

#### How to Construct Geographic Data Structures

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.

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.

1. 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
2. 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'
3. Display the geostruct on a Mercator projection of the Earth's land masses stored in the landareas.shp 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');

4. 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);
5. 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
6. 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" example mapexgshhsmapexgshhs.

#### 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.

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

1. 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:

mininumClass = 4;
minimumLength = 200;
result = (roadclass  >= mininumClass) && ...
end
2. You can then call shaperead like this:

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

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

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

Prior to Version 7 of the Mapping Toolbox software, putting the selector in a file or local function 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 local function thereof rather than in a 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:

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

#### Example 3: Predicate as Anonymous Function

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

s =
115x1 struct array with fields:
Geometry
BoundingBox
X
Y
STREETNAME
RT_NUMBER
CLASS
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 the Mapping Toolbox documentation and function help.

'CLASS', 'LENGTH'})

s =
115x1 struct array with fields:
Geometry
BoundingBox
X
Y
STREETNAME
RT_NUMBER
CLASS
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 a file called constructroadselector.m, which returns a complete selector: a handle to the predicate (named nestedroadfilter) and the two attibute names:

result = (roadclass  >= minimumClass) && ...
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

minimumClass, minimumLength);