{"id":263,"date":"2009-06-05T07:00:53","date_gmt":"2009-06-05T11:00:53","guid":{"rendered":"https:\/\/blogs.mathworks.com\/steve\/2009\/06\/05\/continental-divide-7-all-steps\/"},"modified":"2019-10-28T15:36:25","modified_gmt":"2019-10-28T19:36:25","slug":"continental-divide-7-all-steps","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/steve\/2009\/06\/05\/continental-divide-7-all-steps\/","title":{"rendered":"Locating the US continental divide, part 7 – Putting it all together"},"content":{"rendered":"
\r\n

Today is the final post in my continental divide series<\/a>. In earlier posts I have used the problem of computing the US continental divide as a vehicle for exploring data import,\r\n image display scaling, the watershed transform, label matrices, local and regional minima, binary image manipulation, boundary\r\n tracing, and visualization techniques.\r\n <\/p>\r\n

For this last post, I thought it would useful to gather together all the computational steps in one short script, starting\r\n with data import and finishing with the visualization.\r\n <\/p>

% Import the tiles e10g and f10g.<\/span>\r\ndata_size = [6000 10800 1];  % The data has 1 band.<\/span>\r\nprecision = 'int16=>int16'<\/span>;  % Read 16-bit signed ints into a int16 array.<\/span>\r\nheader_bytes = 0;\r\ninterleave = 'bsq'<\/span>;          % Band sequential. Not critical for 1 band.<\/span>\r\nbyte_order = 'ieee-le'<\/span>;\r\nE = multibandread('e10g'<\/span>, data_size, precision, header_bytes, ...<\/span>\r\n    interleave, byte_order);\r\nF = multibandread('f10g'<\/span>, data_size, precision, header_bytes, ...<\/span>\r\n    interleave, byte_order);\r\nEF = [E, F];\r\n\r\n% Crop the data.<\/span>\r\ndem = EF(1:4000, 6000:14500);\r\n\r\n% Form an ocean mask.<\/span>\r\nocean_mask = dem == -500;\r\n\r\n% Modify the ocean mask so it contains only the two connected components on<\/span>\r\n% the left and right side.<\/span>\r\n[M, N] = size(dem);\r\nocean_mask = bwselect(ocean_mask, [1 N], [1 1]);\r\n\r\n% Modify the DEM so that its only regional minima are the ocean_mask<\/span>\r\n% pixels.<\/span>\r\ndem_modified = imimposemin(dem, ocean_mask);\r\n\r\n% Compute the watershed transform of the modified DEM.<\/span>\r\nL = watershed(dem_modified);\r\n\r\n% Visualize the result. First, convert the output of watershed to a color<\/span>\r\n% image.<\/span>\r\npacific = [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\n\r\n% Superimpose the colored watershed image over the DEM.<\/span>\r\nimshow(dem, [-500 3000], 'InitialMagnification'<\/span>, 'fit'<\/span>)\r\nhold on<\/span>\r\nh = imshow(rgb);\r\nset(h, 'AlphaData'<\/span>, 0.2);\r\n\r\n% Trace the watershed ridge line and superimpose it as a thick red line.<\/span>\r\nb = 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> 
That's it! I hope you enjoyed this series.<\/p>\r\n   
Do you have your own ideas for fun computation and visualization using DEM data?  Please comment.<\/p>\r\n   
About this Series<\/i><\/p>\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   
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