Documentation

Search the WMS Database for the MODIS layer.

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

Construct a WebMapServer object.

server = WebMapServer(modis.ServerURL);

Construct a WMSMapRequest object.

mapRequest = WMSMapRequest(modis, server);

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

Create an array indicating the first seven days.

days = 1:7;

Set the value of startTime to December 01, 2010 and use a serial date number.

time = '2010-12-01';
startTime = datenum(time);

Open a figure window with axes appropriate for the region specified by the modis layer.

hFig = figure('Color', 'white');
worldmap(mapRequest.Latlim, mapRequest.Lonlim);

Save each frame into a video file.

videoFilename = 'modis_dec.avi';
writer = VideoWriter(videoFilename);
writer.FrameRate = 1;
writer.Quality = 100;
writer.open;

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.RasterReference);
        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

Read in all video frames.

v = VideoReader(videoFilename);
vidFrames = read(v);
numFrames = get(v, 'NumberOfFrames');

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

Playback movie once at the video's frame rate.

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

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

iceLayer = wmsfind('Larsen Ice Shelf');

Try the first layer.

Construct a WebMapServer object.

server = WebMapServer(iceLayer(1).ServerURL);

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

View the abstract.

fprintf('%s\n', iceLayer(1).Abstract)

Create the WMSMapRequest object.

request = WMSMapRequest(iceLayer(1), server);

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, ','];

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;

Open a figure window and set up a map axes with appropriate geographic limits.

h = figure;
worldmap(request.Latlim, request.Lonlim)

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

Initialize the value of animated to 0.

animated(1,1,1,numFrames) = 0;

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')
   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

Save and then view the animated GIF file.

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

Courtesy NASA/Goddard Space Flight Center Scientific Visualization Studio

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

mesonet = wmsfind('mesonet*nexrad', 'SearchField', 'serverurl');

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

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

Update your nexrad layer to fill in all fields and obtain most recent data.

nexrad = wmsupdate(nexrad, 'AllowMultipleServers', true); 

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

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'});

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

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

View the animated GIF file.

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