{"id":257,"date":"2009-04-24T15:29:17","date_gmt":"2009-04-24T19:29:17","guid":{"rendered":"https:\/\/blogs.mathworks.com\/steve\/2009\/04\/24\/continental-divide-1-intro\/"},"modified":"2019-10-28T15:30:25","modified_gmt":"2019-10-28T19:30:25","slug":"continental-divide-1-intro","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/steve\/2009\/04\/24\/continental-divide-1-intro\/","title":{"rendered":"Locating the US continental divide, part 1 – Introduction"},"content":{"rendered":"
\r\n

Today I'm introducing a series on computing the location of the US continental divide. (I was inspired to tackle this topic\r\n by Brian Hayes' blog post<\/a> on the divide; thanks very much to Mike<\/a> and Ned<\/a> for sending me the link.)\r\n <\/p>\r\n

The series will start with importing global elevation data for Earth, and it will finish with this visualization of the Atlantic\r\n and Pacific catchment basins for the continental US:\r\n <\/p>\r\n

<\/p>\r\n

My regular blog readers might be happy to know that this time I've planned and written out the whole series in advance. See\r\n below for the series summary.\r\n <\/p>\r\n

In this first post, I'll tackle data import and display.<\/p>\r\n

The data for this project comes from the Global Land One-km Base Elevation Project<\/a>, available from the National Geophysical Data Center. The global elevation data set is available online<\/a> in the form of 16 tiles, labeled A through P. Tile E<\/a> and and tile F<\/a> cover the continental United States, Central America, and portions of Canada and South America.\r\n <\/p>\r\n

After you've downloaded and uncompressed the files, you have to figure out how to import them. As it turns out, these files are\r\n not in a structured format, such as TIFF or HDF, that can be read automatically. Instead, you need to poke around the data\r\n documentation on the web to find out how to read the files. This is fairly common with scientific data sources. The documentation\r\n for this data set includes a section called \"Data Format & Importing GLOBE Data.\"<\/a> By perusing this section I found<\/a> that the files I downloaded contain 16-bit signed integers organized in a raster of 6000 rows and 10800 columns, stored by\r\n row. There are no header bytes, and the integers are stored using little-endian byte order.\r\n <\/p>\r\n

The appropriate MATLAB function to use is multibandread<\/tt>.\r\n <\/p>

data_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\n\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\n\r\ndem = [E, F];  % dem stands for \"Digital Elevation Model\"<\/span><\/pre>
Next I'll use imshow<\/tt>, using some optional inputs, to display the digital elevation model as an image.\r\n   <\/p>
imshow(dem, [], 'InitialMagnification'<\/span>, 'fit'<\/span>)<\/pre> 
About those optional inputs: When displaying gridded data as an image, it is often useful to control the black-to-white display\r\n      range.  If you supply a vector as the second input to imshow<\/tt>, like this:\r\n   <\/p>
   imshow(dem, [a b])<\/pre>
then imshow<\/tt> uses [a b]<\/tt> as the image display range. A value of a<\/tt> gets displayed as black, a value of b<\/tt> gets displayed as white, and values in-between are displayed as shades of gray.\r\n   <\/p>\r\n   
For the special case where the second input is empty:<\/p>
   imshow(dem, [])<\/pre>
imshow<\/tt> uses the data range as the display range. That is, imshow(dem, [])<\/tt> is equivalent to imshow(dem, [min(dem(:)), max(dem(:))])<\/tt>.\r\n   <\/p>\r\n   
The other optional inputs, 'InitialMagnification', 'fit'<\/tt>, tell imshow<\/tt> to magnify or shrink the image display as needed to fit nicely within the current figure window. I'm using that option here\r\n      so that my image displays will fit in the column for this blog post.\r\n   <\/p>\r\n   
For this data, the autoscaling provided by the imshow(dem, [])<\/tt> syntax has given us an image display that is too dark, particularly on the eastern part of the continent.  Let's set a white\r\n      level that's lower than the data maximum; that will have the effect of brightening the image.\r\n   <\/p>
min(dem(:))<\/pre>
\r\nans =\r\n\r\n   -500\r\n\r\n<\/pre>
max(dem(:))<\/pre>
\r\nans =\r\n\r\n   6085\r\n\r\n<\/pre>
display_range = [-500 3000];\r\nimshow(dem, display_range, 'InitialMagnification'<\/span>, 'fit'<\/span>)<\/pre> 
Next I'll crop the image. The Image Processing Toolbox has interactive cropping tools (imcrop<\/tt>, imtool<\/tt>), but I can't use those here.\r\n   <\/p>
dem_cropped = dem(1:4000, 6000:14500);\r\nimshow(dem_cropped, display_range, 'InitialMagnification'<\/span>, 'fit'<\/span>)<\/pre> 
Finally, I'll save dem_cropped<\/tt> to a MAT-file for use in the rest of the series.\r\n   <\/p>
save continental_divide<\/span> dem_cropped<\/span><\/pre>
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