{"id":262,"date":"2011-01-20T13:55:55","date_gmt":"2011-01-20T13:55:55","guid":{"rendered":"https:\/\/blogs.mathworks.com\/loren\/2011\/01\/20\/tracking-a-hurricane-using-web-map-service-wms\/"},"modified":"2011-01-13T20:29:32","modified_gmt":"2011-01-13T20:29:32","slug":"tracking-a-hurricane-using-web-map-service-wms","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/loren\/2011\/01\/20\/tracking-a-hurricane-using-web-map-service-wms\/","title":{"rendered":"Tracking a Hurricane using Web Map Service (WMS)"},"content":{"rendered":"<div xmlns:mwsh=\"https:\/\/www.mathworks.com\/namespace\/mcode\/v1\/syntaxhighlight.dtd\" class=\"content\">\r\n   <introduction>\r\n      <p>Guest blogger, <a href=\"https:\/\/www.mathworks.com\/matlabcentral\/fileexchange\/authors\/18249\">Kelly Luetkemeyer<\/a>, who is a Mapping Toolbox software developer at The MathWorks, returns with another article about using and analyzing data\r\n         obtained from a <a href=\"https:\/\/blogs.mathworks.com\/loren\/2010\/05\/06\/oilslick\/\">Web Map Service (WMS)<\/a> server.\r\n      <\/p>\r\n   <\/introduction>\r\n   <h3>Contents<\/h3>\r\n   <div>\r\n      <ul>\r\n         <li><a href=\"#1\">Introduction<\/a><\/li>\r\n         <li><a href=\"#5\">Data Layer<\/a><\/li>\r\n         <li><a href=\"#6\">Internet Access<\/a><\/li>\r\n         <li><a href=\"#7\">Key Concepts<\/a><\/li>\r\n         <li><a href=\"#8\">Demo Function<\/a><\/li>\r\n         <li><a href=\"#10\">Setup: Define a Data Directory and Filename Utility Function<\/a><\/li>\r\n         <li><a href=\"#13\">Step  1: Find Katrina Layers From Local Database<\/a><\/li>\r\n         <li><a href=\"#18\">Step  2: Synchronize WMSLayer Object with Server<\/a><\/li>\r\n         <li><a href=\"#21\">Step  3: Explore Katrina Layer Details<\/a><\/li>\r\n         <li><a href=\"#27\">Step  4: Create Time Extent Variables<\/a><\/li>\r\n         <li><a href=\"#31\">Step  5: Retrieve Katrina Map from Server<\/a><\/li>\r\n         <li><a href=\"#33\">Step  6: Develop Algorithm to Calculate Region Properties<\/a><\/li>\r\n         <li><a href=\"#49\">Step  7: Display Katrina Map and Vector Overlays<\/a><\/li>\r\n         <li><a href=\"#53\">Step  8: Apply Algorithm to all Frames Using parfor<\/a><\/li>\r\n         <li><a href=\"#61\">Step  9: Create Animation Files<\/a><\/li>\r\n         <li><a href=\"#65\">Step 10: View Animated GIF File<\/a><\/li>\r\n         <li><a href=\"#67\">Step 11: Display Table of Time, Location, and Area<\/a><\/li>\r\n         <li><a href=\"#68\">Step 12: View Video File<\/a><\/li>\r\n         <li><a href=\"#70\">The makeTable Function<\/a><\/li>\r\n         <li><a href=\"#73\">The calculateRegionProperties Function<\/a><\/li>\r\n         <li><a href=\"#74\">The renderFrame Function<\/a><\/li>\r\n         <li><a href=\"#76\">Credits<\/a><\/li>\r\n         <li><a href=\"#77\">Can You Apply Parallel Processing to Your Processing in an Interesting Way?<\/a><\/li>\r\n      <\/ul>\r\n   <\/div>\r\n   <h3>Introduction<a name=\"1\"><\/a><\/h3>\r\n   <p>This demo shows how to track a hurricane using data from a WMS server. Typically web sites containing weather imagery may\r\n      contain a visual animation of GOES imagery for any given hurricane. For example,  a typical time-sequenced imagery of Katrina\r\n      from a web site hosting weather imagery is shown below in an animated GIF.\r\n   <\/p>\r\n   <p><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/262\/wmsgeographic.gif\"> <\/p>\r\n   <p>Note that this animation shows the data in an unprojected geographic coordinate system. However, you may wish to view the\r\n      data in a projected coordinate system, and overlay the hurricane track and the size of the hurricane, at any given point in\r\n      time.\r\n   <\/p>\r\n   <p>You can use MATLAB&reg; to analyze images from a WMS server and track a hurricane's position and area. This demo will show you\r\n      how to obtain the WMS imagery from a web site using features from the Mapping Toolbox&#8482;, analyze the imagery using features\r\n      from the Image Processing Toolbox&#8482;, and use <tt>parfor<\/tt> from the Parallel Computing Toolbox&#8482; to parallelize some computations for performance improvements.\r\n   <\/p>\r\n   <p>In particular, this demo will show you how to:<\/p>\r\n   <div>\r\n      <ul>\r\n         <li>obtain WMS imagery at a particular time-step from a web site,<\/li>\r\n         <li>view the imagery in a projected coordinate system,<\/li>\r\n         <li>overlay the hurricane's track,<\/li>\r\n         <li>use the <tt>parfor<\/tt> command to parallelize computations of the convex hull, boundary, and centroid of the hurricane,\r\n         <\/li>\r\n         <li>overlay those region properties onto the WMS basemap,<\/li>\r\n         <li>and animate the sequence.<\/li>\r\n      <\/ul>\r\n   <\/div>\r\n   <h3>Data Layer<a name=\"5\"><\/a><\/h3>\r\n   <p>The hurricane layer is from the NASA Goddard Space Flight Center's Scientific Visualization Studio (SVS) Image Server. The\r\n      data in this layer shows cloud patterns in the visible wavelength, from 0.52 to 0.72 microns, during hurricane Katrina from\r\n      August 24 through August 30 2005. The cloud data was extracted from the Geostationary Operational Environment Satellite (GOES-12)\r\n      imagery and overlaid on a background color image of the southeast United States.\r\n   <\/p>\r\n   <h3>Internet Access<a name=\"6\"><\/a><\/h3>\r\n   <p>Since WMS servers are located on the Internet, this demo can be set to access the Internet to dynamically render and retrieve\r\n      the maps from the WMS server or it can be set to use data previously retrieved from the Internet using the WMS capabilities\r\n      but now stored in local files.\r\n   <\/p>\r\n   <h3>Key Concepts<a name=\"7\"><\/a><\/h3>\r\n   <p>If you are new to WMS, several key concepts are important to understand and are listed here.<\/p>\r\n   <div>\r\n      <ul>\r\n         <li><i>Web Map Service<\/i>  &#8212; The Open Geospatial Consortium (OGC) defines a   Web Map Service (WMS) to be an entity that \"produces maps of spatially\r\n              referenced data dynamically from geographic information.\"\r\n         <\/li>\r\n         <li><i>WMS server<\/i> &#8212; A server that follows the guidelines of the OGC to   render maps and return them to clients\r\n         <\/li>\r\n         <li><i>map<\/i> &#8212; The OGC definition for map is \"a portrayal of geographic   information as a digital image file suitable for display on\r\n            a computer   screen.\"\r\n         <\/li>\r\n         <li><i>layer<\/i> &#8212; A dataset of a specific type of geographic information, such   as temperature, elevation, weather, orthophotos, boundaries,\r\n              demographics, topography, transportation, environmental measurements,   and various data from satellites\r\n         <\/li>\r\n      <\/ul>\r\n   <\/div>\r\n   <h3>Demo Function<a name=\"8\"><\/a><\/h3>\r\n   <p>The code shown in this demo can be found in this function:<\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\"><span style=\"color: #0000FF\">function<\/span> track_hurricane(useInternet)<\/pre><h3>Setup: Define a Data Directory and Filename Utility Function<a name=\"10\"><\/a><\/h3>\r\n   <p>This demo creates or reads several data files from a directory and uses the variable <tt>datadir<\/tt> to denote their location.\r\n   <\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">datadir = fullfile(pwd, <span style=\"color: #A020F0\">'KatrinaData'<\/span>);\r\n<span style=\"color: #0000FF\">if<\/span> ~exist(datadir, <span style=\"color: #A020F0\">'dir'<\/span>)\r\n   mkdir(datadir)\r\n<span style=\"color: #0000FF\">end<\/span><\/pre><p>Define an anonymous function to prepend <tt>datadir<\/tt> to the input filename:\r\n   <\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">datafile = @(filename)fullfile(datadir, filename);<\/pre><p>This demo uses data obtained from a Web Map Server and thus requires Internet access. However, you can store the data locally\r\n      the first time you run the demo and then set the <tt>useInternet<\/tt> flag to false. If the <tt>useInternet<\/tt> flag is not set, it defaults to false.\r\n   <\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\"><span style=\"color: #0000FF\">if<\/span> ~exist(<span style=\"color: #A020F0\">'useInternet'<\/span>, <span style=\"color: #A020F0\">'var'<\/span>)\r\n    useInternet = false;\r\n<span style=\"color: #0000FF\">end<\/span><\/pre><h3>Step  1: Find Katrina Layers From Local Database<a name=\"13\"><\/a><\/h3>\r\n   <p>One of the more challenging aspects of using WMS is finding a WMS server and then finding the layer that is of interest to\r\n      you. The process of finding a server that contains the data you need and constructing a specific and often complicated URL\r\n      with all the relevant details can be very daunting.\r\n   <\/p>\r\n   <p>The Mapping Toolbox&#8482; simplifies the process of locating WMS servers and layers by providing a local, installed, and pre-qualified\r\n      WMS database, that is search-able, using the function <tt>wmsfind<\/tt>. You can search the database for layers and servers that are of interest to you. Here is how you find layers containing the\r\n      term <tt>katrina<\/tt> in either the <tt>LayerName<\/tt> or <tt>LayerTitle<\/tt> field of the database:\r\n   <\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">katrina = wmsfind(<span style=\"color: #A020F0\">'katrina'<\/span>);\r\nwhos <span style=\"color: #A020F0\">katrina<\/span><\/pre><pre style=\"font-style:oblique\">  Name          Size            Bytes  Class       Attributes\r\n\r\n  katrina      95x1             47468  WMSLayer              \r\n\r\n<\/pre><p>The search for the term <tt>'katrina'<\/tt> returned a <tt>WMSLayer<\/tt> array containing over 90 layers. To inspect information about an individual layer, simply display it like this:\r\n   <\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">katrina(1)<\/pre><pre style=\"font-style:oblique\">ans = \r\n  WMSLayer\r\n\r\n  Properties:\r\n           Index: 1\r\n     ServerTitle: 'NASA SVS Image Server'\r\n       ServerURL: 'http:\/\/aes.gsfc.nasa.gov\/cgi-bin\/wms?'\r\n      LayerTitle: 'GOES-12 Imagery of Hurricane Katrina: Longwave Infrared Close-up (1024x1024 Animation)'\r\n       LayerName: '3216_22510'\r\n          Latlim: [15.0000 45.0000]\r\n          Lonlim: [-100.0000 -70.0000]\r\n\r\n<\/pre><p>If you type, <tt>katrina<\/tt>, in the command window, the entire contents of the array are displayed, with each element's index number included in the\r\n      output. This display makes it easy for you to examine the entire array quickly, searching for a layer of interest. You can\r\n      display only the <tt>LayerTitle<\/tt> property for each element by executing the command:\r\n   <\/p>\r\n   <p><tt>katrina.disp('Properties', 'layertitle', 'Index', 'off', 'Label', 'off');<\/tt><\/p>\r\n   <p>As you've discovered, a search for the generic word <tt>'katrina'<\/tt> returned results of over 90 layers and you need to select only one layer. In general, a search may even return thousands\r\n      of layers, which may be too large to review individually. Rather than searching the database again, you may refine your search\r\n      by using the <tt>refine<\/tt> method of the <tt>WMSLayer<\/tt> class. Using the <tt>refine<\/tt> method is more efficient and returns results faster than <tt>wmsfind<\/tt> since the search has already been narrowed to a smaller set.\r\n   <\/p>\r\n   <p>Supplying the query string, <tt>'goes-12*katrina*visible*close*up*animation'<\/tt>, to the <tt>refine<\/tt> method returns a <tt>WMSLayer<\/tt> array whose elements contain a match of the query string in either the <tt>LayerTitle<\/tt> or <tt>LayerName<\/tt> properties. The <tt>*<\/tt> character indicates a wild-card search. You can use the first layer found.\r\n   <\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">katrina = katrina.refine(<span style=\"color: #A020F0\">'goes-12*katrina*visible*close*up*animation'<\/span>);\r\nwhos <span style=\"color: #A020F0\">katrina<\/span><\/pre><pre style=\"font-style:oblique\">  Name         Size            Bytes  Class       Attributes\r\n\r\n  katrina      2x1              1032  WMSLayer              \r\n\r\n<\/pre><h3>Step  2: Synchronize WMSLayer Object with Server<a name=\"18\"><\/a><\/h3>\r\n   <p>The database only stores a subset of the layer information. For example, information from the layer's abstract, details about\r\n      the layer's attributes and style information, and the coordinate reference system of the layer are not returned by <tt>wmsfind<\/tt>. To return all the information, you need to use the <tt>wmsupdate<\/tt> function. <tt>wmsupdate<\/tt> synchronizes the layer from the database with the server, filling in the missing properties of the layer.\r\n   <\/p>\r\n   <p>Synchronize the first <tt>katrina<\/tt> layer with the server and display the abstract information. Since this action requires access to the Internet, wrap the call\r\n      around the <tt>useInternet<\/tt> flag.\r\n   <\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">cachefile = datafile(<span style=\"color: #A020F0\">'katrina.mat'<\/span>);\r\n<span style=\"color: #0000FF\">if<\/span> useInternet\r\n    katrina = wmsupdate(katrina(1));\r\n    <span style=\"color: #0000FF\">if<\/span> ~exist(cachefile, <span style=\"color: #A020F0\">'file'<\/span>)\r\n        save(cachefile, <span style=\"color: #A020F0\">'katrina'<\/span>)\r\n    <span style=\"color: #0000FF\">end<\/span>\r\n<span style=\"color: #0000FF\">else<\/span>\r\n    load(cachefile)\r\n<span style=\"color: #0000FF\">end<\/span><\/pre>\r\n\r\n<pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">disp([<span style=\"color: #A020F0\">'&lt;html&gt;'<\/span> katrina.Abstract <span style=\"color: #A020F0\">'&lt;\/html&gt;'<\/span>])<\/pre>\r\n\r\n<pre style=\"font-style:oblique\">The GOES-12 satellite sits at 75 degrees west longitude at an altitude of 36,000 kilometers over the equator, in geosynchronous orbit.  At this position its Imager instrument takes pictures of cloud patterns in several wavelengths for all of North and South America, a primary measurement used in weather forecasting.  The Imager takes a pattern of pictures of parts of the Earth in several wavelengths all day, measurements that are vital in weather forecasting.  This animation shows a daily sequence of GOES-12 images in the visible wavelengths, from 0.52 to 0.72 microns, during the period that Hurricane Katrina passed through the Gulf of Mexico.  At one kilometer resolution, the visible band measurement is the highest resolution data from the Imager, which accounts for the very high level of detail in these images.  For this animation, the cloud data was extracted from GOES image and laid over a background color image of the southeast United States.\r\n\r\nAdditional Credit:\r\nB&gt;Please give credit for this item to:&lt;\/b&gt;&lt;br \/&gt;\r\n<\/pre>\r\n\r\n<h3>Step  3: Explore Katrina Layer Details<a name=\"21\"><\/a><\/h3>\r\n   <p>You can find out more information about the <tt>katrina<\/tt> layer by exploring the <tt>Details<\/tt> property of the <tt>katrina<\/tt> layer. The <tt>Details.Attributes<\/tt> field informs you that the layer has fixed width and fixed height attributes, thus the size of the requested map cannot be\r\n      modified.\r\n   <\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">katrina.Details.Attributes<\/pre><pre style=\"font-style:oblique\">ans = \r\n      Queryable: 0\r\n       Cascaded: 0\r\n         Opaque: 1\r\n      NoSubsets: 1\r\n     FixedWidth: 1024\r\n    FixedHeight: 1024\r\n<\/pre><p>The <tt>Details.Dimension<\/tt> field informs you that the layer has a <tt>time<\/tt> dimension\r\n   <\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">katrina.Details.Dimension<\/pre><pre style=\"font-style:oblique\">ans = \r\n              Name: 'time'\r\n             Units: 'ISO8601'\r\n        UnitSymbol: ''\r\n           Default: '2005-08-30T17:45Z'\r\n    MultipleValues: 0\r\n      NearestValue: 0\r\n           Current: 0\r\n            Extent: '2005-08-23T17:45Z\/2005-08-30T17:45Z\/P1D'\r\n<\/pre><p>with an extent from <tt>2005-08-23T17:45Z<\/tt> to <tt>2005-08-30T17:45Z<\/tt> with a period of <tt>P1D<\/tt> (one day), as shown in the <tt>Details.Dimension.Extent<\/tt> field.\r\n   <\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">katrina.Details.Dimension.Extent<\/pre><pre style=\"font-style:oblique\">ans =\r\n2005-08-23T17:45Z\/2005-08-30T17:45Z\/P1D\r\n<\/pre><h3>Step  4: Create Time Extent Variables<a name=\"27\"><\/a><\/h3>\r\n   <p>Since the layer has a time dimension, you can extract the time extent.<\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">extent = katrina.Details.Dimension.Extent;\r\nslash = <span style=\"color: #A020F0\">'\/'<\/span>;\r\nslashIndex = strfind(extent, slash);\r\nstartTime = extent(1:slashIndex(1)-1);\r\nendTime = extent(slashIndex(1)+1:slashIndex(2)-1);<\/pre><p>Calculate numeric values for the start and end days. Note that the time extent is in yyyy-mm-dd format.<\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">hyphen = <span style=\"color: #A020F0\">'-'<\/span>;\r\nhyphenIndex = strfind(startTime, hyphen);\r\ndayIndex = [hyphenIndex(2) + 1, hyphenIndex(2) + 2];\r\nstartDay = str2double(startTime(dayIndex));\r\nendDay = str2double(endTime(dayIndex));<\/pre><p>The hurricane is not well-formed on the first day so start on the next day.<\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">startDay = startDay + 1;<\/pre><p>Create a requestTime variable to contain the start date.<\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">requestTime = startTime;\r\nrequestTime(dayIndex) = num2str(startDay);<\/pre><h3>Step  5: Retrieve Katrina Map from Server<a name=\"31\"><\/a><\/h3>\r\n   <p>Now that you have found a layer of interest, you can retrieve the raster map using the function <tt>wmsread<\/tt> and display the map using the function <tt>geoshow<\/tt>.  You can specify the requested time by setting the <tt>Time<\/tt> optional parameter.\r\n   <\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">cachefile = datafile(<span style=\"color: #A020F0\">'katrina.tif'<\/span>);\r\n<span style=\"color: #0000FF\">if<\/span> useInternet\r\n    [katrinaMap, R] = wmsread(katrina, <span style=\"color: #A020F0\">'Time'<\/span>, requestTime);\r\n    <span style=\"color: #0000FF\">if<\/span> ~exist(cachefile, <span style=\"color: #A020F0\">'file'<\/span>)\r\n        imwrite(katrinaMap, cachefile)\r\n        worldfilewrite(R, getworldfilename(cachefile))\r\n    <span style=\"color: #0000FF\">end<\/span>\r\n<span style=\"color: #0000FF\">else<\/span>\r\n    katrinaMap = imread(cachefile);\r\n    R = worldfileread(getworldfilename(cachefile));\r\n<span style=\"color: #0000FF\">end<\/span><\/pre><p>Display the <tt>katrinaMap<\/tt> and overlay the outlines of land areas worldwide and the United States using data from shape files shipped with the Mapping\r\n      Toolbox&#8482;.\r\n   <\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">states = shaperead(<span style=\"color: #A020F0\">'usastatehi.shp'<\/span>, <span style=\"color: #A020F0\">'UseGeoCoords'<\/span>, true);\r\nland = shaperead(<span style=\"color: #A020F0\">'landareas.shp'<\/span>, <span style=\"color: #A020F0\">'UseGeoCoords'<\/span>, true);\r\nfigure(<span style=\"color: #A020F0\">'Color'<\/span>, <span style=\"color: #A020F0\">'white'<\/span>, <span style=\"color: #A020F0\">'Renderer'<\/span>, <span style=\"color: #A020F0\">'zbuffer'<\/span>)\r\nusamap(katrina.Latlim, katrina.Lonlim)\r\ngeoshow(katrinaMap, R)\r\ngeoshow(states, <span style=\"color: #A020F0\">'FaceColor'<\/span>, <span style=\"color: #A020F0\">'none'<\/span>)\r\ngeoshow(land, <span style=\"color: #A020F0\">'FaceColor'<\/span>, <span style=\"color: #A020F0\">'none'<\/span>)\r\ntitle({katrina.LayerTitle, katrina.Details.Dimension.Default}, <span style=\"color: #0000FF\">...<\/span>\r\n    <span style=\"color: #A020F0\">'Interpreter'<\/span>, <span style=\"color: #A020F0\">'none'<\/span>, <span style=\"color: #A020F0\">'FontWeight'<\/span>, <span style=\"color: #A020F0\">'bold'<\/span>);<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/262\/track_hurricane_01.png\"> <h3>Step  6: Develop Algorithm to Calculate Region Properties<a name=\"33\"><\/a><\/h3>\r\n   <p>Using the imagery from the WMS server, calculate the hurricane's area and location of the center by assuming that the hurricane\r\n      is in the region with the most cloud cover (as defined in pixel units, not map units).\r\n   <\/p>\r\n   <p>The first step is to create a black and white image from the image of Katrina. As an initial guess, try using the <tt>graythresh<\/tt> function to calculate a threshold value.\r\n   <\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">thresh = graythresh(katrinaMap);\r\nkatrina_bw = im2bw(katrinaMap, thresh);\r\nfigure(<span style=\"color: #A020F0\">'Color'<\/span>, <span style=\"color: #A020F0\">'white'<\/span>, <span style=\"color: #A020F0\">'Renderer'<\/span>, <span style=\"color: #A020F0\">'zbuffer'<\/span>)\r\nusamap(katrina.Latlim, katrina.Lonlim)\r\ngeoshow(katrina_bw, R)\r\ntitle(<span style=\"color: #A020F0\">'Black and White Cloud Image Using graythresh'<\/span>)<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/262\/track_hurricane_02.png\"> <p>In the black and white image, you can see that there exists large tracts that are now white even though these areas do not\r\n      have significant cloud cover. One method to remove these regions is to set a higher threshold value.\r\n   <\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">factor = 1.9;\r\nthresh = factor*thresh;\r\nkatrina_bw = im2bw(katrinaMap, thresh);\r\nfigure(<span style=\"color: #A020F0\">'Color'<\/span>, <span style=\"color: #A020F0\">'white'<\/span>, <span style=\"color: #A020F0\">'Renderer'<\/span>, <span style=\"color: #A020F0\">'zbuffer'<\/span>)\r\nusamap(katrina.Latlim, katrina.Lonlim)\r\ngeoshow(katrina_bw, R)\r\ntitle(<span style=\"color: #A020F0\">'Black and White Cloud Image Using 1.9*graythresh'<\/span>)<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/262\/track_hurricane_03.png\"> <p>Calculate the convex hulls and areas of all regions using <tt>regionprops<\/tt>.\r\n   <\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">katrina_props = regionprops(katrina_bw, {<span style=\"color: #A020F0\">'ConvexHull'<\/span>, <span style=\"color: #A020F0\">'Area'<\/span>});<\/pre><p>Find the index that outlines the maximum area (in pixel units).<\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">[~, maxIndex] = max([katrina_props(:).Area]);\r\nmaxIndex = maxIndex(1);<\/pre><p>Convert the row and column indices where the area is at a maximum to latitude and longitude coordinates. Store the coordinates\r\n      into a structure.\r\n   <\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">row = katrina_props(maxIndex).ConvexHull(:,2);\r\ncol = katrina_props(maxIndex).ConvexHull(:,1);\r\nConvexHull = struct(<span style=\"color: #A020F0\">'Lat'<\/span>, [], <span style=\"color: #A020F0\">'Lon'<\/span>, [], <span style=\"color: #A020F0\">'Area'<\/span>, []);\r\n[ConvexHull.Lat, ConvexHull.Lon] = pix2latlon(R, row, col);<\/pre><p>Plot the convex hull onto the black and white cloud image.<\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">figure(<span style=\"color: #A020F0\">'Color'<\/span>, <span style=\"color: #A020F0\">'white'<\/span>, <span style=\"color: #A020F0\">'Renderer'<\/span>, <span style=\"color: #A020F0\">'zbuffer'<\/span>)\r\nusamap(katrina.Latlim, katrina.Lonlim)\r\ngeoshow(katrina_bw, R)\r\ngeoshow(ConvexHull.Lat, ConvexHull.Lon, <span style=\"color: #A020F0\">'Color'<\/span>,  <span style=\"color: #A020F0\">'cyan'<\/span>,  <span style=\"color: #A020F0\">'LineWidth'<\/span>,  2)\r\ntitle(<span style=\"color: #A020F0\">'Convex Hull Image'<\/span>)<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/262\/track_hurricane_04.png\"> <p>Calculate the area of the convex hull in units of square statute miles.<\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">ellipsoid = almanac(<span style=\"color: #A020F0\">'earth'<\/span>,<span style=\"color: #A020F0\">'wgs84'<\/span>,<span style=\"color: #A020F0\">'statutemiles'<\/span>);\r\nConvexHull.Area = areaint(ConvexHull.Lat, ConvexHull.Lon, ellipsoid);<\/pre><p>Calculate the boundaries of all objects using <tt>bwboundaries<\/tt>.\r\n   <\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">[B, L] = bwboundaries(katrina_bw);<\/pre><p>Convert the label matrix to an RGB image and display it.<\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">labelRGB = label2rgb(L, @jet, [.5 .5 .5]);\r\nfigure(<span style=\"color: #A020F0\">'Color'<\/span>, <span style=\"color: #A020F0\">'white'<\/span>, <span style=\"color: #A020F0\">'Renderer'<\/span>, <span style=\"color: #A020F0\">'zbuffer'<\/span>)\r\nusamap(katrina.Latlim, katrina.Lonlim)\r\ngeoshow(labelRGB, R)\r\ntitle(<span style=\"color: #A020F0\">'Boundary Image'<\/span>)<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/262\/track_hurricane_05.png\"> <p>Use a histogram to select the region with the most cloud cover.<\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">clouds = L(L~=0);\r\n[freq, loc] = hist(clouds ,unique(L));\r\n[~, maxIndex] = max(freq);\r\nhurricane = loc(maxIndex);<\/pre><p>Capture indices of the largest cloud cover.<\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">indices = B{hurricane};\r\nrow = indices(:,1);\r\ncol = indices(:,2);<\/pre><p>Convert indices of the largest cloud cover to latitude and longitude.<\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">Boundary = struct(<span style=\"color: #A020F0\">'Lat'<\/span>, [], <span style=\"color: #A020F0\">'Lon'<\/span>, [], <span style=\"color: #A020F0\">'Area'<\/span>, []);\r\n[Boundary.Lat, Boundary.Lon] = pix2latlon(R, row, col);<\/pre><p>Plot the boundary onto the labeled image.<\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">figure(<span style=\"color: #A020F0\">'Color'<\/span>, <span style=\"color: #A020F0\">'white'<\/span>, <span style=\"color: #A020F0\">'Renderer'<\/span>, <span style=\"color: #A020F0\">'zbuffer'<\/span>)\r\nusamap(katrina.Latlim, katrina.Lonlim)\r\ngeoshow(labelRGB, R)\r\ngeoshow(Boundary.Lat, Boundary.Lon, <span style=\"color: #A020F0\">'Color'<\/span>, <span style=\"color: #A020F0\">'white'<\/span>, <span style=\"color: #A020F0\">'LineWidth'<\/span>, 2)\r\ntitle(<span style=\"color: #A020F0\">'Boundary Image with Outline'<\/span>)<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/262\/track_hurricane_06.png\"> <p>Calculate the area of the boundary in units of square statute miles.<\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">Boundary.Area = areaint(Boundary.Lat, Boundary.Lon, ellipsoid);<\/pre><p>Find the center of the hurricane.<\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">katrina_center = regionprops(L, <span style=\"color: #A020F0\">'Centroid'<\/span>);\r\ncenter_row = katrina_center(hurricane).Centroid(:,2);\r\ncenter_col = katrina_center(hurricane).Centroid(:,1);\r\nCentroid = struct(<span style=\"color: #A020F0\">'Lat'<\/span>, [], <span style=\"color: #A020F0\">'Lon'<\/span>, []);\r\n[Centroid.Lat, Centroid.Lon] = pix2latlon(R, center_row, center_col);<\/pre><h3>Step  7: Display Katrina Map and Vector Overlays<a name=\"49\"><\/a><\/h3>\r\n   <p>Import hurricane track vector data. The vector dataset, <tt>katrina.shp<\/tt>, contains the historical track for hurricane Katrina. The track file contains the 6-hourly (0000, 0600, 1200, 1800 UTC) center\r\n      locations and intensities of the hurricane.\r\n   <\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">tracks = shaperead(<span style=\"color: #A020F0\">'KatrinaData\/katrina.shp'<\/span>, <span style=\"color: #A020F0\">'UseGeo'<\/span>, true);\r\ntracks(1)<\/pre><pre style=\"font-style:oblique\">ans = \r\n       Geometry: 'Line'\r\n    BoundingBox: [2x2 double]\r\n            Lon: [-75.1 -75.7 NaN]\r\n            Lat: [23.1 23.4 NaN]\r\n           Year: 2005\r\n          Month: 8\r\n            Day: 23\r\n           Name: 'KATRINA_2005'\r\n     Wind_Knots: 30\r\n       Pressure: 1008\r\n            Cat: 'TD'\r\n<\/pre><p>Display the WMS basemap and vector overlays:<\/p>\r\n   <div>\r\n      <ul>\r\n         <li>Overlay the outlines of worldwide land areas.<\/li>\r\n         <li>Overlay the outline of the United States.<\/li>\r\n         <li>Overlay the outline of the convex hull in red.<\/li>\r\n         <li>Overlay the outline of the boundary of the hurricane in cyan.<\/li>\r\n         <li>Overlay the center of the hurricane in black.<\/li>\r\n         <li>Overlay the eye of the hurricane track in yellow.<\/li>\r\n      <\/ul>\r\n   <\/div><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">figure(<span style=\"color: #A020F0\">'Color'<\/span>, <span style=\"color: #A020F0\">'white'<\/span>, <span style=\"color: #A020F0\">'Renderer'<\/span>, <span style=\"color: #A020F0\">'zbuffer'<\/span>)\r\nusamap(katrina.Latlim, katrina.Lonlim)\r\ngeoshow(katrinaMap, R)\r\ngeoshow(land, <span style=\"color: #A020F0\">'FaceColor'<\/span>, <span style=\"color: #A020F0\">'none'<\/span>)\r\ngeoshow(states, <span style=\"color: #A020F0\">'FaceColor'<\/span>, <span style=\"color: #A020F0\">'none'<\/span>)\r\ngeoshow(ConvexHull.Lat, ConvexHull.Lon, <span style=\"color: #0000FF\">...<\/span>\r\n    <span style=\"color: #A020F0\">'LineWidth'<\/span>, 4, <span style=\"color: #A020F0\">'Color'<\/span>, <span style=\"color: #A020F0\">'red'<\/span>);\r\ngeoshow(Boundary.Lat, Boundary.Lon, <span style=\"color: #0000FF\">...<\/span>\r\n    <span style=\"color: #A020F0\">'LineWidth'<\/span>, 2, <span style=\"color: #A020F0\">'Color'<\/span>, <span style=\"color: #A020F0\">'cyan'<\/span>);\r\ngeoshow(Centroid.Lat, Centroid.Lon, <span style=\"color: #0000FF\">...<\/span>\r\n    <span style=\"color: #A020F0\">'LineWidth'<\/span>, 4, <span style=\"color: #A020F0\">'Color'<\/span>, <span style=\"color: #A020F0\">'black'<\/span>, <span style=\"color: #A020F0\">'Marker'<\/span>, <span style=\"color: #A020F0\">'+'<\/span>)\r\ngeoshow(tracks, <span style=\"color: #A020F0\">'LineWidth'<\/span>, 2, <span style=\"color: #A020F0\">'Color'<\/span>, <span style=\"color: #A020F0\">'yellow'<\/span>)\r\ntitle({katrina.LayerTitle, katrina.Details.Dimension.Default}, <span style=\"color: #0000FF\">...<\/span>\r\n    <span style=\"color: #A020F0\">'Interpreter'<\/span>, <span style=\"color: #A020F0\">'none'<\/span>, <span style=\"color: #A020F0\">'FontWeight'<\/span>, <span style=\"color: #A020F0\">'bold'<\/span>);<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/262\/track_hurricane_07.png\"> <p>Print the area statistics.<\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">fprintf(<span style=\"color: #A020F0\">'Area of the convex hull is %0.2f square statute miles.\\n'<\/span>, <span style=\"color: #0000FF\">...<\/span>\r\n    ConvexHull.Area);\r\nfprintf(<span style=\"color: #A020F0\">'Area of the boundary    is %0.2f square statute miles.\\n'<\/span>, <span style=\"color: #0000FF\">...<\/span>\r\n    Boundary.Area);<\/pre><pre style=\"font-style:oblique\">Area of the convex hull is 81463.28 square statute miles.\r\nArea of the boundary    is 61953.07 square statute miles.\r\n<\/pre><h3>Step  8: Apply Algorithm to all Frames Using parfor<a name=\"53\"><\/a><\/h3>\r\n   <p>Now that you have developed an algorithm to track and calculate the hurricane's size and position, the next step is to apply\r\n      that algorithm to a time sequence of data from the server. Use the <tt>parfor<\/tt> command to distribute the calculations to multiple processors. It is part of the MATLAB&reg; language, and behaves essentially\r\n      like a regular <tt>for<\/tt>-loop if you do not have access to the Parallel Computing Toolbox&#8482; product.\r\n   <\/p>\r\n   <p>Calculate the total number of animation frames.<\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">numFrames = endDay - startDay + 1;<\/pre><p>Assign the initial requestTime.<\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">requestTime = startTime;<\/pre><p>Initialize variables to hold the data calculated from each frame.<\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">time = cell(numFrames, 1);\r\nS = struct(<span style=\"color: #A020F0\">'Lat'<\/span>, [], <span style=\"color: #A020F0\">'Lon'<\/span>, [], <span style=\"color: #A020F0\">'Area'<\/span>, []);\r\nConvexHull = S;\r\nBoundary = S;\r\nConvexHull(numFrames) = ConvexHull;\r\nBoundary(numFrames) = Boundary;\r\nCentroid = struct(<span style=\"color: #A020F0\">'Lat'<\/span>, [], <span style=\"color: #A020F0\">'Lon'<\/span>, []);\r\nCentroid(numFrames) = Centroid;\r\nframeName = datafile(<span style=\"color: #A020F0\">'frame_'<\/span>);<\/pre><p>Check the status of the MATLAB&reg; pool and determine if it is open. Use the MATLAB&reg; pool to allow the body of the <tt>parfor<\/tt> loop to run in parallel.\r\n   <\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">poolSize = matlabpool(<span style=\"color: #A020F0\">'size'<\/span>);\r\n<span style=\"color: #0000FF\">if<\/span> poolSize == 0\r\n    fprintf(<span style=\"color: #A020F0\">'The MATLAB pool is not open, the loop will run sequentially.'<\/span>);\r\n<span style=\"color: #0000FF\">else<\/span>\r\n    fprintf(<span style=\"color: #A020F0\">'Using %d workers.\\n'<\/span>, poolSize);\r\n<span style=\"color: #0000FF\">end<\/span><\/pre><pre style=\"font-style:oblique\">Using 4 workers.\r\n<\/pre><p>For each day, from <tt>startDay<\/tt> to <tt>endDay<\/tt>, perform the following operations:\r\n   <\/p>\r\n   <div>\r\n      <ul>\r\n         <li>obtain the raster map for the time frame from either the web map server   or from the cache file,<\/li>\r\n         <li>calculate the convex hull, boundary, and their areas,<\/li>\r\n         <li>calculate the location of the center of the hurricane,<\/li>\r\n         <li>render the data into a Figure window, and<\/li>\r\n         <li>save the data from the graphics handles to a MAT-file.<\/li>\r\n      <\/ul>\r\n   <\/div>\r\n   <p>To create the animation frames, you can load the rendered data from each saved MAT-file and capture the data to an image.\r\n      You need to do this operation outside of the <tt>parfor<\/tt> loop since each compute lab is started in a headless mode. By storing each rendered frame into a MAT-file within the <tt>parfor<\/tt> loop, you can improve performance since rendering the data into a Figure window from a MAT-file is a faster operation than\r\n      the display operations for the vector and raster data.\r\n   <\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">tworkers = zeros(1, poolSize);\r\ntic;\r\n<span style=\"color: #0000FF\">parfor<\/span> k=1:numFrames\r\n    tstart = tic;\r\n    <span style=\"color: #228B22\">% Update currentDay.<\/span>\r\n    currentDay = startDay + k -1;\r\n\r\n    <span style=\"color: #228B22\">% Assign time value for this frame.<\/span>\r\n    frameTime = requestTime;\r\n    frameTime(dayIndex) = num2str(currentDay);\r\n    time{k} = frameTime;\r\n\r\n    <span style=\"color: #228B22\">% Get the WMS map from the server (or file) for this time period.<\/span>\r\n    cachefile = datafile([<span style=\"color: #A020F0\">'katrina_'<\/span> num2str(currentDay) <span style=\"color: #A020F0\">'.tif'<\/span>]);\r\n    <span style=\"color: #0000FF\">if<\/span> useInternet\r\n         RGB = wmsread(katrina, <span style=\"color: #A020F0\">'Time'<\/span>, frameTime);\r\n         <span style=\"color: #0000FF\">if<\/span> ~exist(cachefile, <span style=\"color: #A020F0\">'file'<\/span>)\r\n            imwrite(RGB, cachefile);\r\n         <span style=\"color: #0000FF\">end<\/span>\r\n    <span style=\"color: #0000FF\">else<\/span>\r\n        RGB = imread(cachefile);\r\n    <span style=\"color: #0000FF\">end<\/span>\r\n\r\n    <span style=\"color: #228B22\">% Calculate the convex hull, boundary, areas, and centroid using the<\/span>\r\n    <span style=\"color: #228B22\">% same threshold value as calculated for the first image.<\/span>\r\n    [ConvexHull(k), Boundary(k), Centroid(k)] <span style=\"color: #0000FF\">...<\/span>\r\n        = calculateRegionProperties(RGB, R, thresh);\r\n\r\n    <span style=\"color: #228B22\">% Render the frame.<\/span>\r\n    hFig = figure(<span style=\"color: #A020F0\">'Color'<\/span>, <span style=\"color: #A020F0\">'white'<\/span>, <span style=\"color: #A020F0\">'Renderer'<\/span>, <span style=\"color: #A020F0\">'zbuffer'<\/span>);\r\n    renderFrame(frameTime, katrina, RGB,  R,  <span style=\"color: #0000FF\">...<\/span>\r\n        ConvexHull(k), Boundary(k), Centroid(k), land, states, tracks);\r\n\r\n    <span style=\"color: #228B22\">% Save the rendered frame.<\/span>\r\n    hAxes = get(hFig, <span style=\"color: #A020F0\">'CurrentAxes'<\/span>);\r\n    frameFilename = [frameName num2str(k) <span style=\"color: #A020F0\">'.fig'<\/span>];\r\n    hgsave(hAxes, frameFilename);\r\n    close(hFig);\r\n\r\n     <span style=\"color: #228B22\">% Save the workers elapsed time.<\/span>\r\n    tworkers(k) = toc(tstart);\r\n\r\n    <span style=\"color: #228B22\">% Print statistics<\/span>\r\n    fprintf(<span style=\"color: #A020F0\">'\\nFrame %d elapsed time: %4.8f\\n'<\/span>, k, tworkers(k));\r\n    fprintf(<span style=\"color: #A020F0\">'Frame %d convex hull area: %0.2f.\\n'<\/span>, <span style=\"color: #0000FF\">...<\/span>\r\n        k, ConvexHull(k).Area);\r\n    fprintf(<span style=\"color: #A020F0\">'Frame %d boundary area:    %0.2f.\\n'<\/span>, <span style=\"color: #0000FF\">...<\/span>\r\n        k, Boundary(k).Area);\r\n<span style=\"color: #0000FF\">end<\/span>\r\ntloop = toc;<\/pre><pre style=\"font-style:oblique\">\r\nFrame 4 elapsed time: 17.75099091\r\nFrame 4 convex hull area: 179144.66.\r\nFrame 4 boundary area:    113185.27.\r\n\r\nFrame 2 elapsed time: 19.09949585\r\nFrame 2 convex hull area: 192105.82.\r\nFrame 2 boundary area:    114881.50.\r\n\r\nFrame 6 elapsed time: 25.17652566\r\nFrame 6 convex hull area: 443757.43.\r\nFrame 6 boundary area:    303722.05.\r\n\r\nFrame 5 elapsed time: 26.48190417\r\nFrame 5 convex hull area: 161777.52.\r\nFrame 5 boundary area:    137392.76.\r\n\r\nFrame 3 elapsed time: 18.83270332\r\nFrame 3 convex hull area: 149839.57.\r\nFrame 3 boundary area:    108983.98.\r\n\r\nFrame 1 elapsed time: 18.44629997\r\nFrame 1 convex hull area: 81463.28.\r\nFrame 1 boundary area:    61953.07.\r\n\r\nFrame 7 elapsed time: 15.59668967\r\nFrame 7 convex hull area: 248729.73.\r\nFrame 7 boundary area:    199408.60.\r\n<\/pre><p>You expect the total time of all the workers to be more than the loop time when the loop is executed in parallel.<\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">fprintf(<span style=\"color: #A020F0\">'Total time for all workers: %4.8f\\n'<\/span>, sum(tworkers));\r\nfprintf(<span style=\"color: #A020F0\">'Total loop time:            %4.8f\\n'<\/span>, tloop);<\/pre><pre style=\"font-style:oblique\">Total time for all workers: 141.38460955\r\nTotal loop time:            43.90086599\r\n<\/pre><h3>Step  9: Create Animation Files<a name=\"61\"><\/a><\/h3>\r\n   <p>An animation can be viewed in the browser when the browser opens an animated GIF file or an AVI video file. To create the\r\n      animation frames of the WMS basemap and vector overlays, render the data from the saved graphics files, capture the rendered\r\n      data into a <tt>frame<\/tt> structure, and save <tt>frame<\/tt> into either an array for the animated GIF file or to an AVI video file.\r\n   <\/p>\r\n   <p>To share with others or to post to web video services, you can create an AVI video file containing all the frames using the\r\n      <tt>VideoWriter<\/tt> class.\r\n   <\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">basename = <span style=\"color: #A020F0\">'wmstrack'<\/span>;\r\nvideoFilename = [<span style=\"color: #A020F0\">'html\/'<\/span> basename <span style=\"color: #A020F0\">'.avi'<\/span>];\r\n<span style=\"color: #0000FF\">if<\/span> exist(videoFilename, <span style=\"color: #A020F0\">'file'<\/span>)\r\n    delete(videoFilename)\r\n<span style=\"color: #0000FF\">end<\/span>\r\nwriter = VideoWriter(videoFilename);\r\nwriter.FrameRate = 1;\r\nwriter.Quality = 100;\r\nwriter.open;<\/pre><p>For each frame, read the graphics from the cache file created in the body of the <tt>parfor<\/tt> loop, and display the WMS basemap and vector overlays. Capture each frame and store into the frames array. The loop cannot\r\n      run in parallel since the video frames must be ordered.\r\n   <\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">hFig = figure(<span style=\"color: #A020F0\">'Color'<\/span>, <span style=\"color: #A020F0\">'white'<\/span>, <span style=\"color: #A020F0\">'Renderer'<\/span>, <span style=\"color: #A020F0\">'zbuffer'<\/span>);\r\n<span style=\"color: #0000FF\">for<\/span> k=1:numFrames\r\n\r\n    <span style=\"color: #228B22\">% Delete the previous frame.<\/span>\r\n    hAxes = get(hFig, <span style=\"color: #A020F0\">'CurrentAxes'<\/span>);\r\n    delete(hAxes)\r\n\r\n    <span style=\"color: #228B22\">% Load the basemap and vector overlays into the axes.<\/span>\r\n    frameFilename = [frameName num2str(k) <span style=\"color: #A020F0\">'.fig'<\/span>];\r\n    hgload(frameFilename);\r\n\r\n    <span style=\"color: #228B22\">% Capture the current frame.<\/span>\r\n    shg\r\n    frame = getframe(hFig);\r\n\r\n    <span style=\"color: #228B22\">% Store the frame into the |animated| array for the GIF file.<\/span>\r\n    <span style=\"color: #0000FF\">if<\/span> k == 1\r\n        [animated, cmap] = rgb2ind(frame.cdata, 256, <span style=\"color: #A020F0\">'nodither'<\/span>);\r\n        animated(:,:,1,numFrames) = animated;\r\n    <span style=\"color: #0000FF\">else<\/span>\r\n        animated(:,:,1,k) = rgb2ind(frame.cdata, cmap, <span style=\"color: #A020F0\">'nodither'<\/span>);\r\n    <span style=\"color: #0000FF\">end<\/span>\r\n\r\n    <span style=\"color: #228B22\">% Write the frame to the AVI file.<\/span>\r\n    writer.writeVideo(frame);\r\n<span style=\"color: #0000FF\">end<\/span>\r\n<span style=\"color: #228B22\">% Close the Figure window and AVI file.<\/span>\r\nclose\r\nwriter.close;<\/pre><p>Write the animated GIF file.<\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">filename = [basename <span style=\"color: #A020F0\">'.gif'<\/span>];\r\ndelayTime = 2.0;\r\nimwrite(animated, cmap, filename, <span style=\"color: #A020F0\">'DelayTime'<\/span>, delayTime, <span style=\"color: #A020F0\">'LoopCount'<\/span>, inf);<\/pre><h3>Step 10: View Animated GIF File<a name=\"65\"><\/a><\/h3>\r\n   <p>An animation can be viewed in the browser when the browser opens an animated GIF file.<\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">web(filename)<\/pre><p><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/262\/wmstrack.gif\"> <\/p>\r\n  \r\n <h3>Step 11: Display Table of Time, Location, and Area<a name=\"67\"><\/a><\/h3><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">headers = {<span style=\"color: #A020F0\">'Time'<\/span>,  <span style=\"color: #A020F0\">'Latitude'<\/span>,  <span style=\"color: #A020F0\">'Longitude'<\/span>,  <span style=\"color: #0000FF\">...<\/span>\r\n    <span style=\"color: #A020F0\">'Convex Hull Area'<\/span>, <span style=\"color: #A020F0\">'Boundary Area'<\/span>};\r\nmakeTable(headers, time, [Centroid.Lat], [Centroid.Lon],  <span style=\"color: #0000FF\">...<\/span>\r\n    [ConvexHull.Area], [Boundary.Area]);<\/pre>\r\n\r\n<table border=1 cellpadding=1><tr> <td align=\"center\"> <b> Time <\/b><\/td><td align=\"center\"> <b> Latitude <\/b><\/td><td align=\"center\"> <b> Longitude <\/b><\/td><td align=\"center\"> <b> Convex Hull Area <\/b><\/td><td align=\"center\"> <b> Boundary Area <\/b><\/td><\/tr> <tr><td align=\"right\"> 2005-08-24T17:45Z <\/td><td align=\"right\"> 25.399 <\/td><td align=\"right\"> -76.144 <\/td><td align=\"right\"> 81463.282 <\/td><td align=\"right\"> 61953.066 <\/td><\/tr><tr><td align=\"right\"> 2005-08-25T17:45Z <\/td><td align=\"right\"> 25.913 <\/td><td align=\"right\"> -79.282 <\/td><td align=\"right\"> 192105.821 <\/td><td align=\"right\"> 114881.500 <\/td><\/tr><tr><td align=\"right\"> 2005-08-26T17:45Z <\/td><td align=\"right\"> 23.568 <\/td><td align=\"right\"> -82.217 <\/td><td align=\"right\"> 149839.573 <\/td><td align=\"right\"> 108983.980 <\/td><\/tr><tr><td align=\"right\"> 2005-08-27T17:45Z <\/td><td align=\"right\"> 23.680 <\/td><td align=\"right\"> -84.421 <\/td><td align=\"right\"> 179144.655 <\/td><td align=\"right\"> 113185.267 <\/td><\/tr><tr><td align=\"right\"> 2005-08-28T17:45Z <\/td><td align=\"right\"> 26.285 <\/td><td align=\"right\"> -88.482 <\/td><td align=\"right\"> 161777.517 <\/td><td align=\"right\"> 137392.758 <\/td><\/tr><tr><td align=\"right\"> 2005-08-29T17:45Z <\/td><td align=\"right\"> 34.192 <\/td><td align=\"right\"> -87.477 <\/td><td align=\"right\"> 443757.435 <\/td><td align=\"right\"> 303722.052 <\/td><\/tr><tr><td align=\"right\"> 2005-08-30T17:45Z <\/td><td align=\"right\"> 37.301 <\/td><td align=\"right\"> -85.280 <\/td><td align=\"right\"> 248729.729 <\/td><td align=\"right\"> 199408.599 <\/td><\/tr><\/table>\r\n\r\n\r\n<h3>Step 12: View Video File<a name=\"68\"><\/a><\/h3>\r\n   <p>The following <a href=\"https:\/\/blogs.mathworks.com\/images\/loren\/262\/wmstrack.avi\">wmstrack.avi<\/a> video is generated and can be viewed using an external viewer.\r\n   <\/p>\r\n  \r\n\r\n\r\n<h3>The makeTable Function<a name=\"70\"><\/a><\/h3><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">dbtype <span style=\"color: #A020F0\">makeTable<\/span><\/pre><pre style=\"font-style:oblique\">\r\n1     function makeTable(headers, time, lat, lon, area1, area2)\r\n2     %makeTable Make HTML table\r\n3     %\r\n4     %   makeTable(headers, time, lat, lon, area1, area2) generates and displays\r\n5     %   (using disp) HTML to visualize a table of input values. HEADERS is a\r\n6     %   cell array of names for the top row. The values of the table are in the\r\n7     %   arrays, TIME, LAT, LON, AREA1, and AREA2.\r\n8     \r\n9     % Copyright 2010-2011 The MathWorks, Inc.\r\n10    \r\n11    header = ' ';\r\n12    for k=1:numel(headers)\r\n13        tableEntry = sprintf( ...\r\n14            '&lt;td align=\"center\"&gt; &lt;b&gt; %s &lt;\/b&gt;&lt;\/td&gt;', ...\r\n15            headers{k});\r\n16        header = [header tableEntry]; %#ok&lt;AGROW&gt;\r\n17    end\r\n18    header = ['&lt;tr&gt;' header '&lt;\/tr&gt;'];\r\n19    \r\n20    table = ' ';\r\n21    for k = 1:numel(time)\r\n22        tableEntry = sprintf( ...\r\n23            ['&lt;tr&gt;', ...\r\n24            '&lt;td align=\"right\"&gt; %s &lt;\/td&gt;', ...\r\n25            '&lt;td align=\"right\"&gt; %0.3f &lt;\/td&gt;', ...\r\n26            '&lt;td align=\"right\"&gt; %0.3f &lt;\/td&gt;', ...\r\n27            '&lt;td align=\"right\"&gt; %0.3f &lt;\/td&gt;', ...\r\n28            '&lt;td align=\"right\"&gt; %0.3f &lt;\/td&gt;&lt;\/tr&gt;'], ...\r\n29            time{k}, lat(k), lon(k), area1(k), area2(k));\r\n30        table = [table tableEntry]; %#ok&lt;AGROW&gt;\r\n31    end\r\n32    \r\n33    disp([ ...\r\n34        '&lt;html&gt;&lt;table border=1 cellpadding=1&gt;', ...\r\n35       header, ...\r\n36        table, ...\r\n37        '&lt;\/table&gt;&lt;\/html&gt;']);\r\n\r\n<\/pre>\r\n<pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\"><span style=\"color: #0000FF\">end<\/span><\/pre><h3>The calculateRegionProperties Function<a name=\"73\"><\/a><\/h3>\r\n   <p>This function calculates the convex hull, BW boundary, area, and centroid properties of an RGB image.<\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\"><span style=\"color: #0000FF\">function<\/span> [ConvexHull, Boundary, Centroid] = <span style=\"color: #0000FF\">...<\/span>\r\n    calculateRegionProperties(RGB, R, thresh)\r\n\r\n<span style=\"color: #228B22\">% Calculate the convex hull and area of the objects in the RGB image.<\/span>\r\nbw = im2bw(RGB, thresh);\r\nprops = regionprops(bw, {<span style=\"color: #A020F0\">'ConvexHull'<\/span>, <span style=\"color: #A020F0\">'Area'<\/span>});\r\n\r\n<span style=\"color: #228B22\">% Find the indices that outline the maximum area.<\/span>\r\n[~, maxIndex] = max([props(:).Area]);\r\n\r\n<span style=\"color: #228B22\">% Convert the row and column indices where the area is at a maximum to<\/span>\r\n<span style=\"color: #228B22\">% latitude and longitude coordinates. Store the coordinates into a<\/span>\r\n<span style=\"color: #228B22\">% structure.<\/span>\r\nrow = props(maxIndex).ConvexHull(:,2);\r\ncol = props(maxIndex).ConvexHull(:,1);\r\nConvexHull = struct(<span style=\"color: #A020F0\">'Lat'<\/span>, [], <span style=\"color: #A020F0\">'Lon'<\/span>, [], <span style=\"color: #A020F0\">'Area'<\/span>, []);\r\n[ConvexHull.Lat, ConvexHull.Lon] = pix2latlon(R, row, col);\r\n\r\n<span style=\"color: #228B22\">% Calculate the area of the convex hull in units of square statute miles.<\/span>\r\nellipsoid = almanac(<span style=\"color: #A020F0\">'earth'<\/span>,<span style=\"color: #A020F0\">'wgs84'<\/span>,<span style=\"color: #A020F0\">'statutemiles'<\/span>);\r\nConvexHull.Area = areaint(ConvexHull.Lat, ConvexHull.Lon, ellipsoid);\r\n\r\n<span style=\"color: #228B22\">% Calculate the boundaries of all objects.<\/span>\r\n[B, L] = bwboundaries(bw);\r\n\r\n<span style=\"color: #228B22\">% Use a histogram to select the largest region.<\/span>\r\nobjects = L(L~=0);\r\n[freq, loc] = hist(objects ,unique(L));\r\n[~, maxIndex] = max(freq);\r\nlargestRegion = loc(maxIndex);\r\n\r\n<span style=\"color: #228B22\">% Capture indices of the largest region.<\/span>\r\nindices = B{largestRegion};\r\nrow = indices(:,1);\r\ncol = indices(:,2);\r\n\r\n<span style=\"color: #228B22\">% Convert indices of the largest region cover to latitude and longitude.<\/span>\r\nBoundary = struct(<span style=\"color: #A020F0\">'Lat'<\/span>, [], <span style=\"color: #A020F0\">'Lon'<\/span>, [], <span style=\"color: #A020F0\">'Area'<\/span>, []);\r\n[Boundary.Lat, Boundary.Lon] = pix2latlon(R, row, col);\r\n\r\n<span style=\"color: #228B22\">% Calculate the area of the boundary in units of square statute miles.<\/span>\r\nBoundary.Area = areaint(Boundary.Lat, Boundary.Lon, ellipsoid);\r\n\r\n<span style=\"color: #228B22\">% Find the center of the largest region.<\/span>\r\ncenter = regionprops(L, <span style=\"color: #A020F0\">'Centroid'<\/span>);\r\ncenter_row = center(largestRegion).Centroid(:,2);\r\ncenter_col = center(largestRegion).Centroid(:,1);\r\nCentroid = struct(<span style=\"color: #A020F0\">'Lat'<\/span>, [], <span style=\"color: #A020F0\">'Lon'<\/span>, []);\r\n[Centroid.Lat, Centroid.Lon] = pix2latlon(R, center_row, center_col);\r\n<span style=\"color: #0000FF\">end<\/span><\/pre><h3>The renderFrame Function<a name=\"74\"><\/a><\/h3><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\"><span style=\"color: #0000FF\">function<\/span> renderFrame(frameTime, layer, layerMap, R,  <span style=\"color: #0000FF\">...<\/span>\r\n    ConvexHull, Boundary, Centroid, land, states, tracks)<\/pre><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\"><span style=\"color: #228B22\">% This function renders the frame data into a Figure window.<\/span>\r\n\r\nusamap(layer.Latlim, layer.Lonlim)\r\ngeoshow(layerMap, R)\r\ngeoshow(land, <span style=\"color: #A020F0\">'FaceColor'<\/span>, <span style=\"color: #A020F0\">'none'<\/span>)\r\ngeoshow(states, <span style=\"color: #A020F0\">'FaceColor'<\/span>, <span style=\"color: #A020F0\">'none'<\/span>)\r\ngeoshow(ConvexHull.Lat, ConvexHull.Lon, <span style=\"color: #A020F0\">'LineWidth'<\/span>, 4, <span style=\"color: #A020F0\">'Color'<\/span>, <span style=\"color: #A020F0\">'red'<\/span>);\r\ngeoshow(Boundary.Lat, Boundary.Lon, <span style=\"color: #A020F0\">'LineWidth'<\/span>, 2, <span style=\"color: #A020F0\">'Color'<\/span>, <span style=\"color: #A020F0\">'cyan'<\/span>);\r\ngeoshow(Centroid.Lat, Centroid.Lon,   <span style=\"color: #A020F0\">'LineWidth'<\/span>,  4,  <span style=\"color: #A020F0\">'Color'<\/span>, <span style=\"color: #A020F0\">'k'<\/span>,  <span style=\"color: #A020F0\">'Marker'<\/span>, <span style=\"color: #A020F0\">'+'<\/span>)\r\ngeoshow(tracks, <span style=\"color: #A020F0\">'LineWidth'<\/span>, 2, <span style=\"color: #A020F0\">'Color'<\/span>,  <span style=\"color: #A020F0\">'yellow'<\/span>)\r\ntitle({layer.LayerTitle, frameTime}, <span style=\"color: #A020F0\">'Interpreter'<\/span>, <span style=\"color: #A020F0\">'none'<\/span>, <span style=\"color: #A020F0\">'FontWeight'<\/span>, <span style=\"color: #A020F0\">'bold'<\/span>);<\/pre><h3>Credits<a name=\"76\"><\/a><\/h3>\r\n   <p><i>Katrina Layer<\/i><\/p>\r\n   <p>The Katrina layer used in the demo is from the NASA Goddard Space Flight Center's SVS Image Server and is maintained by the\r\n      Scientific Visualization Studio.\r\n   <\/p>\r\n   <p>For more information about this server, run:<\/p><pre>  &gt;&gt; wmsinfo('http:\/\/aes.gsfc.nasa.gov\/cgi-bin\/wms?')<\/pre><p><i>usastatehi.shp<\/i><\/p>\r\n   <p>This vector dataset contains moderate-resolution outlines and names for the 50 United States plus the District of Columbia,\r\n      stored in shapefile format with coordinates in latitude (Y) and longitude (X) in degrees. It is based on U.S. Census Bureau's\r\n      TIGER thinned boundary files. The data file is shipped with the Mapping Toolbox&#8482;.\r\n   <\/p>\r\n   <p><i>landareas.shp<\/i><\/p>\r\n   <p>This vector dataset contains low-resolution outlines of land areas worldwide, stored as polygons in shapefile format with\r\n      coordinates in latitude (Y) and longitude (X) in degrees. The raw coordinate data were derived from the Digital Chart of the\r\n      World (DCW) browser layer, published by U.S. National Geospatial-Intelligence Agency (NGA), formerly National Imagery and\r\n      Mapping Agency (NIMA). The data file is shipped with the Mapping Toolbox&#8482;.\r\n   <\/p>\r\n   <p><i>katrina.shp<\/i><\/p>\r\n   <p>This vector dataset contains the historical track for hurricane Katrina. The raw data were obtained from the atl_hurtrack.shp\r\n      file from the National Oceanic and Atmospheric Administration Coastal Services Center. This Historical North Atlantic Tropical\r\n      Cyclone Tracks file contains the 6-hourly (0000, 0600, 1200, 1800 UTC) center locations and intensities for all subtropical\r\n      depressions and storms, extra-tropical storms, tropical lows, waves, disturbances, depressions and storms, and all hurricanes,\r\n      from 1851 through 2009. These data are intended for geographic display and analysis at the national level, and for large regional\r\n      areas. The data should be displayed and analyzed at scales appropriate for 1:2,000,000-scale data. No responsibility is assumed\r\n      by the National Oceanic and Atmospheric Administration in the use of these data. The complete track data was downloaded from\r\n      NOAA. The data file is not shipped with the Mapping Toolbox&#8482;.\r\n   <\/p>\r\n   <h3>Can You Apply Parallel Processing to Your Processing in an Interesting Way?<a name=\"77\"><\/a><\/h3>\r\n   <p>As you can see here, I've used <tt>parfor<\/tt> to independently fetch and process each frame before constructing the final video.  Is that technique something you might\r\n      be able to take advantage of?  Let us know <a href=\"https:\/\/blogs.mathworks.com\/loren\/?p=262\">here<\/a>.\r\n   <\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\"><span style=\"color: #0000FF\">end<\/span><\/pre><script language=\"JavaScript\">\r\n<!--\r\n\r\n    function grabCode_5cb939f7486747759f02d0a247201b95() {\r\n        \/\/ Remember the title so we can use it in the new page\r\n        title = document.title;\r\n\r\n        \/\/ Break up these strings so that their presence\r\n        \/\/ in the Javascript doesn't mess up the search for\r\n        \/\/ the MATLAB code.\r\n        t1='5cb939f7486747759f02d0a247201b95 ' + '##### ' + 'SOURCE BEGIN' + ' #####';\r\n        t2='##### ' + 'SOURCE END' + ' #####' + ' 5cb939f7486747759f02d0a247201b95';\r\n    \r\n        b=document.getElementsByTagName('body')[0];\r\n        i1=b.innerHTML.indexOf(t1)+t1.length;\r\n        i2=b.innerHTML.indexOf(t2);\r\n \r\n        code_string = b.innerHTML.substring(i1, i2);\r\n        code_string = code_string.replace(\/REPLACE_WITH_DASH_DASH\/g,'--');\r\n\r\n        \/\/ Use \/x3C\/g instead of the less-than character to avoid errors \r\n        \/\/ in the XML parser.\r\n        \/\/ Use '\\x26#60;' instead of '<' so that the XML parser\r\n        \/\/ doesn't go ahead and substitute the less-than character. \r\n        code_string = code_string.replace(\/\\x3C\/g, '\\x26#60;');\r\n\r\n        author = 'Kelly Luetkemeyer';\r\n        copyright = 'Copyright 2011 The MathWorks, Inc.';\r\n\r\n        w = window.open();\r\n        d = w.document;\r\n        d.write('<pre>\\n');\r\n        d.write(code_string);\r\n\r\n        \/\/ Add author and copyright lines at the bottom if specified.\r\n        if ((author.length > 0) || (copyright.length > 0)) {\r\n            d.writeln('');\r\n            d.writeln('%%');\r\n            if (author.length > 0) {\r\n                d.writeln('% _' + author + '_');\r\n            }\r\n            if (copyright.length > 0) {\r\n                d.writeln('% _' + copyright + '_');\r\n            }\r\n        }\r\n\r\n        d.write('<\/pre>\\n');\r\n      \r\n      d.title = title + ' (MATLAB code)';\r\n      d.close();\r\n      }   \r\n      \r\n-->\r\n<\/script><p style=\"text-align: right; font-size: xx-small; font-weight:lighter;   font-style: italic; color: gray\"><br><a href=\"javascript:grabCode_5cb939f7486747759f02d0a247201b95()\"><span style=\"font-size: x-small;        font-style: italic;\">Get \r\n            the MATLAB code \r\n            <noscript>(requires JavaScript)<\/noscript><\/span><\/a><br><br>\r\n      Published with MATLAB&reg; 7.11<br><\/p>\r\n<\/div>\r\n<!--\r\n5cb939f7486747759f02d0a247201b95 ##### SOURCE BEGIN #####\r\n%% Tracking a Hurricane using Web Map Service (WMS) \r\n% Guest blogger, <https:\/\/www.mathworks.com\/matlabcentral\/fileexchange\/authors\/18249 Kelly Luetkemeyer>, who is a Mapping Toolbox\r\n% software developer at The MathWorks, returns with another article about\r\n% using and analyzing data obtained from a <https:\/\/blogs.mathworks.com\/loren\/2010\/05\/06\/oilslick\/ Web Map Service (WMS)> server.\r\n\r\n%% Introduction\r\n% This demo shows how to track a hurricane using data from a WMS server.\r\n% Typically web sites containing weather imagery may contain a visual\r\n% animation of GOES imagery for any given hurricane. For example,  a\r\n% typical time-sequenced imagery of Katrina from a web site hosting weather\r\n% imagery is shown below in an animated GIF.\r\n\r\n%%\r\n% <<wmsgeographic.gif>>\r\n\r\n%%\r\n% Note that this animation shows the data in an unprojected geographic\r\n% coordinate system. However, you may wish to view the data in a projected\r\n% coordinate system, and overlay the hurricane track and the size of the\r\n% hurricane, at any given point in time. \r\n\r\n%%\r\n% You can use MATLAB(R) to analyze images from a WMS server and track a\r\n% hurricane's position and area. This demo will show you how to obtain the\r\n% WMS imagery from a web site using features from the Mapping Toolbox\u00e2\u201e\u00a2,\r\n% analyze the imagery using features from the Image Processing Toolbox\u00e2\u201e\u00a2,\r\n% and use |parfor| from the Parallel Computing Toolbox\u00e2\u201e\u00a2 to parallelize some\r\n% computations for performance improvements.\r\n%\r\n% In particular, this demo will show you how to:\r\n%\r\n% * obtain WMS imagery at a particular time-step from a web site,\r\n% * view the imagery in a projected coordinate system, \r\n% * overlay the hurricane's track, \r\n% * use the |parfor| command to parallelize computations of the convex\r\n% hull, boundary, and centroid of the hurricane,\r\n% * overlay those region properties onto the WMS basemap,\r\n% * and animate the sequence.\r\n\r\n%% Data Layer\r\n% The hurricane layer is from the NASA Goddard Space Flight Center's\r\n% Scientific Visualization Studio (SVS) Image Server. The data in this\r\n% layer shows cloud patterns in the visible wavelength, from 0.52 to 0.72\r\n% microns, during hurricane Katrina from August 24 through August 30 2005.\r\n% The cloud data was extracted from the Geostationary Operational\r\n% Environment Satellite (GOES-12) imagery and overlaid on a background\r\n% color image of the southeast United States.\r\n\r\n%% Internet Access\r\n% Since WMS servers are located on the Internet, this demo can be set to\r\n% access the Internet to dynamically render and retrieve the maps from the\r\n% WMS server or it can be set to use data previously retrieved from the \r\n% Internet using the WMS capabilities but now stored in local files.\r\n\r\n%% Key Concepts\r\n% If you are new to WMS, several key concepts are important to understand\r\n% and are listed here.\r\n% \r\n% * _Web Map Service_  \u00e2\u20ac\u201d The Open Geospatial Consortium (OGC) defines a \r\n%   Web Map Service (WMS) to be an entity that \"produces maps of spatially\r\n%   referenced data dynamically from geographic information.\"\r\n% * _WMS server_ \u00e2\u20ac\u201d A server that follows the guidelines of the OGC to \r\n%   render maps and return them to clients\r\n% * _map_ \u00e2\u20ac\u201d The OGC definition for map is \"a portrayal of geographic \r\n%   information as a digital image file suitable for display on a computer \r\n%   screen.\"\r\n% * _layer_ \u00e2\u20ac\u201d A dataset of a specific type of geographic information, such \r\n%   as temperature, elevation, weather, orthophotos, boundaries, \r\n%   demographics, topography, transportation, environmental measurements, \r\n%   and various data from satellites\r\n\r\n%% Demo Function\r\n% The code shown in this demo can be found in this function:\r\nfunction track_hurricane(useInternet)\r\n\r\n%% Setup: Define a Data Directory and Filename Utility Function\r\n% This demo creates or reads several data files from a directory and uses\r\n% the variable |datadir| to denote their location.\r\ndatadir = fullfile(pwd, 'KatrinaData');\r\nif ~exist(datadir, 'dir')\r\n   mkdir(datadir)\r\nend\r\n\r\n%%\r\n% Define an anonymous function to prepend |datadir| to the input filename:\r\ndatafile = @(filename)fullfile(datadir, filename);\r\n\r\n%%\r\n% This demo uses data obtained from a Web Map Server and thus requires\r\n% Internet access. However, you can store the data locally the first time\r\n% you run the demo and then set the |useInternet| flag to false. If the\r\n% |useInternet| flag is not set, it defaults to false.\r\nif ~exist('useInternet', 'var') \r\n    useInternet = false;\r\nend\r\n\r\n%% Step  1: Find Katrina Layers From Local Database\r\n% One of the more challenging aspects of using WMS is finding a WMS server\r\n% and then finding the layer that is of interest to you. The process of\r\n% finding a server that contains the data you need and constructing a\r\n% specific and often complicated URL with all the relevant details can be \r\n% very daunting. \r\n%\r\n% The Mapping Toolbox\u00e2\u201e\u00a2 simplifies the process of locating WMS servers and\r\n% layers by providing a local, installed, and pre-qualified WMS database,\r\n% that is search-able, using the function |wmsfind|. You can search the\r\n% database for layers and servers that are of interest to you. Here is how\r\n% you find layers containing the term |katrina| in either the |LayerName| \r\n% or |LayerTitle| field of the database:\r\nkatrina = wmsfind('katrina');\r\nwhos katrina\r\n\r\n%%\r\n% The search for the term |'katrina'| returned a |WMSLayer| array\r\n% containing over 90 layers. To inspect information about an individual\r\n% layer, simply display it like this:\r\nkatrina(1)\r\n\r\n%%\r\n% If you type, |katrina|, in the command window, the entire contents of the\r\n% array are displayed, with each element's index number included in the\r\n% output. This display makes it easy for you to examine the entire array\r\n% quickly, searching for a layer of interest. You can display only the\r\n% |LayerTitle| property for each element by executing the command:\r\n% \r\n% |katrina.disp('Properties', 'layertitle', 'Index', 'off', 'Label', 'off');|\r\n\r\n%%\r\n% As you've discovered, a search for the generic word |'katrina'| returned\r\n% results of over 90 layers and you need to select only one layer. In\r\n% general, a search may even return thousands of layers, which may be too\r\n% large to review individually. Rather than searching the database again,\r\n% you may refine your search by using the |refine| method of the |WMSLayer|\r\n% class. Using the |refine| method is more efficient and returns results\r\n% faster than |wmsfind| since the search has already been narrowed to a\r\n% smaller set. \r\n\r\n%%\r\n% Supplying the query string, |'goes-12*katrina*visible*close*up*animation'|, \r\n% to the |refine| method returns a |WMSLayer| array whose elements contain\r\n% a match of the query string in either the |LayerTitle| or |LayerName|\r\n% properties. The |*| character indicates a wild-card search. You can use\r\n% the first layer found.\r\nkatrina = katrina.refine('goes-12*katrina*visible*close*up*animation');\r\nwhos katrina\r\n\r\n%% Step  2: Synchronize WMSLayer Object with Server\r\n% The database only stores a subset of the layer information. For example,\r\n% information from the layer's abstract, details about the layer's\r\n% attributes and style information, and the coordinate reference system of\r\n% the layer are not returned by |wmsfind|. To return all the information,\r\n% you need to use the |wmsupdate| function. |wmsupdate| synchronizes the\r\n% layer from the database with the server, filling in the missing\r\n% properties of the layer.\r\n\r\n%%\r\n% Synchronize the first |katrina| layer with the server and display the\r\n% abstract information. Since this action requires access to the Internet,\r\n% wrap the call around the |useInternet| flag.\r\ncachefile = datafile('katrina.mat');\r\nif useInternet\r\n    katrina = wmsupdate(katrina(1));\r\n    if ~exist(cachefile, 'file')\r\n        save(cachefile, 'katrina')\r\n    end\r\nelse\r\n    load(cachefile)\r\nend\r\n\r\n%%\r\ndisp(['<html>' katrina.Abstract '<\/html>'])\r\n\r\n%% Step  3: Explore Katrina Layer Details\r\n% You can find out more information about the |katrina| layer by exploring\r\n% the |Details| property of the |katrina| layer. The |Details.Attributes|\r\n% field informs you that the layer has fixed width and fixed height\r\n% attributes, thus the size of the requested map cannot be modified.\r\n\r\n%%\r\nkatrina.Details.Attributes\r\n\r\n%%\r\n% The |Details.Dimension| field informs you that the layer has a |time|\r\n% dimension \r\n\r\n%%\r\nkatrina.Details.Dimension\r\n\r\n%%\r\n% with an extent from |2005-08-23T17:45Z| to |2005-08-30T17:45Z|\r\n% with a period of |P1D| (one day), as shown in the\r\n% |Details.Dimension.Extent| field. \r\n\r\n%%\r\nkatrina.Details.Dimension.Extent\r\n\r\n%% Step  4: Create Time Extent Variables\r\n% Since the layer has a time dimension, you can extract the time extent.\r\nextent = katrina.Details.Dimension.Extent;\r\nslash = '\/';\r\nslashIndex = strfind(extent, slash);\r\nstartTime = extent(1:slashIndex(1)-1);\r\nendTime = extent(slashIndex(1)+1:slashIndex(2)-1);\r\n\r\n%%\r\n% Calculate numeric values for the start and end days. \r\n% Note that the time extent is in yyyy-mm-dd format.\r\nhyphen = '-';\r\nhyphenIndex = strfind(startTime, hyphen);\r\ndayIndex = [hyphenIndex(2) + 1, hyphenIndex(2) + 2];\r\nstartDay = str2double(startTime(dayIndex));\r\nendDay = str2double(endTime(dayIndex));\r\n\r\n%%\r\n% The hurricane is not well-formed on the first day so start on the next\r\n% day.\r\nstartDay = startDay + 1;\r\n\r\n%%\r\n% Create a requestTime variable to contain the start date.\r\nrequestTime = startTime;\r\nrequestTime(dayIndex) = num2str(startDay);\r\n\r\n\r\n%% Step  5: Retrieve Katrina Map from Server\r\n% Now that you have found a layer of interest, you can retrieve the raster\r\n% map using the function |wmsread| and display the map using the function\r\n% |geoshow|.  You can specify the requested time by setting the |Time| \r\n% optional parameter.\r\ncachefile = datafile('katrina.tif');\r\nif useInternet\r\n    [katrinaMap, R] = wmsread(katrina, 'Time', requestTime);\r\n    if ~exist(cachefile, 'file')\r\n        imwrite(katrinaMap, cachefile)\r\n        worldfilewrite(R, getworldfilename(cachefile))\r\n    end\r\nelse\r\n    katrinaMap = imread(cachefile);\r\n    R = worldfileread(getworldfilename(cachefile));\r\nend\r\n\r\n%%\r\n% Display the |katrinaMap| and overlay the outlines of land areas worldwide\r\n% and the United States using data from shape files shipped with the\r\n% Mapping Toolbox\u00e2\u201e\u00a2.\r\nstates = shaperead('usastatehi.shp', 'UseGeoCoords', true);\r\nland = shaperead('landareas.shp', 'UseGeoCoords', true);\r\nfigure('Color', 'white', 'Renderer', 'zbuffer')\r\nusamap(katrina.Latlim, katrina.Lonlim)\r\ngeoshow(katrinaMap, R)\r\ngeoshow(states, 'FaceColor', 'none')\r\ngeoshow(land, 'FaceColor', 'none')\r\ntitle({katrina.LayerTitle, katrina.Details.Dimension.Default}, ...\r\n    'Interpreter', 'none', 'FontWeight', 'bold');\r\n\r\n%% Step  6: Develop Algorithm to Calculate Region Properties\r\n% Using the imagery from the WMS server, calculate the hurricane's area and\r\n% location of the center by assuming that the hurricane is in the region\r\n% with the most cloud cover (as defined in pixel units, not map units).\r\n\r\n%%\r\n% The first step is to create a black and white image from the image of\r\n% Katrina. As an initial guess, try using the |graythresh| function to\r\n% calculate a threshold value.\r\nthresh = graythresh(katrinaMap);\r\nkatrina_bw = im2bw(katrinaMap, thresh);\r\nfigure('Color', 'white', 'Renderer', 'zbuffer')\r\nusamap(katrina.Latlim, katrina.Lonlim)\r\ngeoshow(katrina_bw, R)\r\ntitle('Black and White Cloud Image Using graythresh')\r\n\r\n%%\r\n% In the black and white image, you can see that there exists large tracts\r\n% that are now white even though these areas do not have significant cloud\r\n% cover. One method to remove these regions is to set a higher threshold\r\n% value.\r\nfactor = 1.9;\r\nthresh = factor*thresh;\r\nkatrina_bw = im2bw(katrinaMap, thresh);\r\nfigure('Color', 'white', 'Renderer', 'zbuffer')\r\nusamap(katrina.Latlim, katrina.Lonlim)\r\ngeoshow(katrina_bw, R)\r\ntitle('Black and White Cloud Image Using 1.9*graythresh')\r\n\r\n%%\r\n% Calculate the convex hulls and areas of all regions using |regionprops|.\r\nkatrina_props = regionprops(katrina_bw, {'ConvexHull', 'Area'});\r\n\r\n%%\r\n% Find the index that outlines the maximum area (in pixel units).\r\n[~, maxIndex] = max([katrina_props(:).Area]);\r\nmaxIndex = maxIndex(1);\r\n\r\n%%\r\n% Convert the row and column indices where the area is at a maximum to\r\n% latitude and longitude coordinates. Store the coordinates into a\r\n% structure.\r\nrow = katrina_props(maxIndex).ConvexHull(:,2);\r\ncol = katrina_props(maxIndex).ConvexHull(:,1);\r\nConvexHull = struct('Lat', [], 'Lon', [], 'Area', []);\r\n[ConvexHull.Lat, ConvexHull.Lon] = pix2latlon(R, row, col);\r\n\r\n%%\r\n% Plot the convex hull onto the black and white cloud image.\r\nfigure('Color', 'white', 'Renderer', 'zbuffer')\r\nusamap(katrina.Latlim, katrina.Lonlim)\r\ngeoshow(katrina_bw, R)\r\ngeoshow(ConvexHull.Lat, ConvexHull.Lon, 'Color',  'cyan',  'LineWidth',  2)\r\ntitle('Convex Hull Image')\r\n\r\n%%\r\n% Calculate the area of the convex hull in units of square statute miles.\r\nellipsoid = almanac('earth','wgs84','statutemiles');\r\nConvexHull.Area = areaint(ConvexHull.Lat, ConvexHull.Lon, ellipsoid);\r\n\r\n%%\r\n% Calculate the boundaries of all objects using |bwboundaries|.\r\n[B, L] = bwboundaries(katrina_bw);\r\n\r\n%%\r\n% Convert the label matrix to an RGB image and display it.\r\nlabelRGB = label2rgb(L, @jet, [.5 .5 .5]);\r\nfigure('Color', 'white', 'Renderer', 'zbuffer')\r\nusamap(katrina.Latlim, katrina.Lonlim)\r\ngeoshow(labelRGB, R)\r\ntitle('Boundary Image')\r\n\r\n%%\r\n% Use a histogram to select the region with the most cloud cover.\r\nclouds = L(L~=0);\r\n[freq, loc] = hist(clouds ,unique(L));\r\n[~, maxIndex] = max(freq);\r\nhurricane = loc(maxIndex);\r\n\r\n%%\r\n% Capture indices of the largest cloud cover.\r\nindices = B{hurricane};\r\nrow = indices(:,1);\r\ncol = indices(:,2);\r\n\r\n%%\r\n% Convert indices of the largest cloud cover to latitude and longitude.\r\nBoundary = struct('Lat', [], 'Lon', [], 'Area', []);\r\n[Boundary.Lat, Boundary.Lon] = pix2latlon(R, row, col);\r\n\r\n%%\r\n% Plot the boundary onto the labeled image.\r\nfigure('Color', 'white', 'Renderer', 'zbuffer')\r\nusamap(katrina.Latlim, katrina.Lonlim)\r\ngeoshow(labelRGB, R)\r\ngeoshow(Boundary.Lat, Boundary.Lon, 'Color', 'white', 'LineWidth', 2)\r\ntitle('Boundary Image with Outline')\r\n\r\n%%\r\n% Calculate the area of the boundary in units of square statute miles.\r\nBoundary.Area = areaint(Boundary.Lat, Boundary.Lon, ellipsoid);\r\n\r\n%%\r\n% Find the center of the hurricane.\r\nkatrina_center = regionprops(L, 'Centroid');\r\ncenter_row = katrina_center(hurricane).Centroid(:,2);\r\ncenter_col = katrina_center(hurricane).Centroid(:,1);\r\nCentroid = struct('Lat', [], 'Lon', []);\r\n[Centroid.Lat, Centroid.Lon] = pix2latlon(R, center_row, center_col);\r\n\r\n%% Step  7: Display Katrina Map and Vector Overlays\r\n\r\n%%\r\n% Import hurricane track vector data. The vector dataset, |katrina.shp|,\r\n% contains the historical track for hurricane Katrina. The track file\r\n% contains the 6-hourly (0000, 0600, 1200, 1800 UTC) center locations and\r\n% intensities of the hurricane.\r\ntracks = shaperead('KatrinaData\/katrina.shp', 'UseGeo', true);\r\ntracks(1)\r\n\r\n%%\r\n% Display the WMS basemap and vector overlays:\r\n% \r\n% * Overlay the outlines of worldwide land areas.\r\n% * Overlay the outline of the United States.\r\n% * Overlay the outline of the convex hull in red.\r\n% * Overlay the outline of the boundary of the hurricane in cyan.\r\n% * Overlay the center of the hurricane in black.\r\n% * Overlay the eye of the hurricane track in yellow.\r\n\r\nfigure('Color', 'white', 'Renderer', 'zbuffer')\r\nusamap(katrina.Latlim, katrina.Lonlim)\r\ngeoshow(katrinaMap, R)\r\ngeoshow(land, 'FaceColor', 'none')\r\ngeoshow(states, 'FaceColor', 'none')\r\ngeoshow(ConvexHull.Lat, ConvexHull.Lon, ...\r\n    'LineWidth', 4, 'Color', 'red');\r\ngeoshow(Boundary.Lat, Boundary.Lon, ...\r\n    'LineWidth', 2, 'Color', 'cyan');\r\ngeoshow(Centroid.Lat, Centroid.Lon, ...\r\n    'LineWidth', 4, 'Color', 'black', 'Marker', '+')\r\ngeoshow(tracks, 'LineWidth', 2, 'Color', 'yellow')\r\ntitle({katrina.LayerTitle, katrina.Details.Dimension.Default}, ...\r\n    'Interpreter', 'none', 'FontWeight', 'bold');\r\n    \r\n%%\r\n% Print the area statistics.\r\nfprintf('Area of the convex hull is %0.2f square statute miles.\\n', ...\r\n    ConvexHull.Area);\r\nfprintf('Area of the boundary    is %0.2f square statute miles.\\n', ...\r\n    Boundary.Area);\r\n\r\n%% Step  8: Apply Algorithm to all Frames Using |parfor|\r\n% Now that you have developed an algorithm to track and calculate the\r\n% hurricane's size and position, the next step is to apply that algorithm\r\n% to a time sequence of data from the server. Use the |parfor| command to\r\n% distribute the calculations to multiple processors. It is part of the\r\n% MATLAB(R) language, and behaves essentially like a regular |for|-loop if\r\n% you do not have access to the Parallel Computing Toolbox(TM) product.\r\n\r\n%%\r\n% Calculate the total number of animation frames.\r\nnumFrames = endDay - startDay + 1;\r\n\r\n%%\r\n% Assign the initial requestTime.\r\nrequestTime = startTime;\r\n\r\n%% \r\n% Initialize variables to hold the data calculated from each frame.\r\ntime = cell(numFrames, 1);\r\nS = struct('Lat', [], 'Lon', [], 'Area', []);\r\nConvexHull = S;\r\nBoundary = S;\r\nConvexHull(numFrames) = ConvexHull;\r\nBoundary(numFrames) = Boundary;\r\nCentroid = struct('Lat', [], 'Lon', []);\r\nCentroid(numFrames) = Centroid;\r\nframeName = datafile('frame_');\r\n\r\n%% \r\n% Check the status of the MATLAB(R) pool and determine if it is open.\r\n% Use the MATLAB(R) pool to allow the body of the |parfor| loop to\r\n% run in parallel. \r\npoolSize = matlabpool('size');\r\nif poolSize == 0\r\n    fprintf('The MATLAB pool is not open, the loop will run sequentially.');\r\nelse\r\n    fprintf('Using %d workers.\\n', poolSize);\r\nend\r\n\r\n%%\r\n% For each day, from |startDay| to |endDay|, perform the following\r\n% operations:\r\n%\r\n% * obtain the raster map for the time frame from either the web map server \r\n%   or from the cache file,\r\n% * calculate the convex hull, boundary, and their areas, \r\n% * calculate the location of the center of the hurricane,\r\n% * render the data into a Figure window, and \r\n% * save the data from the graphics handles to a MAT-file.\r\n\r\n%% \r\n% To create the animation frames, you can load the rendered data from each\r\n% saved MAT-file and capture the data to an image. You need to do this\r\n% operation outside of the |parfor| loop since each compute lab is started\r\n% in a headless mode. By storing each rendered frame into a MAT-file within\r\n% the |parfor| loop, you can improve performance since rendering the data\r\n% into a Figure window from a MAT-file is a faster operation than the\r\n% display operations for the vector and raster data.\r\ntworkers = zeros(1, poolSize);\r\ntic;\r\nparfor k=1:numFrames\r\n    tstart = tic;\r\n    % Update currentDay.\r\n    currentDay = startDay + k -1;\r\n    \r\n    % Assign time value for this frame.\r\n    frameTime = requestTime;\r\n    frameTime(dayIndex) = num2str(currentDay);\r\n    time{k} = frameTime;\r\n\r\n    % Get the WMS map from the server (or file) for this time period.\r\n    cachefile = datafile(['katrina_' num2str(currentDay) '.tif']);\r\n    if useInternet\r\n         RGB = wmsread(katrina, 'Time', frameTime);\r\n         if ~exist(cachefile, 'file')\r\n            imwrite(RGB, cachefile);\r\n         end\r\n    else\r\n        RGB = imread(cachefile);\r\n    end\r\n    \r\n    % Calculate the convex hull, boundary, areas, and centroid using the\r\n    % same threshold value as calculated for the first image.\r\n    [ConvexHull(k), Boundary(k), Centroid(k)] ...\r\n        = calculateRegionProperties(RGB, R, thresh);\r\n    \r\n    % Render the frame.\r\n    hFig = figure('Color', 'white', 'Renderer', 'zbuffer');\r\n    renderFrame(frameTime, katrina, RGB,  R,  ...\r\n        ConvexHull(k), Boundary(k), Centroid(k), land, states, tracks);\r\n    \r\n    % Save the rendered frame.\r\n    hAxes = get(hFig, 'CurrentAxes');\r\n    frameFilename = [frameName num2str(k) '.fig'];\r\n    hgsave(hAxes, frameFilename);\r\n    close(hFig);\r\n    \r\n     % Save the workers elapsed time.\r\n    tworkers(k) = toc(tstart);\r\n     \r\n    % Print statistics\r\n    fprintf('\\nFrame %d elapsed time: %4.8f\\n', k, tworkers(k));\r\n    fprintf('Frame %d convex hull area: %0.2f.\\n', ...\r\n        k, ConvexHull(k).Area);\r\n    fprintf('Frame %d boundary area:    %0.2f.\\n', ...\r\n        k, Boundary(k).Area);\r\nend\r\ntloop = toc;\r\n\r\n%%\r\n% You expect the total time of all the workers to be more than the loop\r\n% time when the loop is executed in parallel.\r\nfprintf('Total time for all workers: %4.8f\\n', sum(tworkers));\r\nfprintf('Total loop time:            %4.8f\\n', tloop);\r\n\r\n%% Step  9: Create Animation Files\r\n% An animation can be viewed in the browser when the browser opens an\r\n% animated GIF file or an AVI video file. To create the animation frames of\r\n% the WMS basemap and vector overlays, render the data from the saved\r\n% graphics files, capture the rendered data into a |frame| structure, and\r\n% save |frame| into either an array for the animated GIF file or to an\r\n% AVI video file.\r\n\r\n%%\r\n% To share with others or to post to web video services, you can create an\r\n% AVI video file containing all the frames using the |VideoWriter| class.\r\nbasename = 'wmstrack';\r\nvideoFilename = ['html\/' basename '.avi'];\r\nif exist(videoFilename, 'file')\r\n    delete(videoFilename)\r\nend\r\nwriter = VideoWriter(videoFilename);\r\nwriter.FrameRate = 1;\r\nwriter.Quality = 100;\r\nwriter.open;\r\n\r\n%%\r\n% For each frame, read the graphics from the cache file created in the body\r\n% of the |parfor| loop, and display the WMS basemap and vector overlays.\r\n% Capture each frame and store into the frames array. The loop cannot run\r\n% in parallel since the video frames must be ordered.\r\nhFig = figure('Color', 'white', 'Renderer', 'zbuffer');\r\nfor k=1:numFrames\r\n \r\n    % Delete the previous frame.\r\n    hAxes = get(hFig, 'CurrentAxes');\r\n    delete(hAxes)\r\n    \r\n    % Load the basemap and vector overlays into the axes.\r\n    frameFilename = [frameName num2str(k) '.fig'];\r\n    hgload(frameFilename);\r\n    \r\n    % Capture the current frame.\r\n    shg\r\n    frame = getframe(hFig);\r\n    \r\n    % Store the frame into the |animated| array for the GIF file.\r\n    if k == 1\r\n        [animated, cmap] = rgb2ind(frame.cdata, 256, 'nodither');\r\n        animated(:,:,1,numFrames) = animated;\r\n    else\r\n        animated(:,:,1,k) = rgb2ind(frame.cdata, cmap, 'nodither');\r\n    end\r\n    \r\n    % Write the frame to the AVI file.\r\n    writer.writeVideo(frame);\r\nend\r\n% Close the Figure window and AVI file.\r\nclose\r\nwriter.close;\r\n\r\n%%\r\n% Write the animated GIF file.\r\nfilename = [basename '.gif'];\r\ndelayTime = 2.0;\r\nimwrite(animated, cmap, filename, 'DelayTime', delayTime, 'LoopCount', inf);\r\n\r\n%% Step 10: View Animated GIF File\r\n% An animation can be viewed in the browser when the browser opens an\r\n% animated GIF file. \r\nweb(filename)\r\n\r\n%%\r\n% <<wmstrack.gif>>\r\n\r\n%% Step 11: Display Table of Time, Location, and Area\r\nheaders = {'Time',  'Latitude',  'Longitude',  ...\r\n    'Convex Hull Area', 'Boundary Area'};\r\nmakeTable(headers, time, [Centroid.Lat], [Centroid.Lon],  ...\r\n    [ConvexHull.Area], [Boundary.Area]);\r\n\r\n%% Step 12: View Video File\r\n% The following <wmstrack.avi wmstrack.avi> video is generated\r\n% and can be viewed using an external viewer.\r\n\r\n%%\r\n% Some browsers support viewing the video in an embedded frame:\r\nwebvideo(videoFilename);\r\n\r\n%% The makeTable Function\r\ndbtype makeTable\r\n\r\n%% The webvideo Function\r\ndbtype webvideo\r\n\r\nend\r\n\r\n%% The calculateRegionProperties Function\r\n% This function calculates the convex hull, BW boundary, area, and centroid\r\n% properties of an RGB image.\r\nfunction [ConvexHull, Boundary, Centroid] = ...\r\n    calculateRegionProperties(RGB, R, thresh)\r\n\r\n% Calculate the convex hull and area of the objects in the RGB image.\r\nbw = im2bw(RGB, thresh);\r\nprops = regionprops(bw, {'ConvexHull', 'Area'});\r\n\r\n% Find the indices that outline the maximum area.\r\n[~, maxIndex] = max([props(:).Area]);\r\n\r\n% Convert the row and column indices where the area is at a maximum to\r\n% latitude and longitude coordinates. Store the coordinates into a\r\n% structure.\r\nrow = props(maxIndex).ConvexHull(:,2);\r\ncol = props(maxIndex).ConvexHull(:,1);\r\nConvexHull = struct('Lat', [], 'Lon', [], 'Area', []);\r\n[ConvexHull.Lat, ConvexHull.Lon] = pix2latlon(R, row, col);\r\n\r\n% Calculate the area of the convex hull in units of square statute miles.\r\nellipsoid = almanac('earth','wgs84','statutemiles');\r\nConvexHull.Area = areaint(ConvexHull.Lat, ConvexHull.Lon, ellipsoid);\r\n\r\n% Calculate the boundaries of all objects.\r\n[B, L] = bwboundaries(bw);\r\n\r\n% Use a histogram to select the largest region.\r\nobjects = L(L~=0);\r\n[freq, loc] = hist(objects ,unique(L));\r\n[~, maxIndex] = max(freq);\r\nlargestRegion = loc(maxIndex);\r\n\r\n% Capture indices of the largest region.\r\nindices = B{largestRegion};\r\nrow = indices(:,1);\r\ncol = indices(:,2);\r\n\r\n% Convert indices of the largest region cover to latitude and longitude.\r\nBoundary = struct('Lat', [], 'Lon', [], 'Area', []);\r\n[Boundary.Lat, Boundary.Lon] = pix2latlon(R, row, col);\r\n\r\n% Calculate the area of the boundary in units of square statute miles.\r\nBoundary.Area = areaint(Boundary.Lat, Boundary.Lon, ellipsoid);\r\n\r\n% Find the center of the largest region.\r\ncenter = regionprops(L, 'Centroid');\r\ncenter_row = center(largestRegion).Centroid(:,2);\r\ncenter_col = center(largestRegion).Centroid(:,1);\r\nCentroid = struct('Lat', [], 'Lon', []);\r\n[Centroid.Lat, Centroid.Lon] = pix2latlon(R, center_row, center_col);\r\nend\r\n\r\n%% The renderFrame Function\r\nfunction renderFrame(frameTime, layer, layerMap, R,  ...\r\n    ConvexHull, Boundary, Centroid, land, states, tracks)\r\n% This function renders the frame data into a Figure window.\r\n\r\nusamap(layer.Latlim, layer.Lonlim)\r\ngeoshow(layerMap, R)\r\ngeoshow(land, 'FaceColor', 'none')\r\ngeoshow(states, 'FaceColor', 'none')\r\ngeoshow(ConvexHull.Lat, ConvexHull.Lon, 'LineWidth', 4, 'Color', 'red');\r\ngeoshow(Boundary.Lat, Boundary.Lon, 'LineWidth', 2, 'Color', 'cyan');\r\ngeoshow(Centroid.Lat, Centroid.Lon,   'LineWidth',  4,  'Color', 'k',  'Marker', '+')\r\ngeoshow(tracks, 'LineWidth', 2, 'Color',  'yellow')\r\ntitle({layer.LayerTitle, frameTime}, 'Interpreter', 'none', 'FontWeight', 'bold');\r\n\r\n%% Credits\r\n%\r\n% _Katrina Layer_ \r\n% \r\n% The Katrina layer used in the demo is from the NASA Goddard Space Flight\r\n% Center's SVS Image Server and is maintained by the Scientific\r\n% Visualization Studio.\r\n%\r\n% For more information about this server, run: \r\n%    \r\n%    >> wmsinfo('http:\/\/aes.gsfc.nasa.gov\/cgi-bin\/wms?')\r\n%\r\n% _usastatehi.shp_\r\n%\r\n% This vector dataset contains moderate-resolution outlines and names for\r\n% the 50 United States plus the District of Columbia, stored in shapefile\r\n% format with coordinates in latitude (Y) and longitude (X) in degrees.\r\n% It is based on U.S. Census Bureau's TIGER thinned boundary files. The\r\n% data file is shipped with the Mapping Toolbox\u00e2\u201e\u00a2.\r\n%\r\n% _landareas.shp_\r\n%\r\n% This vector dataset contains low-resolution outlines of land areas\r\n% worldwide, stored as polygons in shapefile format with coordinates in\r\n% latitude (Y) and longitude (X) in degrees. The raw coordinate data were\r\n% derived from the Digital Chart of the World (DCW) browser layer,\r\n% published by U.S. National Geospatial-Intelligence Agency (NGA), formerly\r\n% National Imagery and Mapping Agency (NIMA). The data file is shipped with\r\n% the Mapping Toolbox\u00e2\u201e\u00a2.\r\n%\r\n% _katrina.shp_\r\n%\r\n% This vector dataset contains the historical track for hurricane Katrina.\r\n% The raw data were obtained from the atl_hurtrack.shp file from the\r\n% National Oceanic and Atmospheric Administration Coastal Services Center.\r\n% This Historical North Atlantic Tropical Cyclone Tracks file contains the\r\n% 6-hourly (0000, 0600, 1200, 1800 UTC) center locations and intensities\r\n% for all subtropical depressions and storms, extra-tropical storms,\r\n% tropical lows, waves, disturbances, depressions and storms, and all\r\n% hurricanes, from 1851 through 2009. These data are intended for\r\n% geographic display and analysis at the national level, and for large\r\n% regional areas. The data should be displayed and analyzed at scales\r\n% appropriate for 1:2,000,000-scale data. No responsibility is assumed by\r\n% the National Oceanic and Atmospheric Administration in the use of these\r\n% data. The complete track data was downloaded from NOAA. The data file is\r\n% not shipped with the Mapping Toolbox\u00e2\u201e\u00a2.\r\n\r\n%% Can You Apply Parallel Processing to Your Processing in an Interesting Way?\r\n% As you can see here, I've used |parfor| to independently fetch and\r\n% process each frame before constructing the final video.  Is that\r\n% technique something you might be able to take advantage of?  Let us know\r\n% <https:\/\/blogs.mathworks.com\/loren\/?p=262 here>.\r\nend\r\n\r\n##### SOURCE END ##### 5cb939f7486747759f02d0a247201b95\r\n-->","protected":false},"excerpt":{"rendered":"<p>\r\n   \r\n      Guest blogger, Kelly Luetkemeyer, who is a Mapping Toolbox software developer at The MathWorks, returns with another article about using and analyzing data\r\n         obtained from a Web... <a class=\"read-more\" href=\"https:\/\/blogs.mathworks.com\/loren\/2011\/01\/20\/tracking-a-hurricane-using-web-map-service-wms\/\">read more >><\/a><\/p>","protected":false},"author":39,"featured_media":0,"comment_status":"open","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":[],"categories":[27,41,34],"tags":[],"_links":{"self":[{"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/posts\/262"}],"collection":[{"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/users\/39"}],"replies":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/comments?post=262"}],"version-history":[{"count":0,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/posts\/262\/revisions"}],"wp:attachment":[{"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/media?parent=262"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/categories?post=262"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/tags?post=262"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}