Display Animation of Radar Images over GOES Backdrop

This example shows how to display NEXRAD radar images. The images cover 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.

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

wmst = wmsfind('nexrad-n0r-wmst', 'SearchField', 'layername'); 
wmst = wmsupdate(wmst); 

Find a generated CONUS composite of GOES IR imagery and update it.

goes = wmsfind('goes*conus*ir', 'SearchField', 'layername');
goes = wmsupdate(goes);

Create a figure with the desired geographic extent.

hfig = figure;
region = 'conus';
usamap(region)
mstruct = gcm;
latlim = mstruct.maplatlimit;
lonlim = mstruct.maplonlimit;

Read the GOES layer as a backdrop image.

cellsize = .1;
[backdrop, R] = wmsread(goes, 'ImageFormat', 'image/png', ...
   'Latlim', latlim, 'Lonlim', lonlim, 'Cellsize', cellsize);

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;

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.

borders = 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);
    uistack(borders,'top')
    title({wmst.LayerName, time})
    drawnow
    frames(k) = getframe(hfig);
    hour = hour + 1/24;
end

View the movie loop.

numTimes = 10;
fps = 1.5;
movie(hfig, frames, numTimes, fps);

See Also

| |

Related Topics