Loren on the Art of MATLAB

May 6th, 2010

Visualizing the Gulf of Mexico Oil Slick Using Web Map Service (WMS)

Today I'd like to introduce a guest blogger, Kelly Luetkemeyer, who is a Mapping Toolbox software developer at The MathWorks. Kelly was the lead developer on the project to incorporate Web Map Service (WMS) into the Mapping Toolbox.

Contents

Introduction

Recently, the tragic accident in the Gulf of Mexico caused me to think if The MathWorks tools could be used to help visualize the oil slick using imagery from Web Map Service (WMS) servers. I'll show you how to use WMS functionality in Mapping Toolbox to obtain MODIS imagery of the region. The MODIS sensor is onboard NASA's Terra satellite. WMS is a standard developed by the Open Geospatial Consortium (OGC) that allows clients to access web servers to obtain maps rendered as an image. We have cataloged over 1000 WMS servers and their layers and stored that information into a database shipped with the Mapping Toolbox. I will use that database to find a server that provides MODIS imagery.

Obtain Daily Planet Layer

The NASA WMS server located at the Jet Propulsion Laboratory (JPL) in Pasadena provides a continuously updating composite of visual images from the MODIS sensor onboard the Terra satellite. I can obtain the stored layer's information with the wmsfind function and narrow the search using the refine method of the WMSLayer object. 

jpl = wmsfind('onearth.jpl', 'SearchField', 'serverurl');
daily_planet = jpl.refine('daily_planet');

Update the Layer and Display the Abstract

The information stored in the database is only a partial subset. I can update the layer with the latest information from the server by using the wmsupdate function. Once the layer is updated, I can display the layer's abstract information. 

daily_planet = wmsupdate(daily_planet);
disp(['<html><bf>' daily_planet.Abstract '</bf></html>'])
A contiunously updating composite of visual images from TERRA MODIS scenes, see http://modis.gsfc.nasa.gov for details about MODIS.
	This dataset is built local on the OnEarth server, it updates as soon as scenes are available, usually with a 6 to 24 hour delay from real time.
        Images are produced from MODIS scenes using the HDFLook application.
	Base resolution is 8 arcseconds per pixel. The WMS "time" dimension can be used to retrieve past data, by using the YYYY-MM-DD notation.

Define the Region-of-Interest

I define the latitude and longitude bounds for a region surrounding the oil slick and request pixels at 8 arcseconds. The time variable is used to request the image from April 29, 2010. 

latlim = [28 31];
lonlim = [-90 -87];
cellSize = 8/3600;  % 8 arcseconds
time  = '2010-04-29';

Obtain the Image

The wmsread function retrieves the image from the server. The output A is the image, R is the referencing matrix, and URL is a WMS URL that can be used in a browser to download the image. The image coordinates are in latitude and longitude. Please note that sometimes this WMS server is busy and may issue an error message due to server overloading. If this happens, please try again at a later time. 

[A, R, URL] = wmsread(daily_planet, 'Time', time, ...
    'CellSize', cellSize, 'Latlim', latlim, 'Lonlim', lonlim);

Display the Caption

I initially found the MODIS image on a NASA website at http://www.nasa.gov/topics/earth/features/oil-creep.html. I copied the caption from the website and placed it into a file. 

fid = fopen('caption.txt');
c = fread(fid,'char=>char');
fclose(fid);
disp(['<html><bf>' c' '</bf></html>'])
 On April 29, the MODIS image on the Terra satellite captured a wide-view natural-color image of the oil slick (outlined in white) just off the Louisiana coast. The oil slick appears as dull gray interlocking comma shapes, one opaque and the other nearly transparent. Sunglint -- the mirror-like reflection of the sun off the water -- enhances the oil slick’s visibility. The northwestern tip of the oil slick almost touches the Mississippi Delta. Credit: NASA/Earth Observatory/Jesse Allen, using data provided courtesy of the University of Wisconsin’s Space Science and Engineering Center MODIS Direct Broadcast system.

Display the Image

Since A is an RGB image, I can use the Image Processing Toolbox function imshow to display it. The coordinates of the image are in longitude (rows) and latitude (columns).

I can use the Mapping Toolbox function geoshow to reproject the image into a Lambert projection. With the Mapping Toolbox functions, I can also overlay geographic data onto the map, in this case, the boundaries of the states in the region-of-interest. The boundaries of the states are obtained from the shapefile included with the Mapping Toolbox. 

S = shaperead('usastatehi.shp', 'UseGeoCoords', true);
figure('Renderer','zbuffer')
usamap(A,R)
geoshow([S.Lat], [S.Lon]);
geoshow(A,R)
title({'MODIS Image of Gulf of Mexico Oil Slick', time, ...
    'Data Courtesy of NASA'});

View the Image at Full Resolution

You can view the full image in your browser by clicking on the following link: http://blogs.mathworks.com/images/loren/231/gom.png

View the Image Using the WMS URL

If you do not have the Mapping Toolbox, you can still use WMS if you know the WMS URL. The third output from wmsread provides this information. If you insert this URL into a browser, the browser will download the image from the WMS server. You can also use the MATLAB function, urlwrite(URL, filename) to download the image. 

disp(URL(1:52)), disp(URL(53:95)), disp(URL(96:140))
disp(URL(141:170)), disp(URL(171:212)), ...
disp(URL(213:262)), disp(URL(263:end))
http://onearth.jpl.nasa.gov/wms.cgi?TIME=2010-04-29&
SERVICE=WMS&LAYERS=daily_planet&EXCEPTIONS=
application/vnd.ogc.se_xml&FORMAT=image/jpeg&
TRANSPARENT=FALSE&HEIGHT=1350&
BGCOLOR=0xFFFFFF&REQUEST=GetMap&WIDTH=1350
&BBOX=-90.0,28.0,-87.0,31.0&STYLES=&SRS=EPSG:4326&
VERSION=1.1.1

Obtain May 1, 2010 Image

The imagery from April 30 through May 6, excluding May 1, is either missing, because the area was not imaged at that time by the sensor, or the region is cloudy. I'll download the May 1 imagery for comparison and see how the shape and size of the oil slick has changed. 

time  = '2010-05-01';
[A2, R] = wmsread(daily_planet, 'Time', time, 'CellSize', cellSize, ...
    'Latlim', latlim, 'Lonlim', lonlim);
figure
imshow(A2)

Enhance the Image

The image can be enhanced by using the adapthisteq function from the Image Processing Toolbox. 

B = A2;
for k=1:3
 B(:,:,k) = adapthisteq(A2(:,:,k));
end
figure
imshow(B)

Display Both Images

For comparison, I'll place both the April 29 and the enhanced May 1 imagery in the same Figure window. You can see how the shape and size of the oil slick has changed. 

figure('Renderer','zbuffer')
subplot(1,2,1)
usamap(A,R)
geoshow(A,R)
geoshow([S.Lat], [S.Lon]);
title({'2010-04-29', 'Data Courtesy of NASA'})

subplot(1,2,2)
usamap(B,R)
geoshow(B,R)
geoshow([S.Lat], [S.Lon]);
title({'2010-05-01', 'Data Courtesy of NASA'})

Using Google™ Earth with WMS

Some WMS servers support non-image outputs. The daily_planet.Details.ImageFormats variable contains all the output types that the JPL server supports. One output type is 'application/vnd.google-earth.kml+xml'. I will modify the URL to specify KML as the output. You can use this file, gom.kml, in Google™ Earth to view the imagery. 

kmlURL = strrep(URL, 'image/jpeg', ...
    'application/vnd.google-earth.kml+xml');
urlwrite(kmlURL, 'gom.kml');

Conclusion

I've shown how MATLAB, combined with functions from the Mapping Toolbox and Image Processing Toolbox, can be used to visualize and enhance data for a current event. Do you have any ideas how you can use these tools to follow or analyze a current event? Let me know here.


Get the MATLAB code

Published with MATLAB® 7.10

Author picture

By Loren Shure

20:46 UTC | Posted in Image Processing, Mapping, News | Permalink | 7 Comments »

You can follow any responses to this entry through the RSS 2.0 feed. Both comments and pings are currently closed.

7 Responses to “Visualizing the Gulf of Mexico Oil Slick Using Web Map Service (WMS)”

  1. Stuart Kozola replied on May 7th, 2010 at 14:30 UTC :

    I liked the post. Amazed you could find a day that wasn’t cloudy so you could see the spill. It also reminded me of Alan’s Webinar: http://www.mathworks.com/company/events/webinars/wbnr41540.html

    It has some cool animations and an oil spill simulation. Thought it was relevant for this topic. Here are the three demos he shows in the webinar:

    • Oil spill simulation: shapefile, map creation, numerical simulation
    • Weather avoidance: WMS, image processing, geographic calculations, navigation
    • Terrain analysis: digital elevation models, KML export, viewshed

  2. Andrew Symington replied on February 21st, 2011 at 16:24 UTC :

    Hey Loren

    Great article. Subsequently, JPL decided to remove the full WMS service and provide a tile serving in its place. Do you have updated example code for the revised service?

    Thanks
    Andrew

  3. Kelly Luetkemeyer replied on February 22nd, 2011 at 21:21 UTC :

    Hi Andrew,

    Thanks for pointing out that the data for this post is no longer available. Unfortunately, using the tiled service from JPL will not provide the spatial resolution needed for this region.

    One source for data replacement can be found from the USGS Emergency Operations Portal website at:

    http:://hdds.usgs.gov/arcgis/services/201004_OilSpill_GulfOfMexico/MapServer/WMSServer

    The data is organized by day listed in the LayerTitle property. The server does not store the data using the Time dimension. The data from the MODIS sensor ranges from 2010-05-01 to 2010-08-31.

    You can retrieve the May 1 layer as follows. I’ve broken up the server URL so that it will fit on the web page conveniently:

    server = 'http://hdds.usgs.gov/arcgis/services/';
product = '201004_OilSpill_GulfOfMexico';
serverURL = [server product '/MapServer/WMSServer'];
info = wmsinfo(serverURL)
layers = info.Layer;
layer = layers.refine('20100501*MODIS');
latlim = [28 31];
lonlim = [-90 -87];
cellSize = 8/3600;
[A, R, URL] = wmsread(layer, 'Latlim', latlim, 'Lonlim', lonlim,  ...
   'CellSize', cellSize);

figure
imshow(A)

    You can click on this link (the contents of the URL variable) to see the image retrieved from the server in your browser:

    http://hdds.usgs.gov/arcgis/services/201004_OilSpill_GulfOfMexico/MapServer/WMSServer?SERVICE=WMS&LAYERS=3&EXCEPTIONS=application/vnd.ogc.se_xml&FORMAT=image/jpeg&TRANSPARENT=FALSE&HEIGHT=1350&BGCOLOR=0xFFFFFF&REQUEST=GetMap&WIDTH=1350&BBOX=-90.0,28.0,-87.0,31.0&STYLES=&SRS=EPSG:4326&VERSION=1.1.1

    Notice that this image is saturated in the blue color. Again, you can enhance the image using the method outlined in the blog (there are of course other ways to enhance the image):

    B = A;
for k=1:3
 B(:,:,k) = adapthisteq(A(:,:,k));
end
figure
imshow(B)

    The enhanced data does not quite look the same as that from the JPL server. My assumption is that the data is not processed to the same level on the USGS server.

    The server provides quite a bit of data to examine. However, as with all WMS servers, if it doesn’t respond, you can either try again or check out the web page:

    http://hdds.usgs.gov/EO/library.php

    -Kelly

  4. Sohrab replied on March 30th, 2011 at 15:07 UTC :

    Dear Kelly

    I was trying to follow your example here for one of my project in which I am keen to use Matlab rather than any other software

    Could you kindly point me out where am i making mistake below.

    I desire to have SRTM data set of a particular region( a small one as described by latlim lonlim ) coordinates

    Also I was trying to use wmsfind but was unsuccessful in that case as well.

    Here is the code

    clear all
    close all
    clc

    dms = [-31 58 00; ...
    -31 57 30; ...
    -70 07 30; ...
    -70 07 00];
    format long g
    angleInDegrees = dms2degrees(dms);
    latlim = [angleInDegrees(1) angleInDegrees(2)];
    lonlim = [angleInDegrees(3) angleInDegrees(4)];

    server = ‘http://dds.cr.usgs.gov/srtm';
    product = ‘Mercedario’;
    serverURL = [server product '/MapServer/WMSServer'];
    info = wmsinfo(serverURL);
    layers = info.Layer;
    layer = layers.refine(‘srtm’);
    cellSize = 8/3600;
    [A, R, URL] = wmsread(layer, ‘Latlim’, latlim, ‘Lonlim’, lonlim, …
    ‘CellSize’, cellSize);

  5. Kelly Luetkemeyer replied on March 30th, 2011 at 17:09 UTC :

    Hi Sohrab,

    The USGS server listed above provides layers for emergeny operations and disaster awareness, but not elevation (SRTM) layers. However, the NASA WorldWind server does provide elevation layers. You can retrieve the data in actual meters. The NASA serverURL is:

    http://www.nasa.network.com/elev?

    Here is how you can modify your code:

    
layers = wmsfind('nasa.network*elev',  ...
   'SearchField', 'serverurl');
layers= wmsupdate(layers);
layer = layers.refine('srtm30', 'SearchField', 'layername');
cellSize =  0.0020833;
 [A, R, URL] = wmsread(layer, 'Latlim', latlim, ...
   'Lonlim', lonlim,  ...
   'CellSize', cellSize);

    A is in meters and is int16 class. Your requested area is quite small, A is only a 4-by-4 matrix.

    Hopefully this modification will help you.

    -Kelly

  6. hussain replied on December 28th, 2011 at 03:48 UTC :

    What about Aster GDEM maps in Matlab having size of 3601-by-3601 16 bit (30 m resolution)?

  7. Kelly Luetkemeyer replied on January 6th, 2012 at 14:37 UTC :

    Hi Hussain,

    You did not provide much information to me, but if you want to retrieve Global ASTER DEM data at 30 meter resolution and create an image that is 3601-by-3601 in size, then here is the code that you can use.

    First, calculate the requested cell size in degrees:

    numberOfCells = 3601;
    metersPerCell = 30;
    kmPerCell = metersPerCell / 1000;
    lengthInMeters = 3601 * 30;
    lengthInKM = lengthInMeters / 1000;
    radiusInKM = 6378.137;
    lengthInDegrees = km2deg(lengthInKM, radiusInKM);
    cellSize = km2deg(kmPerCell, radiusInKM);

    Next, choose a region of interest. I used the southern and western edge of the first region of interest in this article:

    latStart = 28;
    lonStart = -90;

    interval = [0 lengthInDegrees];
    latlim = [latStart latStart] + interval;
    lonlim = [lonStart lonStart] + interval;

    Next, obtain the layers from the server, as outlined above.

    layers = wmsfind(‘nasa.network*elev’, …
    ‘SearchField’, ‘serverurl’);
    layers= wmsupdate(layers);

    Choose any of these ASTER 30 meter layers:

    aster = layers.refine(‘aster’, ‘SearchField’, ‘any’);
    asterLayerNames = {aster.LayerName}’

    asterLayerNames =

    ‘EarthAsterElevations30m’
    ‘aster_30m’
    ‘mergedAsterElevations’

    Obtain and display the image:

    layerName = ‘aster_30m’;
    [A, R, URL] = wmsread(layer, ‘Latlim’, latlim, …
    ‘Lonlim’, lonlim, …
    ‘CellSize’, cellSize);

    figure
    geoshow(A,R,’DisplayType’,'texturemap’)

    Here is the size of A

    Name Size Bytes Class Attributes

    A 3601×3601 25934402 int16

    And here is the URL you can use if you want to download the image:

    http://www.nasa.network.com/elev?SERVICE=WMS&LAYERS=srtm30&CRS=EPSG:4326&FORMAT=image/bil&TRANSPARENT=FALSE&HEIGHT=3601&BGCOLOR=0xFFFFFF&REQUEST=GetMap&WIDTH=3601&BBOX=-90.0,28.0,-89.02954999856568,28.97045000143432&STYLES=&VERSION=1.3.0

    -Kelly
MathWorks

About

Loren Shure works on design of the MATLAB language at MathWorks. She writes here about once a week on MATLAB programming and related topics.

Subscribe

Blog Admin

These postings are the author's and don't necessarily represent the opinions of The MathWorks.