Documentation |
On this page… |
---|
Raster geodata consists of georeferenced data grids and images that in the MATLAB^{®} workspace are stored as matrices or objects. 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. This additional information can be supplied in the form of a referencing object, a referencing matrix, or a referencing vector.
A referencing object is an instance of the map.rasterref.GeographicRasterReference class, for raster data referenced to a geographic latitude-longitude system, or the map.rasterref.MapRasterReference class, for raster data referenced to a planar (projected) map coordinate system. A spatial referencing object encapsulates the relationship between a geographic or planar coordinate system and a system of intrinsic coordinates anchored to the columns and rows of a 2-D spatially referenced raster grid or image. Unlike the older referencing matrix and vector representations (described below), a referencing object is self-documenting, providing a rich set of properties to describe both the intrinsic and extrinsic geometry. The use of referencing objects is preferred, but referencing matrices and vectors continue to be supported for the purpose of compatibility.
A referencing matrix is a 3-by-2 matrix of doubles that 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.
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 may encounter 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 object or 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 objects or matrices.
All regular data grids require a referencing object, 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.
Regular data grids are rectangular, non-sparse, matrices of class double.
Imagine an extremely coarse map of the world in which each cell represents 60º. Such a map matrix would be 3-by-6.
miniZ = [1 2 3 4 5 6; 7 8 9 10 11 12; 13 14 15 16 17 18];
Now make a referencing object:
miniR = georasterref('RasterSize', size(miniZ), ... 'Latlim', [-90 90], 'Lonlim', [-180 180])
Your output appears like this:
miniR = GeographicCellsReference with properties: LatitudeLimits: [-90 90] LongitudeLimits: [-180 180] RasterSize: [3 6] RasterInterpretation: 'cells' ColumnsStartFrom: 'south' RowsStartFrom: 'west' CellExtentInLatitude: 60 CellExtentInLongitude: 60 RasterExtentInLatitude: 180 RasterExtentInLongitude: 360 XIntrinsicLimits: [0.5 6.5] YIntrinsicLimits: [0.5 3.5] CoordinateSystemType: 'geographic' AngleUnit: 'degree'
Set up an equidistant cylindrical map projection:
figure('Color','white') ax = axesm('MapProjection', 'eqdcylin'); axis off setm(ax,'GLineStyle','-', 'Grid','on','Frame','on')
Draw a graticule with parallel and meridian labels at 60º intervals:
setm(ax, 'MlabelLocation', 60, 'PlabelLocation',[-30 30],... 'MLabelParallel','north', 'MeridianLabel','on',... 'ParallelLabel','on','MlineLocation',60,... 'PlineLocation',[-30 30])
Map the data using geoshow and display with a color ramp and legend:
geoshow(miniZ, miniR, 'DisplayType', 'texturemap'); colormap('autumn') colorbar
Note that the first row of the matrix is displayed at the bottom of the map, while the last row is displayed at the top.
The latitude and longitude limits of a regular grid are among the most important properties of its referencing object. Given an older data grid that includes a referencing vector but not a referencing object, a good way to check the limits is to convert the referencing vector to a geographic raster reference object, as in the following simple exercise:
Load the Korea 5-arc-minute elevation grid and inspect the referencing vector, refvec:
korea = load('korea','map','refvec')
korea = map: [180x240 double] refvec: [12 45 115]
The referencing vector, korea.refvec, indicates that there are 12 cells per angular degree. This horizontal resolution is 5 times finer than that of the topo grid, which is one cell per degree. It also indicates that the northwest corner is at 45 degrees North, 115 degrees East, but it needs to be combined with the size of the data grid before the southern and eastern limits can be determined.
Convert the referencing vector to a geographic raster reference object (providing a size vector):
korea.R = refvecToGeoRasterReference(korea.refvec, ...
size(korea.map))
korea = map: [180x240 double] refvec: [12 45 115] R: [1x1 map.rasterref.GeographicCellsReference]
Examine the Latlim and Lonlim properties:
korea.R.LatitudeLimits korea.R.LongitudeLimits
ans = 30 45 ans = 115 135
You can access and manipulate gridded geodata using either geographic or intrinsic raster coordinates. Use the russia.mat file to explore this. As was demonstrated above, the north, south, east, and west limits of the mapped area can be determined as follows:
russia = load('russia','map','refvec'); R = refvecToGeoRasterReference(russia.refvec, size(russia.map)); R.LatitudeLimits R.LongitudeLimits
ans = 35 80 ans = 15 190
Display a map of Russia:
figure('Color','white') worldmap(R.LatitudeLimits,R.LongitudeLimits) cmap = jet(4); geoshow(russia.map,cmap,R)
The map.rasterref.GeographicCellsReference.intrinsicToGeographic method can be used to retrieve the geographic coordinates at the center of a given grid cell. For example, consider the cell in row 23, column 79. In intrinsic raster coordinates, the center of this cell is located at:
xIntrinsic = 79; yIntrinsic = 23;
This corresponds to the following location in latitude-longitude, obtained via the intrinsicToGeographic method:
[lat, lon] = intrinsicToGeographic(R, xIntrinsic, yIntrinsic)
Your output appears like this:
lat = 39.5000 lon = 30.7000
The geographicToIntrinsic method does the reverse, converting from latitude-longitude to intrinsic x and y:
[xIntrinsic, yIntrinsic] = geographicToIntrinsic(R, lat, lon)
Your output appears as follows:
xIntrinsic = 79 yIntrinsic = 23
Finally, if you know the latitude and longitude limits of a region, you can calculate the required matrix size and an appropriate referencing object 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:
Specify latitude limits of 25ºN to 50ºN and longitudes from 60ºW to 130ºW:
latlim = [ 25 50]; lonlim = [-130 -60];
Specify 10 cells per degree, and compute the number of rows and columns that this cell density implies:
cellsPerDegree = 10; numRows = cellsPerDegree * diff(latlim); numCols = cellsPerDegree * diff(lonlim);
Construct a referencing object and verify that its size is a reasonable 250 by 700 cells:
R = georasterref('RasterSize',[numRows numCols], ... 'Latlim', latlim, 'Lonlim', lonlim); R.RasterSize
Your output appears like this:
ans = 250 700
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.
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.
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 the 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:
Z = [ ... 1 2 3 4; ... 5 6 7 8; ... 9 10 11 12; ... 13 14 15 16]; lat = [ ... 30 30 30 30; ... 0 0 0 0; ... -30 -30 -30 -30; ... -60 -60 -60 -60]; lon = [ ... 0 30 60 90;... 0 30 60 90;... 0 30 60 90;... 0 30 60 90];
Display the geolocated data grid with the values of map shown at the associated latitudes and longitudes:
figure('Color','white','Colormap',autumn(64)) set(gcf,'Renderer','zbuffer') axesm('pcarree','Grid','on','Frame','on',... 'PLineLocation',30,'PLabelLocation',30) mlabel; plabel; axis off; tightmap h = geoshow(lat,lon,Z,'DisplayType','surface'); set(h,'ZData',zeros(size(Z))) ht = textm(lat(:),lon(:),num2str(Z(:)), ... 'Color','blue','FontSize',14); colorbar('southoutside')
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 lon 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 lon variables:
delete(h) delete(ht) Z3by3 = [ ... 1 2 3; ... 4 5 6; ... 7 8 9]; h = geoshow(lat,lon,Z3by3,'DisplayType','texturemap'); tlat = [ ... 15 15 15; ... -15 -15 -15; ... -45 -45 -45]; tlon = [ ... 15 45 75; ... 15 45 75; ... 15 45 74]; textm(tlat(:),tlon(:),num2str(Z3by3(:)), ... 'Color','blue','FontSize',14)
Display 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 lon 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.