## Documentation Center |

On this page… |
---|

Vector-to-Raster Data Conversion |

You can convert latitude-longitude vector data to a grid at
any resolution you choose to make a raster base map or grid layer.
Certain Mapping Toolbox™ GUI tools help you do some of this, but
you can also perform vector-to-raster conversions from the command
line. The principal function for gridding vector data is `vec2mtx`,
which allocates lines to a grid of any size you indicate, marking
the lines with 1s and the unoccupied grid cells with 0s. The grid
contains doubles, but if you want a logical grid (see Data Grids as Logical Variables, below) cast the result to
be a logical array.

If the vector data consists of polygons (patches), the gridded
outlines are all hollow. You can differentiate them using the `encodem` function,
calling it with an array of rows, columns, and seed values to produce
a new grid containing polygonal areas filled with the seed values
to replace the binary values generated by `vec2mtx`.

To demonstrate vector-to-raster data conversion, use patch data
for Indiana from the `usastatehi` shapefile:

Use

`shaperead`to get the patch data for the boundary:indiana = shaperead('usastatehi.shp',... 'UseGeoCoords', true,... 'Selector', {@(name)strcmpi('Indiana',name), 'Name'}); inLat = indiana.Lat; inLon = indiana.Lon;

Set the grid density to be 40 cells per degree, and use

`vec2mtx`to rasterize the boundary and generate a referencing vector for it:gridDensity = 40; [inGrid, inRefVec] = vec2mtx(inLat, inLon, gridDensity); whos

Name Size Bytes Class gridDensity 1x1 8 double inGrid 164x137 179744 double inLat 1x626 5008 double inLon 1x626 5008 double inRefVec 1x3 24 double indiana 1x1 10960 struct

The resulting grid contains doubles, and has 80 rows and 186 columns.

Make a map of the data grid in contrasting colors:

figure axesm eqdcyl meshm(inGrid, inRefVec) colormap jet(4)

[latlim, lonlim] = limitm(inGrid, inRefVec); setm(gca, 'Flatlimit', latlim, 'FlonLimit', lonlim) tightmap

To fill (recode) the interior of Indiana, you need a seed point (which must be identified by row and column) and a seed value (to be allocated to all cells within the polygon). Select the middle row and column of the grid and choose an index value of

`3`to identify the territory when calling encodem to generate a new grid:inPt = round([size(inGrid)/2, 3]); inGrid3 = encodem(inGrid, inPt,1);

The last argument (

`1`) identifies the code for boundary cells, where filling should halt.Clear and redraw the map using the filled grid:

meshm(inGrid3, inRefVec)

Plot the original vectors on the grid to see how well data was rasterized:

`plotm(inLat, inLon,'k')`

The resulting map is shown on the left below. Use the

**Zoom**tool on the figure window to examine the gridding results more closely, as the right-hand figure shows.

See the `vec2mtx` and `encodem` reference pages for further information. `imbedm` is a related function for gridding
point values.

In the previous example, had you wanted to include the states
that border Indiana, you could also have extracted Illinois, Kentucky,
Ohio, and Michigan, and then deleted unwanted areas of these polygons
using `maptrimp` (see Trimming Vector Data to a Rectangular Region for specific
details on its use). Use the `seedm` function with
seed points found using the `getseeds` GUI to fill
multiple polygons after they are gridded:

Extract the data for Indiana and its neighbors by passing their names in a cell array to

`shaperead`:pcs = {'Indiana', 'Michigan', 'Ohio', 'Kentucky', 'Illinois'}; centralUS = shaperead('usastatelo.shp',... 'UseGeoCoords', true,... 'Selector',{@(name)any(strcmpi(name,pcs),2), 'Name'}); meLat = [centralUS.Lat]; meLon = [centralUS.Lon];

Rasterize the trimmed polygons at a 1-arc-minute resolution (60 cells per degree), also producing a referencing vector:

[meGrid, meRefVec] = vec2mtx(meLat, meLon, 60);

Set up a map figure and display the binary grid just created:

figure axesm eqdcyl geoshow(meLat, meLon, 'Color', 'r'); meshm(meGrid, meRefVec) colormap jet(8)

Use

`getseeds`to interactively pick seed points for Indiana, Michigan, Ohio, Kentucky, and Illinois, in that order:[row,col,val] = getseeds(meGrid, meRefVec, 5, [3 4 5 6 7]); [row col val]

ans = 207 140 3 219 326 4 212 506 5 56 459 6 393 433 7

The MATLAB

^{®}prompt returns after you pick five locations in the figure window. As you chose them yourself, your`row`and`col`numbers will differ.Use

`encodem`to fill each country with a unique value, producing a new grid:meGrid5 = encodem(meGrid, [row col val], 1);

Clear the display and display

`meGrid5`to see the result:clma meshm(meGrid5, meRefVec)

The rasterized map of Indiana and its neighbors is shown below.

See the `getseeds` reference
page for more information. The `maptrim` and `seedm` GUI tools are also useful in this
context.

You can apply logical criteria to numeric data grids to create *logical
grids*. Logical grids are data grids consisting entirely
of 1s and 0s. You can create them by performing logical tests on data
grid variables. The resulting binary grid is the same size as the
original grid(s) and can use the same referencing vector, as the following
hypothetical data operation illustrates:

logicalgrid = (realgrid > 0);

This transforms all values greater than 0 into 1s and all other values to 0s. You can apply multiple conditions to a grid in one operation:

logicalgrid = (realgrid >- 100)&(realgrid < 100);

If several grids are the same size and share the same referencing vector (i.e., the grids are co-registered), you can create a logical grid by testing joint conditions, treating the individual data grids as map layers:

logicalgrid = (population > 10000)&(elevation < 400)&... (country == nigeria);

Several Mapping Toolbox functions enable the creation of
logical grids using logical and relational operators. Grids resulting
from such operations contain logical rather than numeric values (which
reduce storage by a factor of 8), but might need to be cast to `double` in
order to be used in certain functions. Use the `onem` and `zerom` functions to create grids of all 1s
and all 0s.

You can analyze the results of logical grid manipulations to
determine the area satisfying one or more conditions (either coded
as 1s or an expression that yields a logical value of 1). The `areamat` function
can provide the fractional surface area on the globe associated with
1s in a logical grid. Each grid element is a quadrangle, and the sum
of the areas meeting the logical condition provides the total area:

You can use the topo grid and the greater-than relational operator to determine what fraction of the Earth lies above sea level:

load topo topoR = makerefmat('RasterSize', size(topo), ... 'Latlim', [-90 90], 'Lonlim', [0 360]); a = areamat((topo > 0),topoR)

a = 0.2890

The answer is about 30%. (Note that land areas below sea level are excluded.)

You can include a planetary radius in specified units if you want the result to have those units. Here is the same query specifying units of square kilometers:

`a = areamat((topo > 0),topoR,earthRadius('km'))`

a = 1.4739e+08

Use the

`usamtx`data grid codes to find the area of a specific state within the U.S.A. As an example, determine the area of the state of Texas, which is coded as 46 in the`usamtx`grid:load usamtx a = areamat((map == 46), refvec, earthRadius('km'))

a = 6.2528e+005

The grid codes 625,277 square kilometers of land area as belonging to the U.S.

You can construct more complex queries. For instance, using the last example, compute what portion of the land area of the conterminous U.S. that Texas occupies (water and bordering countries are coded with 2 and 3, respectively):

usaland = areamat((map > 3 | map == 1), maplegend); texasland = areamat((map == 46), maplegend); texasratio = texasland/usaland

texasratio = 0.0735

This indicates that Texas occupies roughly 7.35% of the land area of the U.S.

For further information, see the `areamat` reference
page.

A common application for gridded geodata is to calculate data
values along a path, for example, the computation of terrain height
along a transect, a road, or a flight path. The `mapprofile` function
does this, based on numerical data defining a set of waypoints, or
by defining them interactively via graphic input from a map display.
Values computed for the resulting profile can be displayed in a new
plot or returned as output arguments for further analysis or display.

The following example computes the elevation profile along a straight line:

Load the Korean elevation data:

`figure; load korea`

Get its latitude and longitude limits using

`limitm`and use them to set up a map frame via`worldmap`:[latlim, lonlim] = limitm(map, maplegend); worldmap(latlim, lonlim)

`worldmap`plots only the map frame.Render the map and apply a digital elevation model (DEM) colormap to it:

meshm(map,maplegend,size(map),map) demcmap(map)

Define endpoints for a straight-line transect through the region:

plat = [40.5 30.7]; plon = [121.5 133.5];

Compute the elevation profile, defaulting the track type to

`great circle`and the interpolation type to`bilinear`:[z,rng,lat,lon] = mapprofile(map,maplegend,plat,plon);

Draw the transect in 3-D so it follows the terrain:

plot3m(lat,lon,z,'w','LineWidth',2)

Construct a plot of transect elevation and range:

`figure; plot(rng,z,'r')`

The `mapprofile` function has other useful
options, including the ability to interactively define tracks and
specify units of distance for them. For further information, see the `mapprofile` reference page.

A map profile is often used to determine slopes along a path.
A related application is the calculation of slope at all points on
a matrix. The `gradientm` function uses a finite-difference
approach to compute gradients for either a regular or a georeferenced
data grid. The function returns the components of the gradient in
the north and east directions (i.e., north-to-south, east-to-west),
as well as slope and aspect. The *gradient* components
are the change in the grid variable per meter of distance in the north
and east directions. If the grid contains elevations in meters, the *aspect* and *slope* are
the angles of the surface normal clockwise from north and up from
the horizontal. Slope is defined as the change in elevation per unit
distance along the path of steepest ascent or descent from a grid
cell to one of its eight immediate neighbors, expressed as the arctangent.
The angles are in units of degrees by default.

The following example illustrates computation of gradient, slope,
and aspect data grids for a regular data grid based on the MATLAB `peaks` function:

Construct a 100-by-100 grid using the

`peaks`function and construct a referencing matrix for it:`datagrid = 500*peaks(100); R = makerefmat('RasterSize',size(datagrid));`

Use

`gradientm`to generate grids containing aspect, slope, gradients to north, and gradients to east:[aspect,slope,gradN,gradE] = gradientm(datagrid,R); whos

Name Size Bytes Class aspect 100x100 80000 double datagrid 100x100 80000 double gradE 100x100 80000 double gradN 100x100 80000 double gridrv 1x3 24 double slope 100x100 80000 double

Map the surface data in a cylindrical equal area projection. Start with the original elevations:

figure; axesm eqacyl meshm(datagrid,R) colormap (jet(64)) colorbar('vert') title('Peaks: elevation') axis square

Clear the frame and display the

`slope`grid:figure; axesm eqacyl meshm(slope,R) colormap (jet(64)) colorbar('vert'); title('Peaks: slope')

figure; axesm eqacyl meshm(aspect,R) colormap (jet(64)) colorbar('vert'); title('Peaks: aspect')

Map the gradients to the north:

figure; axesm eqacyl meshm(gradN,R) colormap (jet(64)) colorbar('vert'); title('Peaks: North gradient')

Finally, map the gradients to the east:

figure; axesm eqacyl meshm(gradE,R) colormap (jet(64)) colorbar('vert'); title('Peaks: East gradient')

The maps of the `peaks` surface elevation and
gradient data are shown below. See the `gradientm` reference
page for additional information.

Was this topic helpful?