Retrieve Data from Web Map Server

Typically, a WMS server returns a pictorial representation of a layer (or layers) back to a requester 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 requester. The wmsread function, or the getMap method of the WebMapServer object, 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.

Note

To request actual data, most often you need to create either a Web Coverage Service (WCS) request, for raster data, or a Web Feature Service (WFS) request, for vector data. The Mapping Toolbox™ does not support WCS and WFS requests.

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('https://data.worldwind.arc.nasa.gov/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('data.worldwind*elev', 'SearchField', 'serverurl');
layers = wmsupdate(layers);

  2. Select the 'EarthAsterElevations30m' layer containing SRTM30 data merged with global ASTER data.

    aster = layers.refine('EarthAsterElevations30m', ...
                      'SearchField', 'layername');

  3. Define the region surrounding Afghanistan.

    latlim = [25 40];
lonlim = [55 80];

  4. 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');

  5. Display the elevation data as a texture map. 

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

  6. 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
worldmap('Afghanistan')
demcmap(double(ZA))
geoshow(ZB, RA, 'DisplayType', 'texturemap')
title({'Afghanistan and Country Boundaries', aster.LayerTitle});

Display 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('data.worldwind*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 raster reference object 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
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
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.

    disp(newline + "SRTM30 Plus - " + srtmplus.LayerName ...
+ newline + "Minimum value: " + min(Z1(:))...
+ newline + "Maximum value: " + max(Z1(:)))
disp(newline + "SRTM30 Plus Merged - " + mergedSrtm.LayerName ...
+ newline + "Minimum value: " + min(Z2(:))...
+ newline + "Maximum value: " + 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('/USGSImageryTopo/','SearchField','serverurl'); 

layers = wmsfind('data.worldwind', '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
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);

See Also

| |

Related Topics

Mapping Toolbox Documentation

Support