Mathworks MAPPING TOOLBOX 3 user guide

Mapping Toolbox™ 3
User’s Guide
How to Contact The MathWorks
www.mathworks. comp.soft-sys.matlab Newsgroup www.mathworks.com/contact_TS.html T echnical Support
bugs@mathwo doc@mathworks.com Documentation error reports service@mathworks.com Order status, license renewals, passcodes
info@mathwo
com
rks.com
rks.com
Web
Bug reports
Sales, prici
ng, and general information
508-647-7000 (Phone)
508-647-7001 (Fax)
The MathWorks, Inc. 3 Apple Hill Drive Natick, MA 01760-2098
For contact information about worldwide offices, see the MathWorks Web site.
Mapping Toolbox™ User’s Guide
© COPYRIGHT 1997–20 10 by The MathWorks, Inc.
The software described in this document is furnished under a license agreement. The software may be used or copied only under the terms of the license agreement. No part of this manual may be photocopied or reproduced in any form without prior written consent from The MathW orks, Inc.
FEDERAL ACQUISITION: This provision applies to all acquisitions of the Program and Documentation by, for, or through the federal government of the United States. By accepting delivery of the Program or Documentation, the government hereby agrees that this software or documentation qualifies as commercial computer software or commercial computer software documentation as such terms are used or defined in FAR 12.212, DFARS Part 227.72, and DFARS 252.227-7014. Accordingly, the terms and conditions of this Agreement and only those rights specified in this Agreement, shall pertain to and govern theuse,modification,reproduction,release,performance,display,anddisclosureoftheProgramand Documentation by the federal government (or other entity acquiring for or through the federal government) and shall supersede any conflicting contractual terms or conditions. If this License fails to meet the government’s needs or is inconsistent in any respect with federal procurement law, the government agrees to return the Program and Docu mentation, unused, to The MathWorks, Inc.
Trademarks
MATLAB and Simulink are registered trademarks of The MathWorks, Inc. See
www.mathworks.com/trademarks for a list of additional trademarks. Other product or brand
names may be trademarks or registered trademarks of their respective holders.
Patents
The MathWorks products are protected by one or more U.S. patents. Please see
www.mathworks.com/patents for more information.
Revision History
May 1997 First printing New for Version 1.0 October 1998 Second printing Version 1.1 November 2000 Third printing Version 1.2 (Release 12) July 2002 Online only Revised for Version 1.3 (Release 13) September 2003 Online only Revised for Version 1.3.1 (Release 13SP1) January 2004 Online only Revised for Version 2.0 (Release 13SP1+) April 2004 Online only Revised for Version 2.0.1 (Release 13SP1+) June 2004 Fourth printing Revised for Version 2.0.2 (Release 14) October 2004 Online only R evised for Version 2.0.3 (Release 14SP1) March 2005 Fifth printing Revised for Version 2.1 (Release 14SP2) August 2005 Sixth printing Minor revision for Version 2.1 September 2005 Online only R evised for Version 2.2 (Release 14SP3) March 2006 Online only Revised for Version 2.3 (Release 2006a) September 2006 Seventh printing Revised for Version 2.4 (Release 2006b) March 2007 Online only Revised for Version 2.5 (Release 2007a) September 2007 Eighth printing Revised for Version 2.6 (Release 2007b) March 2008 Online only Revised for Version 2.7 (Release 2008a) October 2008 Online only R evised for Version 2.7.1 (Release 2008b) March 2009 Online only R evised for Version 2.7.2 (Release 2009a) September 2009 Online only Revised for Version 3.0 (Release 2009b) March 2010 Online only Revised for Version 3.1 (Release 2010a)
Getting Started
1
Product Overview ................................. 1-2
Contents
Dedication and Acknowledgment
Your First Maps
See the World Tour Boston with the Map Viewer
Getting Mo re Help
Ways to Get Mapping Toolbox Help Consulting Release Notes
Mapping Toolbox Demos and Data
Available Dem os Locating Geospatial Data
................................... 1-4
.................................... 1-4
................................. 1-26
........................... 1-26
.................................. 1-27
........................... 1-28
................... 1-3
.................... 1-9
................... 1-26
.................. 1-27
Understanding Map Data
2
Maps and Map Data ................................ 2-2
What Is a Map? What Is Geospatial Data?
................................... 2-2
........................... 2-2
Types of Map Data Handled by the Toolbox
Vector Ge odata Raster Geodata Combining Vector and Raster Geodata
Understanding Vector Geodata
................................... 2-4
................................... 2-7
................ 2-10
..................... 2-13
.......... 2-4
v
Points,Lines,andPolygons ......................... 2-13
Segments Versus Polygons Mapping Toolbox Geographic Data Structures Selecting Data to Read with the shaperead Function
.......................... 2-19
.......... 2-21
..... 2-32
Understanding R aster Geodata
Georeferencing R aste r Data Regular Data Grids Geolocated Data Grids
Reading and Writing Geospatial Data
Functions that Read and Write Geospatial Data Exporting Vector Geodata Functions That Read and Write Files in Compressed
Formats
....................................... 2-72
................................ 2-40
............................. 2-49
..................... 2-38
......................... 2-38
............... 2-57
........ 2-57
.......................... 2-62
Understanding Geospatial Geometry
3
Understanding Spherical Coordinates ............... 3-2
Spheres, Spheroids, and Geoids Geoid and Ellipsoid The Ellipsoid Vector
Understanding Latitude and L ongitude
................................ 3-2
............................... 3-4
...................... 3-2
............. 3-11
vi Contents
Understanding Angles, Directions, and Distances
Positions, Azimuths, Headings, Distances, Length, and
Ranges WorkingwithLengthandDistanceUnits Working with Angles: Units and Representations Working with Distances on the Sphere Angles as Binary and Formatted Numbers
Understanding Map Projections
What Is a M ap Projection? Forward a nd Inverse Projection Projection Distortions
........................................ 3-14
.............. 3-15
................ 3-23
............. 3-27
.................... 3-29
.......................... 3-29
...................... 3-30
.............................. 3-30
.... 3-14
....... 3-18
Great Circles, Rhumb Lines, and Sm all Circles ....... 3-32
Great Circles Rhumb Lines Small Circles
..................................... 3-32
..................................... 3-32
..................................... 3-33
Directions an d Areas on the Sphere and Spheroid
About Azimuths Reckoning — The Forward Problem Distance,Azimuth,andBack-Azimuth (the Inverse
Problem)
Measuring Area of Spherical Quadrangles
Planetary Almanac Data
................................... 3-38
.................. 3-38
...................................... 3-41
............. 3-44
........................... 3-46
Creating and Viewing Maps
4
Introduction to Mapping Graphics .................. 4-2
Using worldmap and usamap
Continent, Country, Region, and State Maps Made Easy Using worldmap Using usamap
.................................. 4-5
.................................... 4-7
....................... 4-4
.... 3-38
.. 4-4
Axes for Drawing Maps
What Is a Map Axes? Using Accessing and Manipulating Map Axes Properties Using the Map Limit Properties Switching Between Projections Projected and Unprojected Graphic Objects
Controlling Map Frames and Grids
The Map Frame The Map Grid
Displaying Vector Data with Mapping Toolbox
Functions
axesm ...................................... 4-13
................................... 4-48
.................................... 4-55
....................................... 4-60
............................. 4-12
.............................. 4-12
..................... 4-19
...................... 4-34
............ 4-39
.................. 4-48
....... 4-14
vii
Programming and Scripting Map Construction ......... 4-60
Displaying Vector Data as Points and Lines Displaying Vector Maps as Lines or Patches
............ 4-60
............ 4-63
Displaying Data Grids
Types of Data Grids and Raster Display Functions Fitting Gridded Data to the Graticule Using Raster Data to Create 3-D Displays
Interacting with Displayed Maps
Picking Locations Interactively Defining Small Circles and Tracks Interactively Working with Objects by Name
............................. 4-70
...... 4-70
................. 4-71
............. 4-74
.................... 4-78
...................... 4-78
........ 4-80
...................... 4-83
Making Three-Dimensional Maps
5
Sources of Terrain Data ............................ 5-2
Digital Terrain Elevation Data from NGA Digital Elev ation Model Files from US GS Determining What Elevation Data Exists for a Region
Reading Elevation Data Interactively
Extracting DEM Data with demdataui
............. 5-2
.............. 5-3
............... 5-13
................ 5-13
... 5-3
viii Contents
Determining and Visualizing Visibility Across
Terrain
Computing Line of Sight with los2
Shading and Lighting Terrain Maps
Lighting a Terrain Map Constructed from a DTED File Lighting a Global Terrain M ap with lightm and
lightmui Surface Relief Shading Colored Surface Shaded Relief Relief Mapping with Light Objects
Draping Data on Elevation Maps
......................................... 5-19
.................... 5-19
................. 5-22
....................................... 5-25
............................. 5-29
....................... 5-33
................... 5-36
.................... 5-40
.. 5-22
Draping Geoid Heights over Topography .............. 5-40
Draping Data over Terrain with Different Gridding
..... 5-43
Working with the Globe Display
What Is the Globe Display? The Globe Display Compared with the Orthographic
Projection Using Opacity and Transparency in Globe Displays Over-the-Horizon 3-D Views Using Camera Positioning
Functions Displaying a Rotating Globe
...................................... 5-50
...................................... 5-55
.................... 5-49
......................... 5-49
......................... 5-57
Customizing and Printing Maps
6
Inset M aps ........................................ 6-2
Graphic Scales
North Arrows
Thematic Maps
What Is a Thematic Map? Choropleth Maps Special Thematic Mapping Functions
.................................... 6-8
..................................... 6-14
.................................... 6-17
........................... 6-17
.................................. 6-18
................. 6-23
..... 5-52
Using Colormaps and Colorbars
Colormap for Terrain Data Contour Colo rmaps Colormaps for Political Maps Labeling Colorbars Editing Colorbars
Printing Maps to Scale
................................ 6-30
................................ 6-36
................................. 6-37
.......................... 6-27
........................ 6-32
............................. 6-38
.................... 6-27
ix
Manipulating Geospatial Data
7
Manipulating Vector Geodata ....................... 7-2
Repackaging Vector Objects Matching Line Segments Geographic Interpolation of Vectors Vector Inters ections Polygon Area Overlaying Polygons with Set Logic Cutting Polygons at the Date Line Building Buffer Zones Trimming Vector Data to a Rectangular Region Trimming Vector Data to an Arbitrary Region Simplifying Vector C oordinate Data
..................................... 7-11
............................... 7-8
......................... 7-2
........................... 7-4
.................. 7-5
................... 7-12
.................... 7-17
.............................. 7-19
......... 7-22
.......... 7-25
.................. 7-25
Manipulating Raster Geodata
Vector-to-Raster Data Conversion Data Grids as Logical Variables Data Grid Values Along a Path Data Grid Gradient, Slope , and Aspect
....................... 7-32
.................... 7-32
...................... 7-40
...................... 7-42
................ 7-44
Using Map Projections and Coordinate Systems
8
What Is a Map Projection? .......................... 8-2
Quantitative Properties of Map Projections
The Three Main Families of Map Projections
Unwrapping the Sphere to a Plane Cylindrical Projections Conic Projections Azimuthal Projections
............................. 8-5
.................................. 8-7
.............................. 8-8
................... 8-5
.......... 8-3
......... 8-5
x Contents
Projection Aspect
The Orientation Vector
.................................. 8-10
............................. 8-10
Projection Parameters ............................. 8-18
Projection Characteristics Maps Can Have
............. 8-18
Visualizing and Quantifying Projection Distortions
Displays of Spatial Error in Maps Quantifying Map Distortions at Point Locations
Accessing, Computing, and Inverting Map Projection
Data
Accessing Projected Coordinate Data Projecting Coordinates Without a Map Axes Inverse Map Projection Coordinate Transformation s
Working with the UTM System
What Is the Universal Transverse Mercator System? Understanding UTM Parameters Setting UTM Parameters with a GUI Working in UTM Without a Map Axes Mapping Across UTM Zones
Summary and Guide to Projections
............................................ 8-37
............................. 8-41
.................... 8-27
........ 8-31
................. 8-37
............ 8-39
........................ 8-45
...................... 8-51
.... 8-51
.................... 8-52
................. 8-54
................ 8-59
......................... 8-60
................. 8-63
... 8-27
Creating Web Map Service Maps
9
Introduction to Web Map Service ................... 9-2
What Web Map Service Servers Provide Basic WMS Terminology
Basic Workflow for Creating WMS Maps
Workflow Summary Creating a Map of Elevation in Europe
Searching the WMS Database
Introduction to the WMS Database Finding Temperature Data
............................ 9-4
............................... 9-5
....................... 9-8
.......................... 9-9
............... 9-2
............. 9-5
................ 9-5
................... 9-8
xi
Refining Your Search .............................. 9-11
Refining by Text String Refining by Geographic Limits
............................. 9-11
....................... 9-12
Updating Your Layer
Retrieving Your Map
Ways to Retrieve Your Map Understanding Coordinate Reference Syste m Codes Retrieving Your Map with wmsread Setting Optional Parameters Adding a Legend to Your Map Retrieving Your Map with WebMapServer.getMap
Modifying Your Request
Setting the Geographic Limits and Time Setting the Geographic Limits, Image Dimension, Style,
and Image Format
Manually Editing a URL
Overlaying Multiple Layers
Creating a Composite Map of Multiple Layers from One
Server
Combining Layers from One Server with Data from Other
Sources
DrapingTopographyandOrtho-ImageryLayersovera
DigitalElevationModelLayer
........................................ 9-42
........................................ 9-45
.............................. 9-13
............................... 9-15
......................... 9-15
.................. 9-16
........................ 9-17
....................... 9-19
............................ 9-34
............... 9-34
.............................. 9-36
........................... 9-39
......................... 9-42
..................... 9-46
..... 9-16
...... 9-28
xii Contents
Animating Data Layers
Creating Movie of Daily Planet Images for One Month Creating an Animated GIF File Animating Time-Lapse Radar Observations Displaying Animation of Radar Images over Daily Planet
Backdrop
Saving Favorite Servers
Exploring Other Layers from a Server
Writing a KML File
...................................... 9-59
............................. 9-52
...................... 9-54
............ 9-56
............................ 9-62
............... 9-64
................................ 9-67
... 9-52
Searching for Layers Outside the Database .......... 9-68
10
Hosting Your Own WMS Server
Common Problems with WMS Servers
Connection Errors Wrong Scale Problems with Geographic Limits Problems w ith Server Changing LayerName Non-EPSG:4326 Coordinate Reference Systems Map Not Returned Unsupported WMS Version Other Unrecoverable Server Errors
................................. 9-70
...................................... 9-72
................................ 9-74
..................... 9-69
............... 9-70
.................... 9-72
........... 9-73
......... 9-74
......................... 9-75
................... 9-75
Mapping Applications
Geographic Statistics .............................. 10-2
Statistics for Point Locations on a Sphere Geographic Means Geographic Standard Deviation Equal-Areas in Geographic Statistics
................................. 10-2
...................... 10-4
.............. 10-2
................. 10-7
Navigation
What Is Navigation? Conventions for Navigational Functions Fixing Position Planning the Shortest Path Track Laydown – Displaying Navigational Tracks Dead Reckoning Drift Correction Time Zones
........................................ 10-11
............................... 10-11
............... 10-12
................................... 10-13
......................... 10-25
................................... 10-31
................................... 10-36
....................................... 10-38
....... 10-29
xiii
11
Map Projections Reference
Cylindrical Projections ............................ 11-2
Pseudocylindrical Projections
Conic Projections
Polyconic and Pseudoconic Projections
Azimuthal, Pseudoazimuthal, and Modified Azimuthal
Projections
UTM and UPS Systems
3-D Globe Display
.................................. 11-4
..................................... 11-4
............................. 11-5
.................................. 11-5
...................... 11-2
............. 11-4
xiv Contents
12
A
B
Map Projections — Alphabetical List
Glossary
Bibliography
Examples
Your First Maps ................................... B-2
Understanding Vector Geodata
Understanding R aster Geodata
Combining Vector and Raster Geodata
Geolocated Data Grids
Exporting Vector Geodata
Understanding G eosp atial Geometry
Creating an d Viewing Maps
Making T hree-Dimensional Maps
............................. B-3
..................... B-2
..................... B-2
.......................... B-3
........................ B-3
................... B-4
.............. B-2
................ B-3
xv
Customizing and Printing Maps ..................... B-4
Using Colormaps and Colorbars
Vector Data Manipulation
Raster Data Manipulation
Projections and Transformations
Web M ap Service Maps
Navigation
........................................ B-7
.......................... B-5
.......................... B-5
............................. B-6
.................... B-5
................... B-6
Index
xvi Contents

Getting Started

“Product Overview” on page 1-2
“Dedication and Acknowledgment” on page 1 -3
“Your First Maps” on page 1-4
“Getting More Help” on page 1-26
“Mapping Toolbox Demos and Data” on page 1-27
1
1 Getting Started

Product Overview

Mapping To olbox™ provides tools and utilities for analyzing geog raphic data and creating map displays. You can import ve ctor and raster data from shapefile, GeoTIFF, SDTS DEM, and other file formats, as well as Web-based data from Web Map Service (WMS) servers. The toolbox lets you customize the imported data by subsetting, trimming, intersecting, adjusting spatial resolution, and applying other methods. Geographic data can be combined with base map layers from multiple sources in a single map display. With function-level access to all key features, y ou can automate frequent tasks in your geospatial workflow .
Briefly summarized, the toolbox provides functionality in the following areas:
Vector and raster data im port and export from standard formats and
Data retrieval from Web Map Service (W MS) servers for customized
specific data products
geographic datasets and related metadata
1-2
Digital terrain and elevation model analysis functions, including profile, gradient, line-of-sight, and viewshed calculations
Geometric geodesy, including distance and area calculations, 3D coordinate transformations, and more than 65 map projections
Utilities for converting units, adjusting spatial resolution, wrapping longitudes, and managing spatially referenced images and raster data
2D and 3D map display, customization, and interaction
This chapter provides step-by-step examples of basic Mapping Toolbox capabilities and guides you toward demos and documentation that can help answer your questions. For a complete classified list of Mapping Toolbox functions and features, see “Function Reference”.

Dedication and Acknowledgment

In memory of John P. Snyder (1926–97), whose meticulous studies and systematic des criptions of map projections inspired and enabled the creation of Mapping Too lbox software.
This software was originally developed and maintained through Version 1.3 by Systems Planning and Analysis, Inc. (SPA), of Alexandria, Virginia.
Except where noted, the information contained in demo and sample data files (found in digitaldatasets. ThesedatafilesareprovidedasaconveniencetoMapping Toolbox users. The MathWorks™ makes no claims that any of this data is free of defects or errors, or that the representations of geographic features or names are up to date or authoritative.
toolbox/map/mapdemos) is derived from publ icly available
Dedication and Acknowledgment
1-3
1 Getting Started

You r Fir st Map s

In this section...
“See the World” on page 1-4 “Tour Boston with the Map Viewer” on page 1-9
This section helps you exercise high-level functions and GUIs to map and display geodata. It introduces describes how to use the Map Viewer (
worldmap and other basic functions, and then
mapview).

See the World

Spatial data refers to data describing location,shape,andspatial relationships. Geospatial data is spatial data that is in some way georeferenced, or tied to specific locations on, under, or above the surface of a planet.
Geospatial data can be voluminous, complex, and difficult to process. Mapping Toolbox functions handle many of the details of loading and displaying data, and built-in data structures facilitate data storage. Nevertheless, the more you understand about your data and the capabilities of the toolbox, the more interesting applications you can pursue, and the more useful your results will be to you and others.
1-4
Followthisexampletocreateyourfirstworldmap.
1 In the MATLAB
worldmap world
This cre ates an empty map axes, ready to hold the data of your choice.
®
Command Window, type
Your F i rst M aps
Function worldmap automatically selects a reasonablechoiceforyourmap projection and coordinate limits. In this case, it chose a Robinson projection centered on the prime meridian and the equator (0º latitude, 0º longitude).
If you type you can select a country, continent, or region. The
worldmap without an argument, a list box appears from which
worldmap function then
generates a map axes with appropriate projection and map limits.
2 Import low-resolution world coastlines—a set of discrete vertices that, when
connected in the order given, approximate the coastlines of continents, major islands, and inland seas. The vertex latitudes and longitudes are stored as MATLAB vectors in a MAT-file. First, list the variables in the file:
whos -file coast.mat
The output is shown below:
Name Size Bytes Class Attributes
lat 9865x1 78920 double long 9865x1 78920 double
3 Load and plot the coastlines on the world map:
load coast plotm(lat, long)
1-5
1 Getting Started
The plotm function is a geographic equivalent to the MATLAB plot function. It accepts coordinates in latitude and longitude, transforms them to x and y via a specified map projection, and displays them in a figure axes. In this example,
worldmap specifies the projection.
Note Certain Mapping Toolbox functions that end with m,suchasplotm and textm, are modeled after familiar MATLAB functions that handle non-geographic coordinate data.
1-6
4 Notice how the world coastlines form distinct polygons, even though only
a single vector of latitudes and a corresponding vector of longitudes are provided. The display breaks into separate parts like this because in the vectors
lat and long the vertices of various coastlines are concatenated
together but separated by isolated NaN-valued elements. In other words, “NaN separators” implicitly divideeachvectorintomultipleparts.
long include “NaN terminators” as well as separators, showing that the coast data set is organized into precisely 241 polygons.
lat and
Enter the follow ing code to break out your data into its NaN-separated parts:
[latcells, loncells] = polysplit(lat, long); numel(latcells)
latcells and loncells are cell vectors, with each cell containing the
vertices of just one polygon part. The number of parts appears:
Your F i rst M aps
ans =
241
5 Now create a new map axes for plotting data over Europe, and this time
specify a return argument:
h = worldmap('Europe');
For the map of the world, worldmap chose a pseudocylindrical R obinson projection. For Europe, it chose an Equidistant Conic projection. How can you tell which projection
When you specify a return argument for mapping functions, a handle (e.g.,
worldmap is using?
worldmap and certain other
h) to the figure’s axes is returned.
The axes object on which map data is displayed is called a map axes. In addition to the graphics properties common to any MATLAB ax es object, a map axes object contains additional properties covering map projection type, projection parameters, map limits, etc. The
getm and setm functions
and others allow you to access and modify these properties.
6 To inspect the MapProjection property for the map of Europe, type:
getm(h, 'MapProjection')
1-7
1 Getting Started
Thetypeofprojectionappears:
ans = eqdconic
If you’re not familiar with the abbreviation “eqdconic,” type:
help eqdconic
To see all the map-specific properties for this map axes, type:
getm(h)
7 Add data to the map of Europe by using the geoshow function to import and
display several shapefiles in the
geoshow('landareas.shp', 'FaceColor', [0.15 0.5 0.15]) geoshow('worldlakes.shp', 'FaceColor', 'cyan') geoshow('worldrivers.shp', 'Color', 'blue') geoshow('worldcities.shp', 'Marker', '.',...
toolbox/map/mapdemos directory:
'Color', 'red')
1-8
Note how geoshow can plot data directly from files onto a map axes without first importing it into the MATLAB workspace.
8 Finally, place a label on the map to identify the Mediterranean Sea.
labelLat = 35; labelLon = 14; textm(labelLat, labelLon, 'Mediterranean Sea')
Your F i rst M aps
To learn more about display properties for map axes and how to control them, see “Accessing and Manipulating Map Axes Properties” on page 4-14.

Tour Boston with the Map Viewer

The Map Viewer is an interactive tool for browsing map data. With it you can assemble layers of vector and raster geodata and render them in 2-D. You can import,reorder,symbolize,hide,anddelete data layers; identify coordinate locations; and list data attributes. You can display selected data attributes as datatips (signposts that identify attribute values, such as place names or route numbers). The following exercise shows how the Map Viewer works andwhatitcando.
A Map Viewer Session
1 Start a Map Viewer session by typing
1-9
1 Getting Started
mapview
at the MATLAB prompt. The Map Viewer opens with a blank canvas. (No data is present.) The viewer and its tools are shown below.
Info
Print figure
Select area
Zoom out
Insert text
Insert line
Fit to window
Zoom in
Select annotations
X and Y coordinate readouts
Pan
Insert arrow
Map scale Coordinate unit
Datatips
Prior view
drop-down
Currently active layer drop-down
Most of the tool buttons can also be activated from the Tools menu.
In the Map Viewer graphic, you can see arrows pointing to the X and Y coordinate readouts. The Map Viewer is designed primarily for working with data sets that refer to a projecte d map coordinate system (as opposed to a geographic, latitude-longitude system), so the coordinate axes are named X and Y.
2 For ease in importing Mapping Toolbox demo data, set your working
directory. Type the following in the Command Window:
1-10
cd(fullfile(matlabroot,'toolbox','map','mapdemos'))
Your F i rst M aps
3 You can also navigate to this directory with the Map Viewer Import Data
dialog if you prefer. In the Map Viewer, select the File menu and then choose Im port From File.OpentheGeoTIFFfile
boston.tif,asshown
below.
ThefileopensintheMapViewer.Theimageisavisiblered,green,andblue composite from a georeferenced IKONOS-2 panchromatic/multispectral product created by GeoEye™. Copyright © GeoEye, all rights reserved. For further information about th e image, refer to the text files
boston.txt
and boston_metadata.txt.Toopenboston.txt, type the following at the command line:
open 'boston.txt'
4 To see the map scale in the Map Viewer, set the map distance units. Use
the drop-down Map units menu at the bottom center to select
.
Feet
US Survey
1-11
1 Getting Started
5 Now set the scale to 1:25,000 by typing 1:25000 in the Scale box, which is
above the Map units drop-down. The Map Viewer now looks like this.
Map scale
1-12
The cursor in the picture has been placed so that it points at the front of the Massachusetts State House (capitol building). The map coordinates for this location are shown in the readout at the lower left as
774,114.36
feet easting (X), 2,955,685.56 feet northing (Y), in Massachusetts State Plane coordinates.
6 Next, enter the following code to import a vector data layer, the streets and
highways in the central Boston area:
boston_roads = shaperead('boston_roads.shp');
The shaperead function reads the line shapefile boston_roads.shp into the workspace as a geographic data structure. Asisfrequentlythecase
Your F i rst M aps
when overlaying geodata, the coordinate system used by boston_roads.shp (in units of mete rs) does not completely agree with the one for the satellite image,
boston.tif (inunitsoffeet). Ifyouweretoignorethis,thetwo
data se ts would be out of registration by a large distance.
7 Convert the X and Y coordinate fields of boston_roads.shp from meters to
U.S. survey feet:
surveyFeetPerMeter = unitsratio('survey feet','meter'); for k = 1:numel(boston_roads)
boston_roads(k).X = surveyFeetPerMeter * boston_roads(k).X; boston_roads(k).Y = surveyFeetPerMeter * boston_roads(k).Y;
end
The unitsratio function computes conversion factors between a variety of units of length.
8 In the Map Viewer File menu, select Im port From Workspace > Vector
Data > Geographic Data Structure.Specify
boston_roads as the data
to import from the workspace, and click OK.
1-13
1 Getting Started
1-14
You could clear the workspace now if you wanted, because all the data that the Map Viewer needs is now loaded into it.
9 AftertheMapViewerfinishesimporting the roads layer, it selects a
random color and renders all the shapes with that color as solid lines. The view looks like this.
Active layer
Your F i rst M aps
Being random, the color you see for the road layer may differ.
10 Use the Active layer drop-down menu at the bottom right to select
boston_roads. Changing the active layer has no visual effect. Doing so
allows you to query attributes of the layer you select. You can designate any layer to be the active layer; it does not need to be the topmost layer. By default, the first layer imported is active.
11 OnewaytoseetheattributesforavectorlayeristousetheInfo tool, a
button near the right end of the toolbar. Select the Info tool and click somewhere along the bridge across theCharlesRivernearthelowerleft of the map. This opens a text window displaying the attributes of the selected object.
1-15
1 Getting Started
Info
1-16
The selected road is Massachusetts Avenue (Route 2A). As the above figure shows, the
INDEX attribute added by the Map Viewer.
12 Get information about some other roads. Dismiss open Info windows by
boston_roads vectors have six attributes, including an implicit
clicking their close boxes.
13 Choose an attribute for the Datatip tool to inspect. From the Layers
menu, select boston_roads > Set Label Attribute.Fromthelistinthe Attribute Names dialog, select
CLASS and click OK to dismiss it.
Your F i rst M aps
14 Select the
dialog box dismiss th
15 Use the Datatip tool to identify the administrative class of any road
Datatip tool. The cursor assumes a crosshairs (
appearstoremindyouhowtochangeattributes. ClickOK to
ebox.
+)shape.A
displayed. When you click on a road segment, a datatip is left in that place to indicate the
CLASS attribute of the active layer, as illustrated below.
1-17
1 Getting Started
Datatip
1-18
16 You can change how the roads are rendered by identifying an attribute to
which to ke y line symbology. Color roads according to their
CLASS attribute,
whichtakesonthevalues1:6. Dothisbycreatingasymbolspec in the workspace. A symbolspec is a cell array that associates attribute names and values to graphic properties for a specified geometric class (
'MultiPoint', 'Line', 'Polygon',or'Patch'). To create a symbolspec for
line objects (in this case roads) that have a
roadcolors = makesymbolspec('Line', ... {'CLASS',1,'Color',[1 1 1]}, {'CLASS',2,'Color',[1 1 0]}, ... {'CLASS',3,'Color',[0 1 0]}, {'CLASS',4,'Color',[0 1 1]}, ... {'CLASS',5,'Color',[1 0 1]}, {'CLASS',6,'Color',[0 0 1]})
CLASS attribute, type:
'Point',
The following output appears:
roadcolors =
Your F i rst M aps
ShapeType: 'Line'
Color: {6x3 cell}
17 The Map Viewer recognizes and imports symbolspecs from the workspace.
To apply the one you just created, from the Layers menu, select boston_roads > Set Sym bol Spec. From the Laye r Symbols dialog, sele ct
the
roadcolors symbolspec you just created and click OK.AftertheMap
Viewer has read and applied the symbolspec, the map looks like this.
18 Remove the datatips before going on. To dismiss datatips, right-click one
of them and select Delete all datatips from the pop-up context menu that appears.
19 Add another layer, a set of points that identifies 13 Boston landmarks. As
you d id with the
boston_roads layer, import it from a shapefile:
1-19
1 Getting Started
boston_placenames = shaperead('boston_placenames.shp');
20 The locations for these landmarks are given in meters, so you must convert
their coordinates to units of survey feet before importing them into Map Viewer:
surveyFeetPerMeter = unitsratio('survey feet','meter'); for k = 1:numel(boston_placenames)
boston_placenames(k).X = ...
surveyFeetPerMeter * boston_placenames(k).X;
boston_placenames(k).Y = ...
surveyFeetPerMeter * boston_placenames(k).Y;
end
21 From the File menu, select Impo rt From Workspace > Vector Da ta
> Geographic Data Structure.Choose
boston_placenames as the data
to import from the workspace and click OK.
22 The boston_placenames markers are symbolized as small x markers, but
these markers do not show up over the orthophoto. To solve this problem, create a symbolspec for the markers to represent them as red filled circles. At the MATLAB command line, type:
1-20
places = makesymbolspec('Point',{'Default','Marker','o', ... 'MarkerEdgeColor','r','MarkerFaceColor','r'})
The following output appears:
places =
ShapeType: 'Point'
Marker: {'Default' '' 'o'} MarkerEdgeColor: {'Default' '' 'r'} MarkerFaceColor: {'Default' '' 'r'}
The Default keyword causes the specified symbol to be applied to all point objects in a given layer unless specifically overridden by an attribute-coded symbol in the same or a different symbolspec.
23 To activate this symbolspec, pull down the Layers menu, select
boston_placenames, slide right, and select Set Symbol Spec .Inthe
Layer Symbols dialog that appears, highlight
places and click OK.
Your F i rst M aps
The Map Viewer reads the workspace variable
places; the cross marks
turn into red circles. Note that a layer need not be active in order for you to apply a symbolspec to it.
24 Use the Active layer drop-down menu to make boston_placenames the
currently active layer, and then select the Datatip tool. Click any red circletoseethenameofthefeatureitmarks.Themaplookslikethis (depending on which datatips you show).
25 Zoom in on Beacon Hill for a closer view of the Massachusetts State House
and Boston Common. Select the Zoom in tool; move the (magnifier) cursor until the X readout is approximately
774,011 and the Y readout is roughly
1-21
1 Getting Started
2,955,615;andclickoncetoenlargetheview.Thescalechangestoabout 1:12,500 andthemapappearsasbelow.
Zoom in
1-22
26 From the Tools menu, choose Select Annotations to change from the
Datatip tool back to the original cursor. Right-click any of the data tips
and select Delete all datatips from the pop-up context menu. This clears the place names you added to the map.
27 Select an area of interest to save as a n image file. Click the Select
area tool, and then hold the mouse button down as you draw a selection
rectangle. If you do not like the selection, repeat the operation until you are satisfied. If you know what ground coordinates you want, you can use the coordinate readouts to make a precise selection. The selected area appears as a red rectangle.
Your F i rst M aps
28 In orde
direct
29 Save your selection as an image file. From the File menu, select Save As
rtobeabletosaveafileinthenextstep,changeyourworking
ory to a writable directory, such as
/work.
Raster Map > Selected Area to open an Export to File dialog.
In the Export to File dialog, navigate to a directory where you want to save the map image, and save the selected area’s image as a it
central_boston.tif. (PNG and JPG formats are also available.) A
worldfile,
Whe
An i
central_boston.tfw, is created along with the TIF.
never you save a raster map in this manner, two files are created:
mage file (
file.tif, file.png,orfile.jpg)
.tif file, calling
1-23
1 Getting Started
An accompanying worldfile that georeferences the image (file.tfw,
file.pgw,orfile.jgw)
The following steps show you how to read worldfiles and display a georeferenced image outside of
30 Read in the saved image and its colormap with the MATLAB
function
central_boston.tfw with worldfileread, and display the map with mapshow:
imread, create a referencing matrix for it by reading in
[X cmap] = imread('central_boston.tif'); R = worldfileread('central_boston.tfw'); figure mapshow(X, cmap, R);
mapview.
1-24
See the documentation for mapshow for another example of displaying a georeferenced image.
Your F i rst M aps
31 Experiment with other tools and menu items. For example, you can
annotate the map with lines, arrows, and text; fit the map to the window; draw a bounding box for any layer; and print the current view. You can also spawn a new Map Viewer using New View from the File menu. A new view can duplicate the current view, cover the active layer’s extent, cover all layer extents, or include only the selected area, if any.
When you are through with a viewing session, close the Map Viewer using the window’s close box or select Close from the File menu. For more information about the Map Viewer, see the
32 The two examples in the Getting Started chapter provide an introduction
mapview reference page.
to the Mapping Toolbox. To find out what else the toolbox can do, run the demos described in “Mapping Toolbox Demos and Data” on page 1-27. If you have a specific task in mind and want to see a similar problem solved, search the index of examples.
1-25
1 Getting Started

Getting More Help

In this section...
“Ways to Get Mapping Toolbox Help” on page 1-26 “Consulting Release Notes” on page 1-26

Ways to Get Mapping Toolbox Help

The Mapping Toolbox documentation is available in electronic form as PDF and HTML files through the want to print the chapters to browse through them. This is best done from the PDF version, available at the MathWorks Web site,
http://www.mathworks.com/access/helpdesk/help/pdf_doc/map/map_ug.pdf.
Help is available for individual commands and classes of Mapping Toolbox commands:
helpdesk command. You might
help functionname for help on a specific function, often including exam ples
1-26
doc functionname to read a function’s reference page in the Help browser,
including examples and illustrations
help map for a list of functions by category mapdemos for a list of Mapping Toolbox demos
maps to see a list of all Mapping Toolbox map projections by class, name,
and ID string
maplist to return a structure describing all Mapping Toolbox map
projections
projlist to list map projections supported by projfwd and projinv

Consulting Release Notes

To learn how one version of Mapping Toolbox software differs from the next, read the Mapping Toolbox Release No tes, which include information on enhancements, syntax and GUI changes, known software and documentation problems, and compatibility issues.

Mapping Toolbox Demos and Data

In this section...
“Available Demos” on page 1-27 “Locating Geospatial Data” on page 1-28

Available Demos

Try out some demos to see Mapping Toolbox functions in action. Two of the demos are Web demos, and you cannot view their code. To see the code for any of the other demos, open the demo and click Open <demo>.m in the Editor. This link is located on the left side of the banner at the top of the page.
Creating Map Displays
Importing Geographic Data and Creating Map Displays (Web)
Creating Maps from Geographic (Latitude, Longitude) Data — mapexgeo
Mapping Toolbox™ Demos and Data
Creating Maps from Data in a Projected x-y System — mapexmap
Creating an Interactive Map for Selecting Point Features — mapexfindcity
Using G
Conve
Expor
Plott
Un-P
eospatial Analysis Tools and Formats
rting Coastline Data (GSHHS) to Shapefile Format — mapexgshhs
ting Vector Point Data to KML — mapexkmlexport
ing a 3-D Dome as a Mesh Over a Globe — mapex3ddome
rojecting a Digital Elevation Model (DEM) — mapex unprojectdem
Working with Georeferenced Images
Creating a Half-Resolution Georeferenced Image — mapexrefmat
Georeferencing an Imag e to an Orthotile Base Layer — mapexreg
1-27
1 Getting Started
Interacting with Web Map Service (WMS) Servers
Overlaying Web Map Service (WMS) Weather and Satellite Imagery Layers
to Calculate Storm Area (Web)
Compositin g and Animating Web Map Service (WMS) Meteorological Layers — mapexwmsanimate
To see this list o f links, as well as descriptions of the sample data provi ded in Mapping Toolbox, type the following:
help mapdemos

Locating Geospatial Data

Find sample data sets in the Mapping Toolbox mapdemos directory (
toolbox\map\mapdemos). Most of the sample data sets have .txt files that
provide information or metadata about their source and content.
For information on locating digital map data you can download over the Internet, see the following documentation at the MathWorks Web site.
http://www.mathworks.com/support/tech-notes/2100/2101.html
1-28

Understanding Map Data

“Maps and Map Data” on page 2-2
“Types of Map Data Handled by the Toolbox ” on page 2-4
“Understanding Vector Geodata” on page 2-13
“Understanding Raster Geodata” on page 2-38
“Reading and Writing Geospatial Data” on page 2-57
2
2 Understanding Map Data

Maps and Map Data

In this section...
“What Is a Map?” on page 2 -2 “What Is Geospatial Data?” on page 2-2

What Is a Map?

Mapping Toolbox software manipulates electronic representations of geographic data. It lets you import, create, use, and present g eographic data in a variety of forms and to a variety of ends. In the digital network era, it is easy to think of geospatial data as maps and maps as data, but you should take care to note the differences between these concepts.
The simplest (although perhaps not the most general) definition of a map is a representation of geographic data. Most people today generally think of maps as two-dimensional; to the ancient Egyptians, however, maps first took the form of lists of place names in the order they would be encountered when following a given road. Today such a list would be considered as map data rather than as a map. When most people hear the word “map” they tend to visualize two-dimen si on al renditions such as printed road, political, and topographic maps, but even classroom globes and computer graphic flight simulation scenes are maps under this definition.
2-2
In this toolbox, map data is any variable or set of variables representing a set of geographic locations, properties of a region, or features on a planet’s surface, regardless of how large or complex the data is, or how it is formatted. Such data can be rendered as m aps in a variety of ways using the functions and user interfaces provided.

What Is Geospatial Data?

Geospatial data comes in many forms and formats, and its structure is more complicated than tabular or even nongeographic geometric data. It is, in fact, a subset of spatial data, which is simply data that indicates where things are within a given coordinate system. Mileposts on a highway, an e ngine ering drawing of an automobile part, and a rendering of a building elevation all have coordinate systems, and can be represented as spatial data when
Maps and Map Data
properly quantified (digitized). Such coordinate system s, however, are local and not explicitly tied or oriented to the Earth’s surface; thus, most digital representations of mileposts, machine parts, and buildings do not qualify as geospatial data (also called geodata).
What sets geospatial data apart from other spatial data is that it is absolutely or relatively positioned on a planet, or georeferenced. That is, it has a terrestrial coordinate system that can be shared by other geospatial data. There are many ways to define a terrestrial coordinate system and also to transform it to any number of local coordinate systems, for example, to create a map projection. However, most are based on a framework that represents a planet as a sphere or spheroid that spins on a north-south axis, and which is girded by an equator (an imaginary plane midway between the poles and perpendicular to the rotational axis).
Geodata is coded for computer storage and applications in two principal ways: vector and raster representations. It has been said that “raster is faster but vector is corrector.” There is truth to this, but the situation is more complex. The following discussions explore these two representatio ns: how they differ, what data structures support them, why you would choose one over the other, and how they can work together in the toolbox. The conclude by summarizing the functions available for importing and exporting geospatial data formats.
2-3
2 Understanding Map Data

Types of Map Data Handled by the Toolbox

In this section...
“Vector Geodata” on page 2-4 “Raster Geodata” on page 2-7 “Combining Vector and Raster Geodata” on page 2-10

Vector Geodata

Vector data (in the computer graphics sense rather than the physics sense) can represent a map. Such vectors take the form of sequences of latitude-longitude or projected coordinate pairs representing a point set, a linear map feature, or an areal map feature. For example, points delineating the boundary of the United States, the interstate highway system, the centers of major U.S. cities, or even all three sets taken together, can be used to make a map. In such representations, the geographic data is in vector format and displays of it are referred to as vector maps. Such data consists of lists of specific coordinate locations (which, if describing linear or areal f ea tures, are normally points of inflection where line direction changes), along with some indication of whether each is connected to the points adjacent to it in the list.
2-4
In the Mapping Toolbox environment, vector data consists of sequentially ordered pairs of geographic (latitude, longitude) or projected (x,y) coordinate pairs (also called tuples). Successive pairs are assumed to be connected in sequence; breaks in connectivity must b e delineated by the creation of separate vector variables or by inserting separators (usually NaNs) into the sets at each breakpoint. For vector map data, the connectivity (topological structure) of the data is often only a concern during display, but it also affects the computation of statistics such as length and area.
ALookatVectorData
1 To inspect an example of vector ma p data, enter the following commands:
load coast whos
Name Size Bytes Class Attributes
Types of Map Data Handled by the Toolbox
ans 1x45 90 char lat 9589x1 76712 double long 9589x1 76712 double
The variables lat and long are vectors in the coast MAT-file, which together form a vector map of the coastlines of the world.
2 To view a map of this data, enter these commands:
axesm mercator framem plotm(lat,long)
Inspect the first 20 coordinates of the coastline vector data:
[lat(1:20) long(1:20)]
2-5
2 Understanding Map Data
ans =
-83.83 -180
-84.33 -178
-84.5 -174
-84.67 -170
-84.92 -166
-85.42 -163
-85.42 -158
-85.58 -152
-85.33 -146
-84.83 -147
-84.5 -151
-84 -153.5
-83.5 -153
-83 -154
-82.5 -154
-82 -154
-81.5 -154.5
-81.17 -153
-81 -150
-80.92 -146.5
2-6
Does this give you any clue as to which continent’s coastline these locations represent?
3 To see the coastline these vector points represent, type this command to
display them in red:
plotm(lat(1:20), long(1:20),'r')
As you may have deduced by looking at the first column of the data, there is on ly one continent that lies below -80º latitude, Antarctica.
The above example presents the map in a Mercator projection. A map projection displays the surface of a sphere (or a spheroid) in a two-dimensional plane. As the word “plane” indicates, points on the s phere are geometrically projected to a plane surface. There are many possible ways to project a map, all of which introduce various types of distortions.
For further information on how Mapping Toolbox software manages map projections, see Chapter 8, “Using Map Projections and Coordinate Systems”.
Types of Map Data Handled by the Toolbox
For details on data structures that the toolbox uses to represent vector geodata, see “Mapping Toolbox Geographic Data Structures” on page 2-21.

Raster Geodata

You can also map data represented as a matrix (a 2-D MATLAB array) in which each row-and-column element corresponds to a rectangular patch of a specific geographic area, with implied topological connectivity to adjacent patches. This is commonly referred to as raster data. Raster is actually a hardware term meaning a systematic scan of an image that encodes it into a regular grid of pixel values arrayed in rows and columns.
When data in raster format represents the surface of a planet, it is called a data grid, and the data is stored as an array or matrix. The toolbox leverages the power of MATLAB matrix manipulation in handling this type of map data. This documentation uses the terms raster data and data grid interchangeably to talk about geodata stored in two-dimensional array form.
A raster can encode either an average value across a cell or a value sampled (posted) at the center of that cell. While geolocated data grids explicitly indicate which type of values are present ( see “Geolocated Data Grids” on page 2-49), external metadata/user knowledge is required to be able to specify whether a regular data grid encodes averages or samples of values.
Digital Elevation Data
When raster geodata consists of surface elevations, the map can also be referred to as a digital elevation model/matrix (D EM), and its display is a
topographical map. The DEM is one of the most common forms of digital terrain model (DTM), which can also be represented as contour lines,
triangulated elevation points, quadtrees, octtrees, or otherwise.
The
topo global terrain data is an example of a DEM. In this 180-by-360
matrix, each row represents one degree of latitude, and each column represents one degree of longitude. Each element of this matrix is the average elevation, in meters, for the one-degree-by-one-degree region of the Earth to which its row and column correspond.
2-7
2 Understanding Map Data
Remotely Sensed Image Data
Raster geodata also encompasses georeferenced imagery. Like data grids, images are organized into rows and columns. There are subtle distinctions, however, w hich are important in certain contexts. One distinction is that an image may contain RGB or multispectral channels in a single array, so that it has a third (color or spectral) dimension. In this case a 3-D array is used rather than a 2-D (matrix) array. Another distinction is that while data grids are stored as class double in the toolbox, images may use a range of MATLAB storage classes, with the most common being
logical. Finally, for grayscale and RGB images of class double, the values of
individual array elements are constrained to the interval
In terms of geore feren cing— conv erting between column/row subscripts and 2-D map or geographic c oord in a tes—images and data grids behave the same way (which is why both are considered to be a form of raster geodata). However, when performing operations that process the values raster elements themselves, including most display functions, it is important to be aware of whether you are working with an image or a data grid, and for images, how spectral data is encoded.
uint8, uint16, double,and
[0 1].
2-8
For further details concerning the structure of raster map data, see “Understanding Raster Geodata” on page 2-38.
A Look at Raster Data
1 Load the topo data grid.
load topo topo
2 topo contains raster elevation data. Create a referencing matrix to
georeference
topoR = makerefmat('RasterSize', size(topo), ...
'Latlim', [-90 90], 'Lonlim', [0 360]);
3 Create an equal-area map projection to view the topographic data:
axesm sinusoid
topo.
Types of Map Data Handled by the Toolbox
A figure window is created with map axes set to display a sinusoidal projection.
4 Generate a shaded relief m ap. You can do this in several ways. First use
geoshow and apply a topographic colormap using demcmap:
geoshow(topo,topoR,'DisplayType','texturemap') demcmap(topo)
The geoshow function displays geodata in geographic (unprojected) coordinates. The
geoshow output is shown below:
5 Now create a new figure using a Hammer projection (which, like the
sinusoidal, is also equal-area), and display
topo using meshlsrm,which
enables control of lighting effects:
figure; axesm hammer meshlsrm(topo,topoR)
A colored relief map of the topo data set, illuminated from the east, is rendered in the second figure window.
2-9
2 Understanding Map Data
For additio and Lightin
Note that t completel thesameno symboliz this.
Combini
Vector m togeth displa When se be refe must u same e Geosp
nal details on controlling the illumination of maps, see “Shading
g Terrain Maps” on page 5-22.
he content, symbolization, and the projection of the map are
y independent. The structure and content of the
matter how you display it, although how it is projected and
ed can affect its interpretation. The following exam ple illustrates
topo variable are
ng Vector and Raster Geodata
ap variables and data grid variables are often used or displayed
er. For example, continental coastlines in vector form might be
yed with a grid of temperature data to make the latter more useful.
veral map variables are u sed together, regardless of type, they can
rred to as a single map. To do this, of course, the different data sets
se the same coordinate system (i.e., geographic coordinates on the
llipsoid or an identical map projection). See Chapter 3, “Understanding
atial Geometry” for an introduction to these concepts.
Viewing Raster and Vector Data on the Same Map
gthe
Usin them
coast and topo data from the previous examples, you can combine
in a single map and see how well the two types of data work together:
2-10
ar the current map:
1 Cle
clma
Types of Map Data Handled by the Toolbox
2 Reload the coastline data:
load coast
3 If the topo data is n ot already in the workspace, load it as well:
load topo
4 Set up a Robinson projection:
axesm robinson
5 Plot the raster topographic data with an approp
geoshow(topo,topolegend,'DisplayType','texturemap') demcmap(topo)
6 Plot the coastline data in red on top of the terrain map:
geoshow(lat,long,'Color','r')
riate colormap:
Note that you can use geoshow to display both raster and vector data. Here is the resulting map.
For additional details on how Ma pping Toolbox functions handles raster geodata, see “Understanding Raster Geodata” on page 2-38.
The remainder of this chapter focuses on the fundamental principles of geographic measurement a n d data manipulation tha t are a prerequisite for creating map displays. “Reading and Writing Geospatial Data” on page 2-57
2-11
2 Understanding Map Data
summarizes input functions for importingmanyformatsofgeospatialdata into the toolbox. Chapter 3, “Understanding G eo spatial Geometry” introduces geodetic concepts that underlie all geospatial data and its handling.
2-12

Understanding Vector Geodata

In this section...
“Points, Lines, and Polygons” on page 2-13 “Segments Versus Polygo n s” on page 2-19 “Mapping Toolbox Geographic Data Structures” on page 2-21 “Selecting Data to Read with the shaperead Function” on page 2-32

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 repres ent 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 p airs 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.
Understanding Vector Geodata
1 Display a map of France.
h=wor landa geosh
2 Save the location of the Eiffel Tower in vector form.
Tow Tow
3 Place a red dot on the map to indicate the tower and label it.
ge te
ldmap('France'); reas = shaperead('landareas.shp','UseGeoCoords', true); ow (landareas, 'FaceColor', [1 1 .5]);
erLat = 48.85; erLon = 2.28;
oshow(TowerLat, TowerLon, 'Marker','.','MarkerEdgeColor','red') xtm(TowerLat,TowerLon + 0.5, 'Eiffel Tower');
2-13
2 Understanding Map Data
2-14
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:
Understanding Vector Geodata
h = worldmap('France'); landareas = shaperead('landareas.shp','UseGeoCoords', true); geoshow (landareas, 'FaceColor', [1 1 .5]); rivers = shaperead('worldrivers', 'UseGeoCoords', true); 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:
2-15
2 Understanding Map Data
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'.Datafortheeighthriver, 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.
2-16
1 Copy and paste the following vector variables at the command line:
x=[40553310051040NaN1025302510...
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.
Understanding Vector Geodata
Notice that the NaNs appear in the same location s 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 w ith 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.
2-17
2 Understanding Map Data
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') states = shaperead('usastatelo', 'UseGeoCoords', true); 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')
2-18
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
Understanding Vector Geodata
You can see that Arkansas is of shape type 'Polygon'.Viewthelastentry 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, 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 G reat 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 , When plotted together, the two pairs of variables can form a map. The patch
latcoast and loncoast,that
latborder and lonborder.
2-19
2 Understanding Map Data
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 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 variables that define two polygons each, delineated by NaNs (e.g.,
lonuk).
latengland, lonengland, latscotland,andlonscotland)orintwo
latborder
latuk,
Border line
2-20
Polygon of Great Britain
Polygon and line together
(still one polygon)
The distinction between line and polygon data might not appear to be important,butitcanmakeadifferencewhenyouareperforminggeographic analysis and thematic mapping. For example, polygon data can be treated as line data and displayed w ith functions such as be handled as polygons unless it is restructured to make all objects close on themselves, a s described in “Matching Line Segments” on page 7-4.
linem, but line data cannot
Understanding Vector Geodata
Mapping Toolbox
In examples prov variables. Mapp displaying, ex organized in ge
Ageographicd element per ge and attribut (latitude an coordinates structures (regular or
ided in prior chapters, geodata was in the form of individual
ing Toolbox software also provides an easy means of
tracting, and manipulating collections of vector map features
ographic data structures.
ata structure is a MATLAB structure array that has one
ographic feature. Each feature is represented by coordinates
es. A geographic data structure that holds geographic coordinates
dlongitude)iscalledageostruct, and one that holds map
(projected x and y)iscalledamapstruct. G eographic data
hold only vector features and cannot be used to hold raster data
geolocated data grids or images).
Geographic Data Structures
Shapefiles
Geographi is importe Institut encode co non-geom
Ashapef amainfi base na respec shapef and
c data structures most frequently originate when vector geodata d from a shapefile. The Environmental Systems Research
e designed the shapefile format for vector geodata. Shapefiles
ordinates for points , multipoints, lines, or polygons, along with
etrical attributes.
ile stores attributes and coordinates in separate files; it consists of
le, an index file, and an xBASE file. All three files have the same
me and are distinguished by the extensions
tively. (For example, given the base name
ile filenames would be
ncord_roads.dbf'
'co
'concord_roads.shp', 'concord_roads.shx',
).
.shp, .shx,and.dbf,
'concord_roads' the
The Contents of Geographic Data Structures
haperead
The s
eturns a geographic data structure array. The
and r
rmines the names of the attribute fields at run-time from the shapefile
dete
E table or from optional, user-specified parameters. If a shapefile
xBAS
ribute name cannot be directly used as a field name,
att
field an appropriately modified name, usually by substituting underscores
the
spaces.
for
function reads vector features and attributes from a shapefile
shaperead function
shaperead assigns
2-21
2 Understanding Map Data

Fields in a Geographic Data Structure

Field Name Data Type Description Comments
Geometry
BoundingBox
String One of the following
2-by-2 numerical array
shape types:
'MultiPoint', 'Line',
or
'Polygon'.
Specifies the minimum and maximum feature coordinate values in
'Point',
'PolyLine',
For a the value of the
Geometry field is
simply
'Line'.
Omitted for shape type
'Point'.
eachdimensioninthe following form:
min( ) min( )
XY
⎡ ⎢
max( ) max( )
XY
⎤ ⎥
X, Y, Lon,or Lat
1-by-N array
Coordinate vector.
of class
double
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:
'PolyLineM', 'PolygonZ', 'PolylineM',or'Multipatch'.Also,although 'Null Shape' features can be present in a 'Point', 'Multipoint', 'PolyLine',or'Polygon' shapefile, they are ignored.
'PointZ', 'PointM', 'MultipointZ', 'MultipointM', 'PolyLineZ',
2-22
Understanding Vector Geodata
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 illustra tion of this, see “Displaying a Polygon” on page 2-16. Each disconnected part is separated from the next by a NaN within the You can use the
isShapeMultipart function to determine if a feature has
X and Y (or Lat and Lon)vectors.
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 n ot 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, containing
X and Y fields. This is appropriate if the data set coordinates
shaperead returns a mapstruct
are already projected (in a map coordinate system). O therw ise , if the data set coordinates are unprojected (in a geographic coordinate system), use the parameter-value pair geostruct having
'UseGeoCoords',true to make shaperead return a
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:
Ask your data provider.
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 written in well-know n text, a text m ark-up language. If your contains the term contains the term
.prj file, if one has been provided. The .prj file is
.prj file PROJCS, you have map coordinates. If your .prj file 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
2-23
2 Understanding Map Data
and calls mapshow. If you try to display a geostruct with mapshow,thefunction 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:
S = shaperead('concord_roads.shp')
The output appears as follows:
S= 609x1 struct array with fields:
Geometry BoundingBox X Y STREETNAME RT_NUMBER CLASS ADMIN_TYPE LENGTH
2-24
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,andLENGTH.
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
Understanding Vector Geodata
ADMIN_TYPE: 0
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 have various lengths, but strings, and
CLASS, ADMIN_TYPE and LENGTH must always contain scalar
STREETNAME and RT_NUMBER must always contain
S, X and Y can
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” on page 2-32 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 also might compute vector geodata and attributes by calling various MATL AB 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.
®
Excel®spreadsheet, or by r eading in a delimited text file). Y ou
A one-to-one correspondence e xists between structure elements and geographic features, which extends to the children of hggroups constructed by
mapshow and geoshow.
2-25
2 Understanding Map Data
Achieving these benefits is not dif fi cult. 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 cre ates 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')
2-26
Lat (forpoints,thisisascalardouble)
Lon (forpoints,thisisascalardouble)
% 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
Understanding Vector Geodata
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 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');
2-27
2 Understanding Map Data
2-28
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);
Understanding Vector Geodata
Tracks(3).Name = 'New York-Paris'; [Tracks(3).Lat Tracks(3).Lon] = ...
track2(trackType,latnyc,lonnyc,latparis,lonparis);
5 Compute lengths of the g reat 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]) ...
2-29
2 Understanding Map Data
max([Tracks.Length])],...
'Color',winter(3)});
geoshow(Tracks,'SymbolSpec',colorRange);
2-30
You can save the geostructs you just created as shapefiles by calling
shapewrite with a filename of your choice, for exam p le:
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 a ttribute 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 elemen t 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.
Understanding Vector Geodata
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
mapview. Proper directionality is particularly important if polygons contain
geoshow, mapshow,and
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). W ithin 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 —Convertpolygonalregiontoface-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 f ourth 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, wh en 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 1 display structures and their usage, see “Version 1 Display Structures” in the
updategeostruct. For m ore information about Version
2-31
2 Understanding Map Data
reference page for displaym. Additional information is located in reference pages for
updategeostruct, extractm,andmlayers.

Selecting Data to Read with the shaperead Function

The shaperead function provides you with a powerful method, called a selector,toselectonlythedatafieldsanditemsyouwanttoimportfrom shapefiles.
A selector is a ce ll 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 is not necessarily true: a feature for which the predicate returns
false. The converse
true may be
excluded for other reasons when the selector is used in comb ination 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 the following contents:
function result = roadfilter(roadclass,roadlength) mininumClass = 4; minimumLength = 200; result = (roadclass >= mininumClass) && ...
(roadlength >= minimumLength);
end
roadfilter.m,with
2-32
2 You can then call shaperead like this:
Understanding Vector Geodata
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.
2-33
2 Understanding Map Data
Note that if the call to shaperead took place w ithin a function, then
roadfilter could be defined in a subfunction thereof rather than in a
file of its own.
Example 2: Predicate as Function Handle
As a simple variation on the p revious 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
2-34
Example 3: Predicate as Anonymous Function
Having to define predicate functions in 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', ...
Understanding Vector Geodata
{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
Thereisactuallynoneedtointroduceafunctionhandlevariablewhen 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.
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
2-35
2 Understanding Map Data
Example 5: Parameterizing the Selector; Predicate as Nested Function
In the previous patterns, the predicate invo lves two hard-coded parameters (called the 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
minimumLength at run t ime?
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
minimumClass and minimumLength in roadfilter.m), as well as
roadclass and roadlength input variables. If you use any of these
minimumClass and
nestedroadfilter) and the two attibute names:
2-36
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
Understanding Vector Geodata
2-37
2 Understanding Map Data

Understanding Raster Geodata

In this section...
“Georeferencing Raster Data” on page 2-38 “Regular Data Grids” on page 2-40 “Geolocated Data Grids” on page 2-49

Georeferencing Raster Data

Raster geodata consists of georeferenced data grids and images that in the MATLAB workspace are stored as matrices. While raster geodata looks like any other matrix of real numbers, w hat sets it apart is that it is georeferenced, either to the globe or to a specified map projection, so that each pixel of data occupies a known patch of territory on the planet.
Whether a raster geodata set covers theentireplanetornot,itsplacement and resolution must be specified. Raster geodata is georeferenced in the toolbox throug h a companion data structure called a referencing matrix.This 3-by-2 matrix of doubles describes the scaling, orientation, and placement of the data grid on the globe. For a given referencing matrix, R, one of the following relations holds between rows and columns and coordinates (depending on whether the g rid is based on map coordinates or geographic coordinates, respectively):
2-38
[x y] = [row col 1] * R, or [long lat] = [row col 1] * R
For additional details about and examples of using referencing matrices, see the reference page for
makerefmat.
Referencing Vectors
In many instances (when the data grid or image is based on latitude and longitude and is aligned with the geographic graticule), a referencing matrix has more degrees of freedom than the data requires. In such cases, you can use a more compact representation, a three-element referencing vector.A referencing vector de fine s the pixel siz e and northwest origin for a regular, rectangular data grid:
Understanding Raster Geodata
refvec = [cells-per-degree north-lat west-lon]
In MAT-files, this variable is often called refvec or maplegend.Thefirst element, cells-per-degree, describes the angular extent of each grid cell (e.g., if each cell covers five degrees of latitude and longitude, cells-per-degree would be specified as
0.2). Note that if the latitude extent of cells differs from
their longitude extent you cannot use a referencing vector, and instead must specify a referencing matrix. The second elem ent, north-lat, specifies the northern limit of the data grid (as a latitude), and the third element, west-lon, specifies the western extent of the data grid (as a longitude). In other words, north-lat, west-lon is the northwest corner of the data grid. Note, however, that cell (1,1) is always in the southwest corner of the grid. This need not be the case for grids or images described by referencing matrices, as opposed to referencing vectors.
Note V ersions of Mapping T oolbo x softw are prior to 2.0 did not us e referencing matrices, and called referencing vectors map legend vectors or sometimes just map legends. The current version of the toolbox uses the term legend only to refer to keys to symbolism.
An example of such a grid is the geoid data set (a MAT-file), which represents the shape of the geoid. In the
geoid matrix, each cell represents one degree,
the entire northern edge occupies the north pole, the southern edge occupies the south pole, and the western edge runs down the prime meridian. Thus, the referencing vector for
geoidrefvec = [1 90 0]
geoid is
This structure is stored in the geoid MAT-file (note that it is duplicated by the
geoidlegend referencing vector for backward compatibility). Interpret
this referencing vector as follows:
Each data grid entry represents one degree of latitude and one degree of longitude.
The northern edge of the map is at 90ºN (the North Pole).
The western edge of the map is at 0º (the prime meridian).
2-39
2 Understanding Map Data
All regular data grids require a referencing matrix or vector, even if they cover the entire planet . Geolocated data grids do not, as they explicitly identify the geographic coordinates of all rows and columns. For details on geolocated grids, see “Geolocated Data Grids” on page 2-49. For additional information on referencing matrices and vectors, see the reference pages for
makerefmat, limitm,andsizem.

Regular Data Grids

Regular data grids are rectangular, non-sparse, matrices of class double.
Constructing a Global Data Grid
Imagine an extremely coarse map of the world in which each cell represents 60º. Such a map matrix would be 3-by-6.
1 First create data for this , starting with the data grid:
miniZ = [1 2 3 4 5 6; 7 8 9 10 11 12; 13 14 15 16 17 18];
2-40
2 Now make a referencing matrix:
miniR = makerefmat('RasterSize', size(miniZ), ...
'Latlim', [-90 90], 'Lonlim', [-180 180])
miniR =
060
60 0
-210 -120
3 Set up an equidistant cylindrical map projection:
axesm('MapProjection', 'eqdcylin') setm(gca,'GLineStyle','-', 'Grid','on','Frame','on')
4 Draw a graticule with parallel and meridian labels at 60º intervals:
setm(gca, 'MlabelLocation', 60, 'PlabelLocation',[-30 30],...
'MLabelParallel','north', 'MeridianLabel','on',... 'ParallelLabel','on','MlineLocation',60,... 'PlineLocation',[-30 30])
Understanding Raster Geodata
5 Map the data using meshm and display with a color ramp and legend:
meshm(miniZ, miniR); colormap('autumn'); colorbar
Note that the first row of the matrix is displayed as the bottom of the map, while the last row is displayed as the top.
Computing Map Limits from Referencing Vectors
Given a regular data grid and its referencing vector, the full extent of the g rid can be computed using the a data grid that does not encompass the entire world, do the following exercise:
1 Load the Korea 5-arc-minute elevation grid and inspect the referencing
vector,
refvec:
load korea refvec refvec =
limitm function. To understand how this works for
2-41
2 Understanding Map Data
12 45 115
The refvec referencing vector indicates that there are 12 cells p er angular degree. This horizontal resolution is 5 times finer than that of the
topo
data grid, which is one cell per degree.
2 Use limitm to determine that the korea region extends from 30ºN to 45ºN
and from 115ºW to 135ºW:
[latlimits,longlimits] = limitm(map,refvec)
latlimits =
30 45
longlimits =
115 135
3 Verify this computation manually by getting the dimensions of the
elevation array and computing the eastern and southern map limits from the reference vector:
2-42
[rows cols] = size(map)
rows =
180
cols =
240
southlat = refvec(2) - rows/refvec(1)
southlat =
30
eastlon = refvec(3) + cols/refvec(1)
eastlon =
135
The results match latlimits(1) and longlimits(2).Thetwoformulas use different signs because latitudes decrease southwards and longitudes increase eastward.
Understanding Raster Geodata
Geographic Interpretation of Matrix Elements
You can access and manipulate gridded geodata and its associated referencing vector by either geographic or matrix coordinates. Use the to explore this. As was demonstrated above, the north, south, east, and west limits of the mapped area can be determined as follows:
clear; load russia [latlim,longlim] = limitm(map,refvec)
The limits are:
latlim =
35 80
longlim =
15 190
Display a map of Russia:
worldmap(latlim,longlim); cmap = jet(4); geoshow(map,cmap,maplegend)
russia data set
The data grid in the russia MAT-file ex tends over the international date line (180º longitude). You could use the function
wrapTo180 to rename the
eastern limit to be -170, or 170ºW.
2-43
2 Understanding Map Data
The function setltln retrieves the geographic coordinates of a particular matrix element. The returned coordinates actually show the center of the geographic area represented by the matrix entry:
row = 23; col = 79; [lat,long] = setltln(map,refvec,row,col)
lat =
39.5
long =
30.7
setpostn
does the reverse of this, determining the row and column of the data
grid element containing a given geographic point location:
[r,c] = setpostn(map,maplegend,lat,long)
r=
23
c=
79
The Geography of Gridded Geodata
Each matrix element (analogous to a pixel) can be thought of as a spheroidal quadrangle, which includes its northern and eastern edges, but not its western edge or southern edge.
2-44
Understanding Raster Geodata
East and North edges
N
are included
West and South edges
are excluded
An Element in a Data Grid
The exceptions to this are that the southernmost row (row 1) also contains its southern edge, and the westernmost column (column 1) contains its western edge, except when the map encompasses the entire 360º of longitude. In that case, the westernmost edge of the first column is not included, because it is identical to the easternmost edge of the last column. These exceptions ensure that all points on the globe can be represented exactly once in a regular data grid.
Although each data grid element represents an area, not a point, it is often useful to assign singular coordinates to provide a point of reference. The
setltln function does this. It geolocates an element by the point in the center
of the area represented by the element. The following code references the center cell coordinate for the row 3, column 17 of the Russia map:
clear; load russia row = 3; col = 17; [lat,long] = setltln(map,refvec,row,col)
lat =
35.5
long =
18.3
2-45
2 Understanding Map Data
Because the cells in the russia matrix represent 0.2º squares (5 cells per degree), the cell in question extends from north of 35.4ºS to exactly 35.6ºS, and from e a st of 18.2ºE to e xa ctly 18.4ºE.
Accessing Data Grid Elements
The actual values contained within the map matrix entries are important as well. Several Mapping Toolbox functions can access and alter the values of data grid elements.
If the actual row and column of a desired entry are known, then a simple matrix index can return the appropriate value:
1 Use the row and column from the previous example (row 3, column 17) to
determine the value of t h at cell simply by querying the matrix:
value = map(row,col)
value =
2
2-46
2 More often, the geographic coordinates are known, and the value can be
retrieved with
value = ltln2val(map,maplegend,lat,long)
value =
2
3 The latitude-longitude coordinates associate d with particular v al u es in a
data grid can be found with
find. Here the coordinates of elements in the topo matrix have values
ltln2val:
findm, analogous to the MATLAB function
greater than 5,500 meters:
load topo [lats,longs] = findm(topo>5500,topolegend); [lats longs]
ans =
34.5000 79.5000
34.5000 80.5000
Understanding Raster Geodata
30.5000 84.5000
28.5000 86.5000
4 To get the row and colum n indices instead, simply use the Mapping Toolbox
function
find:
[i,j]=find(topo>5500)
i=
125 125 121 119
j=
80 81 85 87
5 To recode a specific matrix value to some other value, use changem.Loador
reload the
russia M AT-file , and then change all instances of a given value
inadatagridtoanewvalueinonestep:
oldcode = ltln2val(map,maplegend,37,79)
oldcode =
4
newmap = changem(map,5,oldcode); newcode = ltln2val(newmap,maplegend,37,79)
newcode =
5
All entries in newmap corres ponding to 4s in map now have the value 5.
Using a Mask to Recode a Data Grid
You can also define a logical mask to identify the map entries to change. A mask is a matrix the same size as the map matrix, with 1s everywhere that
2-47
2 Understanding Map Data
values are to change. A mask is often generated by a logical operation on a map variable, a topic that is described in greater detail below:
1 The russia data grid contains 3 for each cell covering Russia. To set every
non-Russia m atrix entry to zero, use the following commands:
clear; load russia nonrussia = map; nonrussia(map~=3) = 0;
2 Verify the data that results from these operations:
whos
Name Size Bytes Class clrmap 4x3 96 double description 5x69 690 char map 225x875 1575000 double maplegend 1x3 24 double nonrussia 225x875 1575000 double refvec 1x3 24 double source 1x68 136 char
2-48
newcode = ltln2val(nonrussia,refvec,37,79)
newcode =
0
Precomputing the Size of a Data Grid
Finally, if you know the latitude and longitude li mits of a region, you can calculate the required matrix size and an appropriate referencing vector for any desired map resolution and scale. However, before making a large, memory-taxing data grid, you should first determine what its size will be. For a map of the continental U.S. at a scale of 10 cells per degree, do the following:
1 Compute the matrix dimensions using sizem, specifying latitude limits of
25ºN to 50ºN and longitudes from 60ºW to 130ºW:
cellsperdeg = 10; [r,c,maplegend] = sizem([25 50],[-130 -60],cellsperdeg)
Understanding Raster Geodata
r=
250
c=
700
maplegend =
10 50 -130
msize = r * c * 8
msize =
1400000
This data grid would be 250-by-700, and consume 1,400,000 bytes.
2 Now determine what the storage requirements would be if the scale were
reduced to 5 rows/columns per deg ree:
cellsperdeg2 = 5; [r,c,maplegend] = sizem([25 50],[-130 -60],cellsperdeg2)
r=
125
c=
350
maplegend =
5 50 -130
msize = r * c * 8
msize =
350000
A 125-by-300 matrix that used 350,000 bytes might be more manageable, if it had sufficient resolution at its intended publication scale.

Geolocated Data Grids

In addition to regular data grids, the toolbox provides another format for geodata: geolocated data grids. These multivariate data sets can be displayed, and their values and coordinates can be queried, but unfortunately much of
2-49
2 Understanding Map Data
the functionality supporting regular data grids is not available for geolocated data grids.
The examples thus far have shown maps that covered simple, re gular quadrangles, that is, geographically rectangular and aligned with parallels and meridians. Geolocated data grids, in addition to these rectangular orientations, can have other shapes as well.
Geolocated Grid Format
To define a geolocated data grid, you must define three variables:
A matrix of indices or values associated with the mapped region
A matrix giving cell-by-cell latitude coordinates
A matrix giving cell-by-cell longitude co ordinates
The following exercise demonstrates this data representation:
1 Load the MAT-file example of an irregularly shaped geolocated data grid
called
mapmtx:
2-50
load mapmtx whos
Name Size Bytes Class Attributes
description 1x54 108 char lg1 50x50 20000 double lg2 50x50 20000 double lt1 50x50 20000 double lt2 50x50 20000 double map1 50x50 20000 double map2 50x50 20000 double source 1x43 86 char
Two geolocated data grids are in this data set, each requiring three variables. The values contained in longitude coordinates, respectively, in matrices are the same size. Similarly,
map1 correspond to the latitude and
lt1 an d lg1.Noticethatallthree map2, lt2,andlg2 together form a
second geolocated data grid. These data sets were extracted from the
topo
Understanding Raster Geodata
datagridshowninpreviousexamples. Neither of these maps is regular, because their columns do not run north to south.
2 To see their geography, display t he grids one after another:
close all axesm mercator gridm on framem on h1=surfm(lt1,lg1,map1); h2=surfm(lt2,lg2,map2);
3 Showing coastlines will help to orient you to these skewed grids:
load coast plotm(lat,long,'r')
Notice that neither topo matrix is a regular rectangle. One looks like a diamond geographically, the other lik e a trapezoid. The trapezoid is displayed in two pieces because it crosses the edge of the map. These shapes can be thought of as the geographic organization of the data, just as rectangles are for regular data grids. But, just as for regular data grids,
2-51
2 Understanding Map Data
this organizational logic does not mean that displays of these maps are necessarily a specific shape.
4 Now change the view to a polyconic projection with an origin at 0ºN, 90ºE:
setm(gca,'MapProjection','polycon', 'Origin',[0 90 0])
2-52
As the porti
polyconic projection is limited to a 150º range in longitude, those
ons of the maps outside this region are automatically trimmed.
Geographic Interpretations of Geolocated Grids
ing Toolbox software supports three different interpretations o f
Mapp
located data grids:
geo
Understanding Raster Geodata
First, a map matrix having the same number of rows and columns as the latitude and longitude coordinate matrices represents the val ues of the map d ata at the corresponding geographic points (centers of data cells).
Next, a map matrix having one fewer row and one fewer column than the geographic coordinate matrices represents the values of the map data within the area formed by the four adjacent latitudes and longitudes.
Finally, if the latitude and longitude matrices have smaller dimensions than the map matrix, you can interpret them as describing a coarser graticule, or mesh of latitude and longitude cells, into which the blocks of map data are warped.
This section discusses the first two interpretations of geolocated data grids. For more information on the use of graticules, see “The Map Grid” on page 4-55.
Type 1: Values associated with upper left grid coordinate. As an example of the first interpretation, considera4-by-4mapmatrixwhosecell size is 30-by-30 degrees, along with its corresponding 4-by-4 latitude and longitude matrices:
map = [ 1 2 3 4;...
5 6 7 8;... 9 10 11 12;... 3141516];
lat = [ 30 30 30 30;...
0 0 0 0;...
-30 -30 -30 -30;...
-60 -60 -60 -60];
long = [0 30 60 90;...
0 30 60 90;... 0 30 60 90;... 0306090];
This geolocated data grid is displayed with the values of map shown at the associated latitudes and longitudes.
2-53
2 Understanding Map Data
Notice that each cell is are given b last row an in the
Type 2: Va
second in
long var
map = [1 2 3;...
4 5 6;... 7 8 9];
Here is
r of the associated cells:
cente
only 9 of the 16 total cells are displayed. The value displayed for
the value at the upper left corner of that cell, whose coordinates y the corresponding d column of the map matrix are not displayed, although they exist
ta
CDa
property of the surface object.
lues centered within four adjacent coordinates. For the
terpretation, consider a 3-by-3 map matrix with the same
iables:
a surface plot of the map matrix, with the values of
lat and long elements. By convention, the
lat and
map shown at the
2-54
Understanding Raster Geodata
All the map data is displayed for this geolocated data grid. The value of each cell is the value at the center of the cell, and the latitudes and longitudes in the coordinate matrices are the boundaries for the cells.
Ordering of Cells. You may have noticed that the first row of the matrix is displayed as the top of the map, whereas for a regular data grid, the opposite was true: the f irst row corresponded to the bottom of the map. This difference is entirely due to how the
lat and long matrices are ordered. In a geolocated
data grid, the order of values in the two coordinate matrices determines the arrangement of the displayed values.
Transforming Regular to Geolocated Grids. When required, a regular data grid can be transformed into a geolocated data grid. This simply requires that a pair of coordinates matrices be computed at the desired spatial resolution from the regular grid. Do this with the
meshgrat function, a s
follows:
load topo [lat,lon] = meshgrat(topo,topolegend);
Name Size Bytes Class Attributes
lat 180x360 518400 double lon 180x360 518400 double topo 180x360 518400 double topolegend 1x3 24 double topomap1 64x3 1536 double topomap2 128x3 3072 double
2-55
2 Understanding Map Data
Transforming Geolocated to Regular Grids. Conversely, a regular data
grid can also be constructed from a geolocated data grid. The coordinates and values can be embedded in a new regular data grid. The function that performs this conversion is a cell size as inputs.
geoloc2grid;ittakesageolocateddatagridand
2-56
Loading...