Main Content

Animate Data Layers

Create Movie of Terra/MODIS Maps

You can create maps of the same geographic region at different times and view them as a movie. For a period of seven days, read and display a daily composite of visual images from NASA's Moderate Resolution Imaging Spectroradiometer (MODIS) scenes captured during the month of December 2010.

  1. Search the WMS Database for the MODIS layer.

    neo = wmsfind('neowms*nasa', 'SearchField', 'serverurl'); 
    modis = neo.refine('true*color*terra*modis'); 
    modis = wmsupdate(modis); 
  2. Construct a WebMapServer object.

    server = WebMapServer(modis.ServerURL);
  3. Construct a WMSMapRequest object.

    mapRequest = WMSMapRequest(modis, server);
  4. The Extent field provides the information about how to retrieve individual frames. You can request a single day since the extent is defined by day ('/P1D'). Note that for December 2010, the frames for December 8 and December 31 are not available.

  5. Create an array indicating the first seven days.

    days = 1:7;
  6. Set the value of startTime to December 01, 2010 and use a serial date number.

    time = '2010-12-01';
    startTime = datenum(time);
  7. Open a figure window with axes appropriate for the region specified by the modis layer.

    hFig = figure('Color', 'white');
    worldmap(mapRequest.Latlim, mapRequest.Lonlim);
  8. Save each frame into a video file.

    videoFilename = 'modis_dec.avi';
    writer = VideoWriter(videoFilename);
    writer.FrameRate = 1;
    writer.Quality = 100;;
  9. Retrieve a map of the modis layer for each requested day. Set the Time property to the day number. When obtaining the data from the server, use a try/catch statement to ignore either data not found on the server or any error issued by the server. Set startTime to one day less for correct indexing.

    startTime = startTime - 1;
    for k = days
            mapRequest.Time = startTime + k;
            timeStr = datestr(mapRequest.Time);
            dailyImage = server.getMap(mapRequest.RequestURL);
            geoshow(dailyImage, mapRequest.RasterReference);
            title({mapRequest.Layer.LayerTitle, timeStr}, ...
                'Interpreter', 'none', 'FontWeight', 'bold')
            frame = getframe(hFig);
        catch e
            fprintf(['Server error: %s.\n', ...
                'Ignoring frame number %d on day %s.\n'], ...
                e.message, k, timeStr)
  10. Read in all video frames.

    v = VideoReader(videoFilename);
    vidFrames = read(v);
    numFrames = get(v, 'NumberOfFrames');
  11. Create a MATLAB® movie structure from the video frames.

    frames = struct('cdata', [], 'colormap', []);
    frames(numFrames) = frames(1);
    for k = 1 : numFrames
        frames(k).cdata = vidFrames(:,:,:,k);
        frames(k).colormap = [];
  12. Playback movie once at the video's frame rate.

    movie(hFig, frames, 1, v.FrameRate)

Create Animated GIF File of WMS Maps

Read and display an animation of the Larsen Ice Shelf experiencing a dramatic collapse between January 31 and March 7, 2002.

  1. Search the WMS Database for the phrase "Larsen Ice Shelf."

    iceLayer = wmsfind('Larsen Ice Shelf');

    Try the first layer.

  2. Construct a WebMapServer object.

    server = WebMapServer(iceLayer(1).ServerURL);
  3. Use the WebMapServer.updateLayers method to synchronize the layer with the WMS source server. Retrieve the most recent data and fill in the Abstract, CoordRefSysCodes, and Details fields.

    iceLayer = server.updateLayers(iceLayer(1));
  4. View the abstract.

    fprintf('%s\n', iceLayer(1).Abstract)
  5. Create the WMSMapRequest object.

    request = WMSMapRequest(iceLayer(1), server);
  6. Because you have updated your layer, the Details field now has content. Click Details in the MATLAB Variables editor. Then, click Dimension. The name of the dimension is 'time'. Click Extent. The Extent field provides the available values for a dimension, in this case time. Save this information by entering the following at the command line:

    extent = [',' iceLayer.Details.Dimension.Extent, ','];
  7. Calculate the number of required frames. (The extent contains a comma before the first frame and after the last frame. To obtain the number of frames, subtract 1.)

    frameIndex = strfind(extent, ',');
    numFrames = numel(frameIndex) - 1;
  8. Open a figure window and set up a map axes with appropriate geographic limits.

    h = figure;
    worldmap(request.Latlim, request.Lonlim)
  9. Set the map axes properties. MLineLocation establishes the interval between displayed grid meridians. MLabelParallel determines the parallel where the labels appear.

    setm(gca,'MLineLocation', 1, 'MLabelLocation', 1, ...
       'MLabelParallel',-67.5, 'LabelRotation', 'off');
  10. Initialize the value of animated to 0.

    animated(1,1,1,numFrames) = 0;
  11. Display the image of the Larsen Ice Shelf on different days.

    for k=1:numFrames
       request.Time = extent(frameIndex(k)+1:frameIndex(k+1)-1);
       iceImage = server.getMap(request.RequestURL);
       geoshow(iceImage, request.RasterReference)
       title(request.Time, 'Interpreter', 'none')
       frame = getframe(h);
       if k == 1
          [animated, cmap] = rgb2ind(frame.cdata, 256, 'nodither');
          animated(:,:,1,k) = rgb2ind(frame.cdata, cmap, 'nodither');
  12. Save and then view the animated GIF file.

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

    Snapshot from Animation of Larsen Ice Shelf

    Courtesy NASA/Goddard Space Flight Center Scientific Visualization Studio

Animate Time-Lapse Radar Observations

Display Next-Generation Radar (NEXRAD) images for the United States using data from the Iowa Environmental Mesonet (IEM) Web map server. The server stores layers covering the past 50 minutes up to the present time in increments of 5 minutes. Read and display the merged layers.

  1. Find layers in the WMS Database that include 'mesonet' and 'nexrad' in their ServerURL fields.

    mesonet = wmsfind('mesonet*nexrad', 'SearchField', 'serverurl');
  2. NEXRAD Base Reflect Current ('nexrad-n0r') measures the intensity of precipitation. Refine your search to include only layers with this phrase in one of the search fields.

    nexrad = mesonet.refine('nexrad-n0r', 'SearchField', 'any');
  3. Remove the 900913 layers because they are intended for Google Maps™ overlay. Also remove the WMST layer because it contains data for different times.

    layers_900913 = nexrad.refine('900913', 'SearchField', ...
    layer_wmst = nexrad.refine('wmst', 'SearchField', 'layername'); 
    rmLayerNames = {layers_900913.LayerName layer_wmst.LayerName}; 
    index = ismember({nexrad.LayerName}, rmLayerNames); 
    nexrad = nexrad(~index); 
  4. Update your nexrad layer to fill in all fields and obtain most recent data.

    nexrad = wmsupdate(nexrad, 'AllowMultipleServers', true); 
  5. 'conus' represents the conterminous 48 U.S. states (all except Hawaii and Alaska). Use the usamap function to construct a map axes for the conterminous states. Read in the nexrad layers.

    region = 'conus';
    mstruct = gcm;
    latlim = mstruct.maplatlimit;
    lonlim = mstruct.maplonlimit;
    [A, R] = wmsread(nexrad, 'Latlim', latlim, 'Lonlim', lonlim);
  6. Display the NEXRAD merged layers map. Overlay with United States state boundary polygons.

    geoshow(A, R);
    geoshow('usastatehi.shp', 'FaceColor', 'none');
    title({'NEXRAD Radar Map', 'Merged Layers'});

  7. Loop through the sequence of time-lapse radar observations.

    hfig = figure;
    hstates = geoshow('usastatehi.shp', 'FaceColor', 'none');
    numFrames = numel(nexrad);
    frames = struct('cdata', [], 'colormap', []);
    frames(numFrames) = frames;
    hmap = [];
    frameIndex = 0;
    for k = numFrames:-1:1
       frameIndex = frameIndex + 1;
       [A, R] = wmsread(nexrad(k), 'Latlim', latlim, 'Lonlim', lonlim);
       hmap = geoshow(A, R);
       frames(frameIndex) = getframe(hfig);
  8. Create an array to write out as an animated GIF file.

    animated(1,1,1,numFrames) = 0;
    for k=1:numFrames
       if k == 1
          [animated, cmap] = rgb2ind(frames(k).cdata, 256, 'nodither');
          animated(:,:,1,k) = ...
             rgb2ind(frames(k).cdata, cmap, 'nodither');
  9. View the animated GIF file.

    filename = 'wmsnexrad.gif';
    imwrite(animated, cmap, filename, 'DelayTime', 1.5, ...
       'LoopCount', inf);

See Also

| |

Related Topics