Compositing and Animating Web Map Service (WMS) Meteorological Layers

Understanding Basic WMS Terminology

If you are new to WMS, several key concepts are important to understand and are listed here.

Internet Access

Since WMS servers are located on the Internet, this example can be set to access the Internet to dynamically render and retrieve maps from WMS servers or it can be set to use data previously retrieved from the Internet using the WMS capabilities but now stored in local files. You can use a variable, useInternet, to determine whether to read data from locally stored files, or retrieve the data from the Internet.

If the useInternet flag is set to true, then an Internet connection must be established to run the example. Note that the WMS servers may be unavailable, and several minutes may elapse before the maps are returned. One of the challenges of working with WMS servers is that sometimes you will encounter server errors. A function, such as wmsread, may time out if a server is unavailable. Often, this is a temporary problem and you will be able to connect to the server if you try again later.

This example writes data to files if useInternet is true or reads data from files if useInternet is false. You can store the data locally the first time you run the example and then set the useInternet flag to false.

useInternet = false;

Step 1: Find Katrina Layers From Local Database

One of the more challenging aspects of using WMS is finding a WMS server and then finding the layer that is of interest to you. The process of finding a server that contains the data you need and constructing a specific and often complicated URL with all the relevant details can be very daunting.

The Mapping Toolbox™ simplifies the process of locating WMS servers and layers by providing a local, installed, and pre-qualified WMS database, that is searchable, using the function wmsfind. You can search the database for layers and servers that are of interest to you. Here is how you find layers containing the term katrina in either the LayerName or LayerTitle field of the database:

katrina = wmsfind('katrina');
whos katrina

  Name          Size            Bytes  Class       Attributes

  katrina      34x1             16754  WMSLayer              

The search for the term 'katrina' returned a WMSLayer array containing multiple layers. To inspect information about an individual layer, simply display it like this:

katrina(1)

ans = 
  WMSLayer

  Properties:
           Index: 1
     ServerTitle: 'NASA SVS Image Server'
       ServerURL: 'https://svs.gsfc.nasa.gov/cgi-bin/wms?'
      LayerTitle: 'GOES-12 Imagery of Hurricane Katrina: Longwave Infrared Close-up (1024x1024 Animation)'
       LayerName: '3216_22510'
          Latlim: [15.0000 45.0000]
          Lonlim: [-100.0000 -70.0000]

  Methods

If you type, katrina, in the command window, the entire contents of the array are displayed, with each element's index number included in the output. This display makes it easy for you to examine the entire array quickly, searching for a layer of interest. You can display only the LayerTitle property for each element by executing the command:

disp(katrina,'Properties','layertitle','Index','off','Label','off');

As you've discovered, a search for the generic word 'katrina' returned results of many layers and you need to select only one layer. In general, a search may even return thousands of layers, which may be too large to review individually. Rather than searching the database again, you may refine your search by using the refine method of the WMSLayer class. Using the refine method is more efficient and returns results faster than wmsfind since the search has already been narrowed to a smaller set. Supplying the query string, 'goes-12*katrina*visible*close*up*animation', to the refine method returns a WMSLayer array whose elements contain a match of the query string in either the LayerTitle or LayerName properties. The * character indicates a wild-card search. If multiple entries are returned, select only the first one from the svs.gsfc.nasa.gov server.

katrina = refine(katrina,'goes-12*katrina*visible*close*up*animation');
katrina = refine(katrina,'svs.gsfc.nasa.gov','Searchfield','serverurl');
katrina = katrina(1);
whos katrina

  Name         Size            Bytes  Class       Attributes

  katrina      1x1               466  WMSLayer              

Step 2: Synchronize WMSLayer Object with Server

The database only stores a subset of the layer information. For example, information from the layer's abstract, details about the layer's attributes and style information, and the coordinate reference system of the layer are not returned by wmsfind. To return all the information, you need to use the wmsupdate function. wmsupdate synchronizes the layer from the database with the server, filling in the missing properties of the layer.

Synchronize the first katrina layer with the server in order to obtain the abstract information. Since this action requires access to the Internet, call wmsupdate only if the useInternet flag is true.

if useInternet
    katrina = wmsupdate(katrina);
    save('katrina.mat','katrina')
else
    load('katrina.mat')
end

Display the abstract information of the layer. Use isspace to help determine where to line wrap the text.

abstract = katrina.Abstract;
endOfLine = find(isstrprop(abstract,'cntrl'),1);
abstract = abstract(1:endOfLine);
numSpaces = 60;
while(~isempty(abstract))
    k = find(isspace(abstract));
    n = find(k > numSpaces,1);
    if ~isempty(n)
        fprintf('%s\n',abstract(1:k(n)))
        abstract(1:k(n)) = [];
    else
        fprintf('%s\n',abstract)
        abstract = '';
    end
end

The GOES-12 satellite sits at 75 degrees west longitude at an 
altitude of 36,000 kilometers over the equator, in geosynchronous 
orbit.  At this position its Imager instrument takes pictures 
of cloud patterns in several wavelengths for all of North and 
South America, a primary measurement used in weather forecasting. 
 The Imager takes a pattern of pictures of parts of the Earth 
in several wavelengths all day, measurements that are vital in 
weather forecasting.  This animation shows a daily sequence of 
GOES-12 images in the visible wavelengths, from 0.52 to 0.72 
microns, during the period that Hurricane Katrina passed through 
the Gulf of Mexico.  At one kilometer resolution, the visible 
band measurement is the highest resolution data from the Imager, 
which accounts for the very high level of detail in these images. 
 For this animation, the cloud data was extracted from GOES image 
and laid over a background color image of the southeast United 

States.

Note that this abstract information, including any typographical issues and incomplete fragments, was obtained directly from the server.

Step 3: Explore Katrina Layer Details

You can find out more information about the katrina layer by exploring the Details property of the katrina layer. The Details.Attributes field informs you that the layer has fixed width and fixed height attributes, thus the size of the requested map cannot be modified.

katrina.Details.Attributes

ans = struct with fields:
      Queryable: 0
       Cascaded: 0
         Opaque: 1
      NoSubsets: 1
     FixedWidth: 1024
    FixedHeight: 1024

The Details.Dimension field informs you that the layer has a time dimension

katrina.Details.Dimension

ans = struct with fields:
              Name: 'time'
             Units: 'ISO8601'
        UnitSymbol: ''
           Default: '2005-08-30T17:45Z'
    MultipleValues: 0
      NearestValue: 0
           Current: 0
            Extent: '2005-08-23T17:45Z/2005-08-30T17:45Z/P1D'

with an extent from 2005-08-23T17:45Z to 2005-08-30T17:45Z with a period of P1D (one day), as shown in the Details.Dimension.Extent field.

katrina.Details.Dimension.Extent

ans = 
'2005-08-23T17:45Z/2005-08-30T17:45Z/P1D'

Step 4: Retrieve Katrina Map from Server

Now that you have found a layer of interest, you can retrieve the raster map using the function wmsread and display the map using the function geoshow. Since Time is not specified when reading the layer, the default time, 2005-08-30T17:45Z, is retrieved as specified by the Details.Dimension.Default field. If the useInternet flag is set to true, then cache the image and reference object in a GeoTIFF file.

if useInternet
    [katrinaMap,R] = wmsread(katrina);
    geotiffwrite('katrina.tif',katrinaMap,R)
else
    [katrinaMap,R] = readgeoraster('katrina.tif');
end

Display the katrinaMap and overlay the data from the usastatehi.shp file.

states = readgeotable('usastatehi.shp');
figure
usamap(katrina.Latlim, katrina.Lonlim)
geoshow(katrinaMap,R)
geoshow(states,'FaceColor','none')
title({katrina.LayerTitle, katrina.Details.Dimension.Default}, ...
    'Interpreter','none')

Step 5: Find NEXRAD Radar Layer

NEXRAD radar images for the United States are stored on the Iowa State University's IEM Web map server. The server conveniently stores NEXRAD images in five minute increments from 1995-01-01 to the present time. You can find the layer by first searching for the term IEM WMS Service in the ServerTitle field of the WMS database, then refining the search by requesting the layer of interest, nexrad-n0r-wmst.

iemLayers = wmsfind('IEM WMS Service','SearchField','servertitle');
nexrad = refine(iemLayers,'nexrad-n0r-wmst');

Synchronize the layer with the server.

if useInternet
    nexrad = wmsupdate(nexrad);
    save('nexrad.mat','nexrad')
else
    load('nexrad.mat');
end

Step 6: Obtain Extent Parameters

To composite the nexrad layer with the katrina layer, you need to obtain the nexrad layer at coincidental time periods, and concurrent geographic and image extents. The Details.Dimension field informs you that the layer has a time dimension,

nexrad.Details.Dimension

ans = struct with fields:
              Name: 'time'
             Units: 'ISO8601'
        UnitSymbol: ''
           Default: '2006-06-23T03:10:00Z'
    MultipleValues: 0
      NearestValue: 0
           Current: 0
            Extent: '1995-01-01/2022-12-31/PT5M'

and the Details.Dimension.Default field informs you that the layer's time extent includes seconds.

nexrad.Details.Dimension.Default

ans = 
'2006-06-23T03:10:00Z'

Obtain a time value coincidental with the katrina layer, and add seconds to the time specification.

nexradTime = [katrina.Details.Dimension.Default(1:end-1) ':00Z'];

Assign latlim and lonlim variables to specify the limits for the nexrad layer. Set the values to the limits of the katrina layer so that the geographic areas match. Note that the nexrad layer's southern latitude limit does not extend as far south as the katrina layer's southern latitude limit. The values that lie outside the geographic bounding quadrangle of the nexrad layer are set to the background color.

fprintf('%s%d\n','Southern latitude limit of NEXRAD  layer: ',nexrad.Latlim(1))

Southern latitude limit of NEXRAD  layer: 24

fprintf('%s%d\n','Southern latitude limit of Katrina layer: ',katrina.Latlim(1))

Southern latitude limit of Katrina layer: 10

latlim = katrina.Latlim;
lonlim = katrina.Lonlim;

Assign imageHeight and imageWidth variables.

imageHeight = katrina.Details.Attributes.FixedHeight;
imageWidth  = katrina.Details.Attributes.FixedWidth;

Step 7: Retrieve NEXRAD Radar Map from Server

You can retrieve the nexradMap from the server, specified at the same time as the katrinaMap and for the same geographic and image extents, by supplying parameter/value pairs to the wmsread function. To accurately retrieve the radar signal from the map, set the ImageFormat parameter to the image/png format. In order to easily retrieve the signal from the background, set the background color to black ([0 0 0]).

Retrieve the nexradMap.

black = [0 0 0];
if useInternet
    [nexradMap,R] = wmsread(nexrad, ...
        'Latlim',latlim,'Lonlim',lonlim,'Time',nexradTime, ...
        'BackgroundColor',black,'ImageFormat','image/png', ...
        'ImageHeight',imageHeight,'ImageWidth',imageWidth);
    geotiffwrite('nexrad.tif',nexradMap,R)
else
    [nexradMap,R] = readgeoraster('nexrad.tif');
end

Display the nexradMap.

figure
usamap(latlim,lonlim)
geoshow(nexradMap,R)
geoshow(states,'FaceColor','none','EdgeColor',[0.9 0.9 0.9])
title({nexrad.LayerTitle, nexradTime},'Interpreter','none');

Step 8: Composite NEXRAD Radar Map with Katrina Map

To composite the nexradMap with a copy of the katrinaMap, you need to identify the non-background pixels in the nexradMap. The nexradMap data is returned as an image with class double, because of how this web map server handles PNG format, so you need convert it to uint8 before merging.

Identify the pixels of the nexradMap image that do not contain the background color.

threshold = 0;
index = any(nexradMap > threshold, 3);
index = repmat(index,[1 1 3]);

Composite the nexradMap with the katrinaMap.

combination = katrinaMap;
combination(index) = uint8(nexradMap(index)*255);

Display the composited map.

figure
usamap(latlim,lonlim)
geoshow(combination,R)
geoshow(states,'FaceColor','none')
title({'GOES 12 Imagery of Hurricane Katrina', ...
    'Composited with NEXRAD Radar',nexradTime})

Step 9: Initialize Variables to Animate the Katrina and NEXRAD Maps

The next step is to initialize variables in order to animate the composited katrina and nexrad maps.

Create variables that contain the time extent of the katrina layer.

extent = katrina.Details.Dimension.Extent;
slash = '/';
slashIndex = strfind(extent,slash);
startTime = extent(1:slashIndex(1)-1);
endTime = extent(slashIndex(1)+1:slashIndex(2)-1);

Calculate numeric values for the start and end days. Note that the time extent is in yyyy-mm-dd format.

hyphen = '-';
hyphenIndex = strfind(startTime,hyphen);
dayIndex = [hyphenIndex(2) + 1, hyphenIndex(2) + 2];
startDay = str2double(startTime(dayIndex));
endDay = str2double(endTime(dayIndex));

Assign the initial katrinaTime.

katrinaTime = startTime;

Since multiple requests to a server are required for animation, it is more efficient to use the WebMapServer and WMSMapRequest classes.

Construct a WebMapServer object for each layer's server.

nasaServer = WebMapServer(katrina.ServerURL);
iemServer  = WebMapServer(nexrad.ServerURL);

Create WMSMapRequest objects.

katrinaRequest = WMSMapRequest(katrina, nasaServer);
nexradRequest  = WMSMapRequest(nexrad, iemServer);

Assign properties.

nexradRequest.Latlim = latlim;
nexradRequest.Lonlim = lonlim;
nexradRequest.BackgroundColor = black;
nexradRequest.ImageFormat = 'image/png';
nexradRequest.ImageHeight = imageHeight;
nexradRequest.ImageWidth  = imageWidth;

Step 10: Create Animation Files

An animation can be viewed in the browser when the browser opens an animated GIF file or an AVI video file. To create the animation frames of the WMS basemap and vector overlays, create a loop through each day, from startDay to endDay, and obtain the katrinaMap and the nexradMap for that day. Composite the maps into a single image, display the image, retrieve the frame, and store the results into a frame of an AVI file and a frame of an animated GIF file.

To share with others or to post to web video services, create an AVI video file containing all the frames using the VideoWriter class.

videoFilename = fullfile(pwd,'wmsanimated.avi');
if exist(videoFilename,'file')
    delete(videoFilename)
end
writer = VideoWriter(videoFilename);
writer.FrameRate = 1;
writer.Quality = 100;
open(writer)

The animation is viewed in a single map display. Outside the animation loop, create a map display. Initialize hmap, used in the loop as the return handle from the function geoshow, so it can be deleted on the first pass through the loop. Loop through each day, retrieve and display the WMS map, and save the frame.

fig = figure;
usamap(latlim,lonlim)
hstates = geoshow(states,'FaceColor','none');
hmap = [];

for k = startDay:endDay
    
    % Update the time values and assign the Time property for each server.
    currentDay = num2str(k);
    katrinaTime(dayIndex) = currentDay;
    nexradTime = [katrinaTime(1:end-1) ':00Z'];
    katrinaRequest.Time = katrinaTime;
    nexradRequest.Time = nexradTime;
    
    % Retrieve the WMS map of Katrina from the server (or file)
    % for this time period.
    nmKatrina = ['katrina_' num2str(currentDay) '.tif'];
    if useInternet
        katrinaMap = getMap(nasaServer, katrinaRequest.RequestURL);
        geotiffwrite(nmKatrina, katrinaMap, katrinaRequest.RasterRef)
    else
        katrinaMap = readgeoraster(nmKatrina);
    end
    
    % Retrieve the WMS map of the NEXRAD imagery from the server (or file)
    % for this time period.
    nmNEXRAD = ['nexrad_' num2str(currentDay) '.tif'];
    if useInternet
        nexradMap = getMap(iemServer, nexradRequest.RequestURL);
        geotiffwrite(nmNEXRAD, nexradMap, nexradRequest.RasterRef)
    else
        nexradMap = readgeoraster(nmNEXRAD);
    end
    
    % Identify the pixels of the nexradMap image that do not contain the
    % background color.
    index = any(nexradMap > threshold, 3);
    index = repmat(index,[1 1 3]);
    
    % Composite nexradMap with katrinaMap.
    combination = katrinaMap;
    combination(index) = uint8(nexradMap(index)*255);
    
    % Delete the old map and display the new composited map.
    delete(hmap)
    hmap = geoshow(combination, katrinaRequest.RasterRef);
    uistack(hstates,'top')
    title({'GOES 12 Imagery of Hurricane Katrina', ...
        'Composited with NEXRAD Radar',nexradTime})
    drawnow
    
    % Save the current frame as an RGB image.
    currentFrame = getframe(fig);
    RGB = currentFrame.cdata;
    
    % Create an indexed image for each RGB frame in order to display an
    % animated GIF. 
    if k == startDay
        % The first time through the loop, convert the RGB image to
        % an indexed image and save the colormap into the
        % variable, cmap. Use cmap to convert later frames.
        [frame,cmap] = rgb2ind(RGB,256,'nodither');
        
        % Use the size of the first frame and the total
        % number of frames to initialize animated with
        % a size large enough to contain all the frames.
        frameSize = size(frame);
        numFrames = endDay - startDay + 1;
        animated = zeros([frameSize 1 numFrames],'like',frame);
    else
        % Use the colormap from the first frame conversion and
        % convert this frame to an indexed image.
        frame = rgb2ind(RGB,cmap,'nodither');
    end
   
    % Store the frame into the animated array for the GIF file.
    frameCount = k - startDay + 1;
    animated(:,:,1,frameCount) = frame;
    
    % Write the RGB frame to the AVI file.
    writeVideo(writer,RGB);
end

Close the Figure window and the AVI file.

close(fig)
close(writer)

Write the animated GIF file.

filename = 'wmsanimated.gif';
delayTime = 2.0;
loopCount = inf;
imwrite(animated,cmap,filename, ...
    'DelayTime',delayTime,'LoopCount',loopCount);

Step 11: View Animated GIF File

An animation can be viewed in the browser when the browser opens an animated GIF file.

Credits

Katrina Layer

The Katrina layer used in the example is from the NASA Goddard Space Flight Center's SVS Image Server and is maintained by the Scientific Visualization Studio.

For more information about this server, use the command wmsinfo('http://svs.gsfc.nasa.gov/cgi-bin/wms?').

NEXRAD Layer

The NEXRAD layer used in the example is from the Iowa State University's IEM WMS server and is a generated CONUS composite of National Weather Service (NWS) WSR-88D level III base reflectivity.

For more information about this server, use the command wmsinfo('http://mesonet.agron.iastate.edu/cgi-bin/wms/nexrad/n0r-t.cgi?').