Quantcast

Documentation Center

  • Trial Software
  • Product Updates

Manipulating Raster Geodata

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.

Creating Data Grids from Vector Data

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

  1. 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;
  2. 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.

  3. Make a map of the data grid in contrasting colors:

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

  4. Set up the map limits:

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

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

  6. Clear and redraw the map using the filled grid:

    meshm(inGrid3, inRefVec)

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

Using a GUI to Rasterize Polygons

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:

  1. 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];
  2. 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);
  3. 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)

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

  5. Use encodem to fill each country with a unique value, producing a new grid:

    meGrid5 = encodem(meGrid, [row col val], 1);
  6. 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.

Data Grids as Logical Variables

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.

Obtaining the Area Occupied by a Logical Grid Variable

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:

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

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

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

Data Grid Values Along a Path

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.

Using the mapprofile Function

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

  1. Load the Korean elevation data:

    figure;
    load korea
  2. 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.

  3. Render the map and apply a digital elevation model (DEM) colormap to it:

    meshm(map,maplegend,size(map),map)
    demcmap(map)
  4. Define endpoints for a straight-line transect through the region:

    plat = [40.5 30.7];
    plon = [121.5 133.5];
  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);
  6. Draw the transect in 3-D so it follows the terrain:

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

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

Data Grid Gradient, Slope, and Aspect

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.

Computing Gradient Data from a Regular Data Grid

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

  1. 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));
  2. 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
    
  3. 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
  4. Clear the frame and display the slope grid:

    figure; axesm eqacyl
    meshm(slope,R)
    colormap (jet(64))
    colorbar('vert');
    title('Peaks: slope')
  5. Map the aspect grid:

    figure; axesm eqacyl
    meshm(aspect,R)
    colormap (jet(64))
    colorbar('vert');
    title('Peaks: aspect')
  6. Map the gradients to the north:

    figure; axesm eqacyl
    meshm(gradN,R)
    colormap (jet(64))
    colorbar('vert');
    title('Peaks: North gradient')
  7. 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?