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.

A WMS server renders one or more layers and stores the results in a file, which is streamed to the requester. The wmsread and getMap functions make the request on your behalf, capture the stream in a temporary file, and import the file content into a variable into a workspace variable.

Different file formats might be available, such as JPEG, PNG, GIF, and GeoTIFF. Almost all WMS servers support use of the JPEG format, and many support more than one standard graphics format.

The format might affect the quality of your results. For example, the PNG format avoids the tile-related artifacts that are common with JPEG. If a server supports multiple formats, you can specify the format by using the wmsread function and the ImageFormat name-value argument; such as 'ImageFormat','image/png'.

The format also determines whether you 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 GeoTIFF and BIL formats.

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.

To specify the GeoTIFF or BIL formats, use the name-value arguments 'ImageFormat','image/tiff' or 'ImageFormat','image/bil', respectively. Although GeoTIFF and BIL files might contain multiple co-registered bands (channels), the GeoTIFF and BIL files returned by WMS servers include only a single band and the output of wmsread is a 2-D array.

For example, you can obtain signed, quantitative elevation data, rather than an RGB image, from the NASA WorldWind WMS server. See the output from the command:

wmsinfo('https://worldwind26.arc.nasa.gov/wms/elev?')

After retrieving data in a quantitative format, 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.

Get information about the NASA WorldWind server. Select a layer containing SRTM data.

info = wmsinfo('https://worldwind26.arc.nasa.gov/wms/elev?');
layers = info.Layer;
srtm = refine(layers,'SRTM','SearchField','layername');

Specify the latitude and longitude limits for a region surrounding Afghanistan.
```
latlim = [25 40];
lonlim = [55 80];
```

Read the data using a 1-minute sampling interval. Specify the 'ImageFormat' as 'image/tiff'.

cellSize = dms2degrees([0,1,0]);
[ZA,RA] = wmsread(srtm,'Latlim',latlim,'Lonlim',lonlim, ...
   'CellSize',cellSize,'ImageFormat','image/tiff');

Note that ZA is an array of elevation data rather than an RGB image. Find the minimum and maximum elevations.
```
elevationLimits = [min(min(ZA)) max(max(ZA))]
```
```
elevationLimits = 1×2 int16 row vector

   -1415    7715
```

Display the elevation data as a texture map.

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

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);
ZB = ZA;
ZB(B(:,:,1) < 250) = min(ZA(:));
figure
worldmap('Afghanistan')
demcmap(double(ZA))
geoshow(ZB,RA,'DisplayType','texturemap')
title({'Afghanistan and Country Boundaries',srtm.LayerTitle}, ...
    'Interpreter','none');

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.

For this example, 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.

Get information about the NASA WorldWind server. Find the 'srtm30' layer, which is the SRTM30 Plus data set.

info = wmsinfo('https://worldwind26.arc.nasa.gov/wms/elev?');
layers = info.Layer;
srtm = refine(layers,'srtm30','SearchField','layername');

Specify the latitude and longitude limits for a region surrounding the Gulf of Maine. Set the sampling interval to 30 arc-seconds and the image format to 'image/tiff'. Request the map from the server.
```
latlim = [40 46];
lonlim = [-71 -65];
samplesPerInterval = dms2degrees([0 0 30]);
imageFormat = 'image/tiff';
[Z,R] = wmsread(srtm,'Latlim',latlim, ...
   'Lonlim',lonlim,'ImageFormat',imageFormat, ...
   'CellSize',samplesPerInterval);
```
Note that Z is an array of elevation data rather than an RGB image. Find minimum and maximum values of Z.
```
elevationLimits = [min(min(Z)) max(max(Z))]
```
```
elevationLimits = 1×2 int16 row vector

   -4540    1469
```
Create a map. The raster reference object R1 ties the intrinsic coordinates of the raster map to the EPSG:4326 geographic coordinate system. Apply a colormap appropriate for elevation data. Display and contour the map at sea level.
```
figure
worldmap(Z,R)
geoshow(Z,R,'DisplayType','texturemap')
demcmap(double(Z))
contourm(double(Z),R,[0 0],'Color','k')
colorbar
title ({'Gulf of Maine',srtm.LayerTitle}, ...
   'Interpreter','none')
```

Drape WMS Imagery onto Elevation Data

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

Find layers containing the orthoimage and elevation data.

ortho = wmsfind('/USGSImageryTopo/','SearchField','serverurl'); 
ortho = wmsupdate(ortho);

info = wmsinfo('https://worldwind26.arc.nasa.gov/wms/elev?');
layers = info.Layer;
us_ned = refine(layers,'usgs-ned');

Specify the latitude and longitude limits for an area around the Grand Canyon. Specify an image size.
```
latlim = [36 36.23];
lonlim = [-113.36 -113.13];
imageHeight = 575;
imageWidth = 575;
```

Read the layers. Specify the 'ImageFormat' of the elevation layer as 'image/tiff'.

A = wmsread(ortho,'Latlim',latlim,'Lonlim',lonlim, ...
   'ImageHeight',imageHeight,'ImageWidth',imageWidth);

[Z,R] = wmsread(us_ned,'ImageFormat','image/tiff', ...
    'Latlim',latlim,'Lonlim',lonlim, ...
    'ImageHeight',imageHeight,'ImageWidth',imageWidth);

Note that Z is an array of elevation data rather than an RGB image.

Drape the orthoimage 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

Adjust the camera view and lighting.

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

Documentation