Understanding Raster Geodata :: Understanding Map Data (Mapping Toolbox™)

Understanding Raster Geodata

On this page…
Georeferencing Raster Data Regular Data Grids Geolocated Data Grids

On this page…

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, what 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 the entire planet or not, its placement and resolution must be specified. Raster geodata is georeferenced in the toolbox through 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 grid is based on map coordinates or geographic coordinates, respectively):

[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 defines the pixel size and northwest origin for a regular, rectangular data grid:

refvec = [cells-per-degree north-lat west-lon]

In MAT-files, this variable is often called refvec or maplegend. The first 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 element, 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 Versions of Mapping Toolbox software prior to 2.0 did not use 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 geoid is

geoidrefvec = [1 90 0]

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

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. For additional information on referencing matrices and vectors, see the reference pages for makerefmat, limitm, and sizem.

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.

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

Now make a referencing matrix:

miniR = makerefmat('RasterSize', size(miniZ), ...
   'Latlim', [-90 90], 'Lonlim', [-180 180])

miniR =

     0    60
    60     0
  -210  -120

Set up an equidistant cylindrical map projection:

axesm('MapProjection', 'eqdcylin')
setm(gca,'GLineStyle','-', 'Grid','on','Frame','on')

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

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 grid can be computed using the limitm function. To understand how this works for a data grid that does not encompass the entire world, do the following exercise:

Load the Korea 5-arc-minute elevation grid and inspect the referencing vector, refvec:
```
load korea
refvec
refvec =
           12           45          115
```
The refvec referencing vector indicates that there are 12 cells per angular degree. This horizontal resolution is 5 times finer than that of the topo data grid, which is one cell per degree.
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
```

Verify this computation manually by getting the dimensions of the elevation array and computing the eastern and southern map limits from the reference vector:

[rows cols] = size(map)

rows =
   180
cols =
   240

southlat = refvec(2) - rows/refvec(1)

southlat =
    30

eastlon = refvec(3) + cols/refvec(1)

eastlon =
   135

The results match latlimits(1) and longlimits(2). The two formulas use different signs because latitudes decrease southwards and longitudes increase eastward.

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 russia data set 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)

latlim =
	35    80
longlim =
	15   190

The data grid in the russia MAT-file extends over the international date line (180º longitude). You could use the function wrapTo180 to rename the eastern limit to be -170, or 170ºW.

The function setltln retrieves the geographic coordinates of a particular matrix element. The returned coordinates actually show the center of the geographic area represented by the matrix entry:

row = 23; col = 79;
[lat,long] = setltln(map,refvec,row,col)

lat =
	39.5
long =
	30.7

setpostn does the reverse of this, determining the row and column of the data grid element containing a given geographic point location:

[r,c] = setpostn(map,maplegend,lat,long)

r =
	23
c =
	79

The Geography of Gridded Geodata

Each matrix element (analogous to a pixel) can be thought of as a spheroidal quadrangle, which includes its northern and eastern edges, but not its western edge or southern edge.

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

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 east of 18.2ºE to exactly 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:

Use the row and column from the previous example (row 3, column 17) to determine the value of that cell simply by querying the matrix:
```
value = map(row,col)

value =
	2
```
More often, the geographic coordinates are known, and the value can be retrieved with ltln2val:
```
value = ltln2val(map,maplegend,lat,long)

value =
	2
```
The latitude-longitude coordinates associated with particular values in a data grid can be found with findm, analogous to the MATLAB function find. Here the coordinates of elements in the topo matrix have values greater than 5,500 meters:
```
load topo
[lats,longs] = findm(topo>5500,topolegend);
[lats longs]

ans =
	34.5000   79.5000
	34.5000   80.5000
	30.5000   84.5000
	28.5000   86.5000
```
To get the row and column indices instead, simply use the Mapping Toolbox function find:
```
[i,j]=find(topo>5500)

i =
   125
   125
   121
   119
j =
    80
    81
    85
    87
```
To recode a specific matrix value to some other value, use changem. Load or reload the russia MAT-file, and then change all instances of a given value in a data grid to a new value in one step:
```
oldcode = ltln2val(map,maplegend,37,79)

oldcode =
	4

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

newcode =
	5
```
All entries in newmap corresponding 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 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:

The russia data grid contains 3 for each cell covering Russia. To set every non-Russia matrix entry to zero, use the following commands:
```
clear; load russia
nonrussia = map;
nonrussia(map~=3) = 0;
```

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

newcode = ltln2val(nonrussia,refvec,37,79)

newcode =
	0

Precomputing the Size of a Data Grid

Finally, if you know the latitude and longitude limits 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:

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)

r =
   250
c =
   700
maplegend =
    10    50  -130

msize = r * c * 8

msize =
     1400000
```
This data grid would be 250-by-700, and consume 1,400,000 bytes.

Now determine what the storage requirements would be if the scale were reduced to 5 rows/columns per degree:

cellsperdeg2 = 5;
[r,c,maplegend] = sizem([25 50],[-130 -60],cellsperdeg2)

r =
   125
c =
   350
maplegend =
     5    50  -130

msize = r * c * 8

msize =
      350000

A 125-by-300 matrix that used 350,000 bytes might be more manageable, if it had sufficient resolution at its intended publication scale.

Geolocated Data Grids

In addition to regular data grids, the toolbox provides another format for geodata: geolocated data grids. These multivariate data sets can be displayed, and their values and coordinates can be queried, but unfortunately much of the functionality supporting regular data grids is not available for geolocated data grids.

The examples thus far have shown maps that covered simple, regular 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 coordinates

The following exercise demonstrates this data representation:

Load the MAT-file example of an irregularly shaped geolocated data grid called mapmtx:

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 map1 correspond to the latitude and longitude coordinates, respectively, in lt1 and lg1. Notice that all three matrices are the same size. Similarly, map2, lt2, and lg2 together form a second geolocated data grid. These data sets were extracted from the topo data grid shown in previous examples. Neither of these maps is regular, because their columns do not run north to south.

To see their geography, display the grids one after another:

close all
axesm mercator
gridm on
framem on
h1=surfm(lt1,lg1,map1);
h2=surfm(lt2,lg2,map2);

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 like 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, this organizational logic does not mean that displays of these maps are necessarily a specific shape.
Now change the view to a polyconic projection with an origin at 0ºN, 90ºE:
```
setm(gca,'MapProjection','polycon', 'Origin',[0 90 0])
```

As the polyconic projection is limited to a 150º range in longitude, those portions of the maps outside this region are automatically trimmed.

Geographic Interpretations of Geolocated Grids

Mapping Toolbox software supports three different interpretations of geolocated data grids:

First, a map matrix having the same number of rows and columns as the latitude and longitude coordinate matrices represents the values of the map data 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.

Type 1: Values associated with upper left grid coordinate. As an example of the first interpretation, consider a 4-by-4 map matrix whose cell 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;...
	3 14 15 16];

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;...
	0 30 60 90];

This geolocated data grid is displayed with the values of map shown at the associated latitudes and longitudes.

Notice that only 9 of the 16 total cells are displayed. The value displayed for each cell is the value at the upper left corner of that cell, whose coordinates are given by the corresponding lat and long elements. By convention, the last row and column of the map matrix are not displayed, although they exist in the CData property of the surface object.

Type 2: Values centered within four adjacent coordinates. For the second interpretation, consider a 3-by-3 map matrix with the same lat and long variables:

map = [1 2 3;...
		4 5 6;...
		7 8 9];

Here is a surface plot of the map matrix, with the values of map shown at the center of the associated cells:

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 first 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, as 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

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 geoloc2grid; it takes a geolocated data grid and a cell size as inputs.

Understanding Vector Geodata

Reading and Writing Geospatial Data

Recommended Products

Get MATLAB trial software

Includes the most popular MATLAB recorded presentations with Q&A sessions led by MATLAB experts.

Get the Interactive Kit