Animating Data Layers

Creating Movie of Terra/MODIS Images

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.

    modis.Details.Dimension.Extent
  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;
    writer.open;
    
  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
        try
            mapRequest.Time = startTime + k;
            timeStr = datestr(mapRequest.Time);
            dailyImage = server.getMap(mapRequest.RequestURL);
            geoshow(dailyImage, mapRequest.RasterRef);
            title({mapRequest.Layer.LayerTitle, timeStr}, ...
                'Interpreter', 'none', 'FontWeight', 'bold')
            shg
            frame = getframe(hFig);
            writer.writeVideo(frame);
        catch e
            fprintf(['Server error: %s.\n', ...
                'Ignoring frame number %d on day %s.\n'], ...
                e.message, k, timeStr)
        end
        drawnow
        shg
    end
    writer.close
    
  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 = [];
    end
    
  12. Playback movie once at the video's frame rate.

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

Creating an Animated GIF File

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.RasterRef)
       title(request.Time, 'Interpreter', 'none')
       drawnow
       shg
       frame = getframe(h);
       if k == 1
          [animated, cmap] = rgb2ind(frame.cdata, 256, 'nodither');
       else
          animated(:,:,1,k) = rgb2ind(frame.cdata, cmap, 'nodither');
       end
    end
  12. Save and then view the animated GIF file.

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

    Snapshot from Animation of Larsen Ice Shelf

Animating 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', ...
       'layername'); 
    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';
    figure
    usamap(region)
    mstruct = gcm;
    latlim = mstruct.maplatlimit;
    lonlim = mstruct.maplonlimit;
    [A, R] = wmsread(nexrad, 'Latlim', latlim, 'Lonlim', lonlim, ...
       'Transparent', true);
    
  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;
    usamap(region)
    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;
       delete(hmap)
       [A, R] = wmsread(nexrad(k), 'Latlim', latlim, 'Lonlim', lonlim);
       hmap = geoshow(A, R);
       uistack(hstates,'top')
       title(nexrad(k).LayerName)
       drawnow
       frames(frameIndex) = getframe(hfig);
    end
  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');
       else
          animated(:,:,1,k) = ...
             rgb2ind(frames(k).cdata, cmap, 'nodither');
       end     
    end
    
  9. View the animated GIF file.

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

Displaying Animation of Radar Images over GOES Backdrop

Display NEXRAD radar images for the past 24 hours, sampled at one-hour intervals, for the United States using data from the IEM WMS server. Use the JPL Daily Planet layer as the backdrop.

  1. Find the 'nexrad-n0r-wmst' layer and update it.

    wmst = wmsfind('nexrad-n0r-wmst', 'SearchField', 'layername'); 
    wmst = wmsupdate(wmst); 
    
  2. Find a generated CONUS composite of GOES IR imagery and update it.

    goes = wmsfind('goes*conus*ir', 'SearchField', 'layername');
    goes = wmsupdate(goes);
    
  3. Create a figure with the desired geographic extent.

    hfig = figure;
    region = 'conus';
    usamap(region)
    mstruct = gcm;
    latlim = mstruct.maplatlimit;
    lonlim = mstruct.maplonlimit;
    
  4. Read the GOES layer as a backdrop image.

    cellsize = .1;
    [backdrop, R] = wmsread(goes, 'ImageFormat', 'image/png', ...
       'Latlim', latlim, 'Lonlim', lonlim, 'Cellsize', cellsize);
    
  5. Calculate current time minus 24 hours and set up frames to hold the data from getframe.

    now_m24 = datestr(now-1);
    hour_m24 = [now_m24(1:end-5) '00:00'];
    hour = datenum(hour_m24);
    hmap = [];
    numFrames = 24;
    frames = struct('cdata', [], 'colormap', []);
    frames(numFrames) = frames;
    
  6. For each hour, obtain the hourly NEXRAD map data and combine it with a copy of the backdrop. Because of how this Web server handles PNG format, the resulting map data has an image with class double. Thus, you must convert it to uint8 before merging.

    geoshow('usastatehi.shp', 'FaceColor', 'none');
    black = [0,0,0];
    threshold = 0;
    for k=1:numFrames
        time = datestr(hour);
        [A, R] = wmsread(wmst, 'Latlim', latlim, 'Lonlim', lonlim, ...
            'Time', time, 'CellSize', cellsize, ...
            'BackgroundColor', black, 'ImageFormat', 'image/png');
        delete(hmap)
        index = any(A > threshold, 3);
        combination = backdrop;
        index = cat(3,index,index,index);
        combination(index) = uint8(255*A(index));
        hmap = geoshow(combination, R);
        title({wmst.LayerName, time})
        drawnow
        frames(k) = getframe(hfig);
        hour = hour + 1/24;
    end
    
  7. View the movie loop.

    numTimes = 10;
    fps = 1.5;
    movie(hfig, frames, numTimes, fps);
    
Was this topic helpful?