Draping Data on Elevation Maps :: Making Three-Dimensional Maps (Mapping Toolbox™)

Draping Data on Elevation Maps

On this page…
Draping Geoid Heights over Topography Draping Data over Terrain with Different Gridding

On this page…

Draping Data over Terrain with Different Gridding

Draping Geoid Heights over Topography

Lighting effects can provide important visual cues when elevation maps are combined with other kinds of data. The shading resulting from lighting a surface makes it possible to "drape" satellite data over a grid of elevations. It is common to use this kind of display to overlay georeferenced land cover images from Earth satellites such as LANDSAT and SPOT on topography from digital elevation models. Mapping Toolbox displays use variations of techniques described in the previous section.

When the elevation and image data grids correspond pixel-for-pixel to the same geographic locations, you can build up such displays using the optional altitude arguments in the surface display functions. If they do not, you can interpolate one or both source grids to a common mesh.

The following example shows the figure of the Earth (the geoid data set) draped on topographic relief (the topo data set). The geoid data is shown as an attribute (using a color scale) rather than being depicted as a 3-D surface itself. The two data sets are both 1-by-1-degree meshes sharing a common origin.

Note The geoid can be described as the surface of the ocean in the absence of waves, tides, or land obstructions. It is influenced by the gravitational attraction of denser or lighter materials in the Earth's crust and interior and by the shape of the crust. A model of the geoid is required for converting ellipsoidal heights (such as might be obtained from GPS measurements) to orthometric heights. Geoid heights vary from a minimum of about 105 meters below sea level to a maximum of about 85 meters above sea level.

Begin by loading the topo and geoid regular data grids:
```
load topo
load geoid
```
Create a map axes using a Gall stereographic cylindrical projection (a perspective projection):
```
axesm gstereo
```
Use meshm to plot a colored display of the geoid's variations, but specify topo as the final argument, to give each geoid grid cell the height (z-value) of the corresponding topo grid cell:
```
meshm(geoid,geoidrefvec,size(geoid),topo)
```
Low geoid heights are shown as blue, high ones as red.
For reference, plot the world coastlines in black, raise their elevation to 1000 meters (high enough to clear the surface in their vicinity), and expand the map to fill the frame:
```
load coast
plotm(lat,long,'k')
zdatam(handlem('allline'),1000)
tightmap
```
At this point the map looks like this.
Due to the vertical view and lack of lighting, the topographic relief is not visible, but it is part of the figure's surface data. Bring it out by exaggerating relief greatly, and then setting a view from the south-southeast:
```
daspectm('m',200); tightmap
view(20,35)
```
Remove the bounding box, shine a light on the surface (using the default position, offset to the right of the viewpoint), and render again with Phong shading:
```
set(gca,'Box','off')
camlight;
lighting phong
```
Finally, set the perspective to converge slightly (the default perspective is orthographic):
```
set(gca,'projection','perspective')
```
The final map is shown below. From it, you can see that the geoid mirrors the topography of the major mountain chains such as the Andes, the Himalayas, and the Mid-Atlantic Ridge. You can also see that large areas of high or low geoid heights are not simply a result of topography.

Draping Data over Terrain with Different Gridding

If you want to combine elevation and attribute (color) data grids that cover the same region but are gridded differently, you must resample one matrix to be consistent with the other. It helps if at least one of the grids is a geolocated data grid, because their explicit horizontal coordinates allow them to be resampled using the ltln2val function. To combine dissimilar grids, you can do one of the following:

Construct a geolocated grid version of the regular data grid values.
Construct a regular grid version of the geolocated data grid values.

The following two examples illustrate these closely related approaches.

Draping via Converting a Regular Grid to a Geolocated Data Grid

This example drapes slope data from a regular data grid on top of elevation data from a geolocated data grid. Although the two data sets actually have the same origin (the geolocated grid derives from the topo data set), this approach works with any dissimilar grids. The example uses the geolocated data grid as the source for surface elevations and transforms the regular data grid into slope values, which are then sampled to conform to the geolocated data grid (creating a set of slope values for the diamond-shaped grid) and color-coded for surface display.

Note When you use ltln2val to resample a regular data grid over an irregular area, make sure that the regular data grid completely covers the area of the geolocated data grid.

Begin by loading the geolocated data grids from mapmtx, which contains two regions. You will only use the diamond-shaped portion of mapmtx (lt1, lg1, and map1) centered on the Middle East, not the lt2, lg2, and map1 data:
```
load mapmtx lt1
load mapmtx lg1
load mapmtx map1
```
Load the topo global regular data grid:
```
load topo
```
Compute surface aspect, slope, and gradients for topo. You use only the slopes in subsequent steps:
```
[aspect,slope,gradN,gradE] = gradientm(topo,topolegend);
```
Use ltln2val to interpolate slope values to the geolocated grid specified by lt1, lg1:
```
slope1 = ltln2val(slope,topolegend,lt1,lg1);
```
The output is a 50-by-50 grid of elevations matching the coverage of the map1 variable.
Set up a figure with a Miller projection and use surfm to display the slope data. Specify the z-values for the surface explicitly as the map1 data, which is terrain elevation:
```
figure; axesm miller
surfm(lt1,lg1,slope1,map1)
```
The map mainly depicts steep cliffs, which represent mountains (the Himalayas in the northeast), and continental shelves and trenches.
The coloration depicts steepness of slope. Change the colormap to make the steepest slopes magenta, the gentler slopes dark blue, and the flat areas light blue:
```
colormap cool;
```
Use view to get a southeast perspective of the surface from a low viewpoint:
```
view(20,30); daspectm('m',200)
```
In 3-D, you immediately see the topography as well as the slope.
The default rendering uses faceted shading (no smooth interpolation). Render the surface again, this time making it shiny with Phong shading and lighting from the east (the default of camlight lights surfaces from over the viewer's right shoulder):
```
material shiny;camlight;lighting phong
```
Finally, remove white space and re-render the figure in perspective mode:
```
axis tight; set(gca,'Projection','Perspective')
```
Here is the mapped result.

Draping a Geolocated Grid on Regular Data Grid via Texture Mapping

The second way to combine a regular and a geolocated data grid is to construct a regular data grid of your geolocated data grid's z-data. This approach has the advantage that more computational functions are available for regular data grids than for geolocated ones. Another aspect is that the color and elevation grids do not have to be the same size. If the resolutions of the two are different, you can create the surface as a three-dimensional elevation map and later apply the colors as a texture map. You do this by setting the surface Cdata property to contain the color matrix, and setting the surface face color to 'TextureMap'.

In the following steps, you create a new regular data grid that covers the region of the geolocated data grid, then embed the color data values into the new matrix. The new matrix might need to have somewhat lower resolution than the original, to ensure that every cell in the new map receives a value.

Load the topo and terrain data from mapmtx:

load topo; 
load mapmtx lt1
load mapmtx lg1
load mapmtx map1

Identify the geographic limits of one of the mapmtx geolocated grids:

latlim = [min(lt1(:)) max(lt1(:))];
lonlim = [min(lg1(:)) max(lg1(:))];

Trim the topo data to the rectangular region enclosing the smaller grid:
```
[topo1,topo1ref] = maptrims(topo,topolegend,latlim,lonlim);
```
Create a regular grid filled with NaNs to receive texture data:
```
[curve1,curve1ref] = nanm(latlim,lonlim,.5);
```
The last parameter establishes the grid at 1/.5 cells per degree.
Use imbedm to embed values from map1 into the curve1 grid; the values are the discrete Laplacian transform (the difference between each element of the map1 grid and the average of its four orthogonal neighbors):
```
curve1 = imbedm(lt1,lg1,del2(map1),curve1,curve1ref);
```
Set up a map axes with the Miller projection and use meshm to draw the topo1 extract of the topo DEM:
```
figure; axesm miller
h = meshm(topo1,topo1ref,size(topo1),topo1);
```
Render the figure as a 3-D view from a 20º azimuth and 30º altitude, and exaggerate the vertical dimension by a factor of 200:
```
view(20,30); daspectm('m',200)
```
Light the view and render with Phong shading in perspective:
```
material shiny; camlight; lighting phong
axis tight; set(gca,'Projection','Perspective')
```
So far, both the surface relief and coloring represent topographic elevation.
Apply the curve1 matrix as a texture map directly to the figure using the set function:
```
set(h,'Cdata',curve1,'FaceColor','TextureMap')
```
The area originally covered by the [lt1, lg1, map1] geolocated data grid, and recoded via the Laplacian transform as curve1, now controls color symbolism, with the NaN-coded outside cells rendered in black.

Shading and Lighting Terrain Maps

Working with the Globe Display

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