Retrieving Elevation Data

Typically, a WMS server returns a pictorial representation of a layer (or layers) back to a requestor rather than the actual data. However, in some rare cases, you can request actual data from specific WMS servers, by using a certain option with wmsread.

A WMS server renders one or more layers and stores the results in a file, which is streamed to the requestor. The wmsread function, or the getMap method of the WebMapServer class, makes the request on your behalf, captures the stream in a temporary file, and imports the file content into a variable in your MATLAB® workspace. The format of the file may be a standard graphics format, such as JPEG, PNG, or GIF, or it may be the band-interleaved by line (BIL) format, which is popular in remote sensing. Almost all WMS servers support use of the JPEG format, and many support more than one standard graphics format. Only a very few WMS servers support the BIL format, although it is very useful.

The choice of format can affect the quality of your results. For example, PNG format avoids the tile-related artifacts that are common with JPEG. Format choice also determines whether you will get a pictorial representation, which is the case with any of the standard graphics formats, or an absolutely quantitative data grid (possibly including negative as well as positive values). Quantitative data sets are provided via the BIL format.

With a server that supports multiple formats, you can control which format is used by specifying an ImageFormat name-value pair when calling wmsread. For example, if the server supports the PNG format you would choose PNG by specifying 'ImageFormat','image/png', thus avoiding the possibility of JPEG artifacts.

With a server that supports it, you can obtain an absolutely quantitative data grid by specifying BIL format when calling wmsread. To do this, use the name-value pair 'ImageFormat','image/bil'. Although a BIL file typically contains multiple, co-registered bands (channels), the BIL files returned by a WMS server include only a single band. In this case, the output of wmsread enters the MATLAB workspace as a 2-D array.

For example, you can obtain signed, quantitative elevation data, rather than an RGB image, from the NASA WorldWind WMS server (the only server in the Mapping Toolbox™ WMS database known to support the 'image/bil' option). See the output from the command:

wmsinfo('http://www.nasa.network.com/elev?')

In fact, because the NASA WorldWind WMS server returns rendered layers only in the 'image/bil' format, you do not need to provide an 'ImageFormat' name-value pair when using this server. However, it is good practice to specify the image format, in case the server is ever updated to provide rendered layers in other image formats.

After retrieving data using the 'ImageFormat','image/bil' option, display it as a surface or a texture-mapped surface, rather than as an image, as shown in the examples below.

Merge Elevation Data with Rasterized Vector Data

The NASA WorldWind WMS server contains a wide selection of layers containing elevation data. Follow this example to merge elevation data with a raster map containing national boundaries.

  1. Find the layers from the NASA WorldWind server.

    layers = wmsfind('nasa.network*elev', 'SearchField', 'serverurl');
    layers = wmsupdate(layers);
    
  2. Display the name and title of each layer.

    disp(layers,'Properties',{'LayerTitle','LayerName'})
    
     11x1 WMSLayer
    
      Properties:
               Index: 1
          LayerTitle: 'SRTM30 with Bathymetry (900m) merged with global ASTER (30m)'
           LayerName: 'EarthAsterElevations30m'
    
               Index: 2
          LayerTitle: 'USGS NED 30m'
           LayerName: 'NED'
    
               Index: 3
          LayerTitle: 'ScankortElevationsDenmarkDSM'
           LayerName: 'ScankortElevationsDenmarkDSM'
                    .
                    .
                    .
               Index: 10
          LayerTitle: 'SRTM30 Plus'
           LayerName: 'srtm30'
    
               Index: 11
          LayerTitle: 'USGS NED 10m'
           LayerName: 'usgs_ned_10m'
    
  3. Select the 'EarthAsterElevations30m' layer containing SRTM30 data merged with global ASTER data.

    aster = layers.refine('earthaster', 'SearchField', 'layername');
    
  4. Define the region surrounding Afghanistan.

    latlim = [25 40];
    lonlim = [55 80];
    
  5. Obtain the data at a 1-minute sampling interval.

    cellSize = dms2degrees([0,1,0]);
    [ZA, RA] = wmsread(aster, 'Latlim', latlim, 'Lonlim', lonlim, ...
       'CellSize', cellSize, 'ImageFormat', 'image/bil');
    
  6. Display the elevation data as a texture map.

    figure('Renderer','zbuffer')
    worldmap('Afghanistan')
    geoshow(ZA, RA, 'DisplayType', 'texturemap')
    demcmap(double(ZA))
    title({'Afghanistan and Surrounding Region', aster.LayerTitle});
    

  7. Embed national boundaries from the VMAP0 WMS server into the elevation map.

    vmap0 = wmsfind('vmap0.tiles', 'SearchField', 'serverurl');
    boundaries = refine(vmap0, 'country_02');
    B = wmsread(boundaries, 'Latlim', latlim, ...
        'Lonlim', lonlim, 'CellSize', cellSize, 'ImageFormat','image/png');
    ZB = ZA;
    ZB(B(:,:,1) < 250) = min(ZA(:));
    figure('Renderer','zbuffer')
    worldmap('Afghanistan')
    demcmap(double(ZA))
    geoshow(ZB, RA, 'DisplayType', 'texturemap')
    title({'Afghanistan and Country Boundaries', aster.LayerTitle});
    

Display a Merged Elevation and Bathymetry Layer (SRTM30 Plus)

The Shuttle Radar Topography Mission (SRTM) is a project led by the U.S. National Geospatial-Intelligence Agency (NGA) and NASA. SRTM has created a high-resolution, digital, topographic database of Earth. The SRTM30 Plus data set combines GTOPO30, SRTM-derived land elevation and Sandwell bathymetry data from the University of California at San Diego.

Follow this example to read and display the SRTM30 Plus layer for the Gulf of Maine at a 30 arc-second sampling interval using data from the WorldWind server.

  1. Find and update the 'srtm30' layer in the WMS Database. The 'srtm30' layer name from NASA WorldWind is the name for the SRTM30 Plus data set.

    wldwind = wmsfind('nasa.network*elev', 'SearchField', 'serverurl');
    wldwind = wmsupdate(wldwind);
    srtmplus = wldwind.refine('srtm30', 'SearchField', 'layername');
    
  2. Set the desired geographic limits.

    latlim = [40 46];
    lonlim = [-71 -65];
    
  3. Set the sampling interval to 30 arc-seconds.

    samplesPerInterval = dms2degrees([0 0 30]);
    
  4. Set the ImageFormat to image/bil.

    imageFormat = 'image/bil';
  5. Request the map from the NASA server.

    [Z1, R1] = wmsread(srtmplus, 'Latlim', latlim, ...
       'Lonlim', lonlim, 'ImageFormat', imageFormat, ...
       'CellSize', samplesPerInterval);
    
  6. Open a figure window and set up a map axes with geographic limits that match the desired limits. The referencing matrix R1 ties the intrinsic coordinates of the raster map to the EPSG:4326 geographic coordinate system. Create a colormap appropriate for elevation data. Then, display and contour the map at sea level (0 m).

    figure('Renderer','zbuffer')
    worldmap(Z1, R1)
    geoshow(Z1, R1, 'DisplayType', 'texturemap')
    demcmap(double(Z1))
    contourm(double(Z1), R1, [0 0], 'Color', 'black')
    colorbar
    title ({'Gulf of Maine', srtmplus.LayerTitle}, 'Interpreter','none')
    

  7. Compare the NASA WorldWind SRTM30 Plus layer with the SRTM30 with Bathymetry (900m) merged with SRTM3 V4.1 (90m) and USGS NED (30m) (mergedSrtm) layer.

    mergedSrtm = wldwind.refine('mergedSrtm');
  8. Request the map from the NASA WorldWind server.

    [Z2, R2] = wmsread(mergedSrtm, 'Latlim', latlim, 'Lonlim', lonlim, ...
        'CellSize', samplesPerInterval, 'ImageFormat', 'image/bil');
    
  9. Display the data.

    figure('Renderer','zbuffer')
    worldmap(Z2, R2)
    geoshow(Z2, R2, 'DisplayType', 'texturemap')
    demcmap(double(Z2))	
    contourm(double(Z2), R2, [0 0], 'Color', 'black')
    colorbar
    title ({'Gulf of Maine', mergedSrtm.LayerTitle})
    

  10. Compare the results.

    fprintf(...
       '\nSRTM30 Plus - %s\nMinimum value: %d\nMaximum value: %d\n', ...
       srtmplus.LayerName, min(Z1(:)), max(Z1(:)));
    fprintf(...
       '\nSRTM30 Plus Merged - %s\nMinimum value: %d\nMaximum value: %d\n', ...
       mergedSrtm.LayerName, min(Z2(:)), max(Z2(:)));
    
  11. The output appears as follows:

    SRTM30 Plus - srtm30 
    Minimum value: -4543
    Maximum value: 1463
    
    Merged SRTM30 Plus - mergedSrtm
    Minimum value: -4543 
    Maximum value: 1463
    

Drape WMS Imagery onto Elevation Data

This example shows how to drape WMS imagery onto elevation data from the USGS National Elevation Dataset (NED).

  1. Obtain the layers of interest.

    ortho = wmsfind('usgs*topographic*small*scale','SearchField','any');
    
    layers = wmsfind('nasa.network', 'SearchField', 'serverurl');
    us_ned = layers.refine('usgs ned 30');
    
  2. Assign geographic extent and image size.

    latlim = [36 36.23];
    lonlim = [-113.36 -113.13];
    imageHeight = 575;
    imageWidth = 575;
    
  3. Read the ortho layer.

    A = wmsread(ortho, 'Latlim', latlim, 'Lonlim', lonlim, ...
       'ImageHeight', imageHeight, 'ImageWidth', imageWidth);
    
  4. Read the USGS us_ned layer.

    [Z, R] = wmsread(us_ned, 'ImageFormat', 'image/bil', ...
        'Latlim', latlim, 'Lonlim', lonlim, ...
        'ImageHeight', imageHeight, 'ImageWidth', imageWidth);
    
  5. Drape the ortho image onto the elevation data.

    figure('Renderer','opengl')
    usamap(latlim, lonlim)
    framem off; mlabel off; plabel off; gridm off
    geoshow(double(Z), R, 'DisplayType', 'surface', 'CData', A);
    daspectm('m',1)
    title({'Grand Canyon', 'USGS NED and Ortho Image'}, ...
       'FontSize',8);
    axis vis3d
    

  6. Assign camera parameters.

    cameraPosition = [96431 4.2956e+06 -72027];
    cameraTarget = [-82.211 4.2805e+06 3054.6];
    cameraViewAngle = 8.1561;
    cameraUpVector = [3.8362e+06 5.9871e+05 5.05123e+006];
  7. Set camera and light parameters.

    set(gca,'CameraPosition', cameraPosition, ...
        'CameraTarget', cameraTarget, ...
        'CameraViewAngle', cameraViewAngle, ...
        'CameraUpVector', cameraUpVector);
    lightHandle = camlight;
    camLightPosition = [7169.3 1.4081e+06 -4.1188e+006];
    set(lightHandle, 'Position', camLightPosition);
    

Was this topic helpful?