Mathworks MAPPING TOOLBOX 3 user guide

Download

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

suggest@mathworks.com Product enhancement suggestions

bugs@mathwo doc@mathworks.com Documentation error reports service@mathworks.com Order status, license renewals, passcodes

info@mathwo

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

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

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

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

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

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

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

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

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

Manipulating Geospatial Data

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

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

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

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

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

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

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-3

........................ B-3

................... B-4

.............. B-2

................ B-3

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-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 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 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',

'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( )

⎡ ⎢

max( ) max( )

⎣

⎤ ⎥

⎦

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

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 =

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)

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

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-46

2 More often, the geographic coordinates are known, and the value can be

retrieved with

value = ltln2val(map,maplegend,lat,long)

value =

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)

125 125 121 119

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 =

newmap = changem(map,5,oldcode); newcode = ltln2val(newmap,maplegend,37,79)

newcode =

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 =

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

250

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)

125

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

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

Mathworks MAPPING TOOLBOX 3 user guide

Specifications and Main Features

Frequently Asked Questions

User Manual

Getting Started

Product Overview

Dedication and Acknowledgment

You r Fir st Map s

See the World

Tour Boston with the Map Viewer

A Map Viewer Session

Getting More Help

Ways to Get Mapping Toolbox Help

Consulting Release Notes

Mapping Toolbox Demos and Data

Available Demos

Creating Map Displays

Working with Georeferenced Images

Interacting with Web Map Service (WMS) Servers

Locating Geospatial Data

Understanding Map Data

Maps and Map Data

What Is a Map?

What Is Geospatial Data?

Types of Map Data Handled by the Toolbox

Vector Geodata

ALookatVectorData

Raster Geodata

Digital Elevation Data

Remotely Sensed Image Data

A Look at Raster Data

Viewing Raster and Vector Data on the Same Map

Understanding Vector Geodata

Points, Lines, and Polygons

Displaying a Point

Displaying a Line

Displaying a Polygon

Segments Versus Polygons

Shapefiles

The Contents of Geographic Data Structures

Fields in a Geographic Data Structure

Examining a Geographic Data Structure

How to Construct Geographic Data Structures

Mapping Toolbox Version 1 Display Structures

Selecting Data to Read with the shaperead Function

Example 1: Predicate Function in Separate File

Example 2: Predicate as Function Handle

Example 3: Predicate as Anonymous Function

Understanding Raster Geodata

Georeferencing Raster Data

Referencing Vectors

Regular Data Grids

Constructing a Global Data Grid

Computing Map Limits from Referencing Vectors

Geographic Interpretation of Matrix Elements

The Geography of Gridded Geodata

Accessing Data Grid Elements

Using a Mask to Recode a Data Grid

Precomputing the Size of a Data Grid

Geolocated Data Grids

Geolocated Grid Format

Geographic Interpretations of Geolocated Grids