{"id":262,"date":"2009-05-29T07:00:09","date_gmt":"2009-05-29T11:00:09","guid":{"rendered":"https:\/\/blogs.mathworks.com\/steve\/2009\/05\/29\/continental-divide-6-watershed\/"},"modified":"2019-10-28T15:35:12","modified_gmt":"2019-10-28T19:35:12","slug":"continental-divide-6-watershed","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/steve\/2009\/05\/29\/continental-divide-6-watershed\/","title":{"rendered":"Locating the US continental divide, part 6 – Final computation and visualization"},"content":{"rendered":"
\r\n

Today's post in my continental divide series<\/a> demonstrates the final computational step and shows one procedure for visualizing the Pacific and Atlantic catchment basins.\r\n <\/p>\r\n

In my previous post, I modified the DEM (digital elevation model) data set so that its only regional minima were the two oceans.<\/p>

s = load('continental_divide'<\/span>);\r\ndem_modified = s.dem_modified;\r\nimshow(dem_modified, [-500 3000], 'InitialMagnification'<\/span>, 'fit'<\/span>);<\/pre> 
Since the modified DEM has only two regional minima, we expect to get only two catchment basins when we compute the watershed\r\n      transform.\r\n   <\/p>
L = watershed(dem_modified);\r\nmax(L(:))<\/pre>
\r\nans =\r\n\r\n     2\r\n\r\n<\/pre>
I'll use label2rgb<\/tt> to display the Pacific catchment basin in green, the Atlantic catchment basin in blue, and the watershed ridge pixels (the\r\n      continental divide) in red.\r\n   <\/p>
pacific = [0.2 1 0.2];\r\natlantic = [0.3 0.3 1];\r\nridge = [1 0 0];\r\nrgb = label2rgb(L, [pacific; atlantic], ridge);\r\nimshow(rgb, 'InitialMagnification'<\/span>, 'fit'<\/span>)<\/pre> 
In a February blog post<\/a>, I showed how to overlay one image transparently over another.  We can use that technique here to achieve an effect similar\r\n      to the visualization that Brian Hayes used in his blog post on the continental divide<\/a>, where the ocean catchment basins are shown shaded with two different colors.\r\n   <\/p>
imshow(dem, [-500 3000], 'InitialMagnification'<\/span>, 'fit'<\/span>)\r\nhold on<\/span>\r\nh = imshow(rgb);\r\nset(h, 'AlphaData'<\/span>, 0.2);<\/pre> 
It would be nice if we could make the watershed ridge line itself more visible.  For example, if we could use plot<\/tt>, then we could use a thick line. But right now we just have the ridge line represented as zero-valued elements of a raster,\r\n      instead of the x and y vectors we would need to pass to plot<\/tt>.  We could use find<\/tt> to get the coordinates of the ridge pixels, but we wouldn't get them in the right order.  One solution is to use bwboundaries<\/tt> to \"trace\" the zero-valued ridge pixels.  Note that I specify 4-connectivity in the call to bwboundaries<\/tt> below. I've found that I get better results using 4-connectivity when tracing a single-pixel-wide stroke that's not a closed\r\n      loop.\r\n   <\/p>
b = bwboundaries(L == 0, 4);\r\nb1 = b{1};  % There's only 1 boundary found here - the ridge line.<\/span>\r\nx = b1(:,2);\r\ny = b1(:,1);\r\nplot(x, y, 'r'<\/span>, 'LineWidth'<\/span>, 2);<\/pre> 
I have one more post planned for this series.  Next time I'll gather together all the computational steps in a single, short\r\n      script.\r\n   <\/p>\r\n   
About this Series<\/i><\/p>\r\n   
As an algorithm development problem, computing the divide lets us explore aspects of data import and visualization, manipulating\r\n      binary image masks, label matrices, regional minima, and the watershed transform.\r\n   <\/p>\r\n   
\r\n      
This series<\/a> explores the problem of computing the location of the continental divide for the United States.  The divide separates the\r\n      Atlantic and Pacific Ocean catchment basins for the North American continent.\r\n   <\/p>\r\n