Digital Elevation Model
|
URI: |
http://herbert.gandraxa.com/herbert/dem.asp |
|
Link template: |
<a href="http://herbert.gandraxa.com/herbert/dem.asp">Digital Elevation Model</a> |
|
Link symbols: |
|
On current page |
ZIP archive |
PDF |
E-Mail
Table of Contents
Article
Organization
Home »
Articles »
Digital Terrain Analysis »
Digital Elevation Model
Scope
This article is the first one in a series dealing with Digital Terrain Analysis. It explains the process of obtaining elevation data and verifying what we obtained. The next article provides the theoretical background to render the obtained data with isometrical projections.
Author
Published
2007-Feb-24
External Links
Digital Elevation Models on Wikipedia- NASA Jet Propulsion Laboratory, Homepage of the
Shuttle Radar Topography Mission - National Geospatial-Intelligence Agency (NGA), Homepage
- U.S. Geological Survey (USGS), Homepage
Basics
Terminology
Digital Elevation Models are collections of data to represent the different elevations of surface area. In their most intuitive form they come as a two-dimensional matrix, where the rows represent latitudes and the columns longitudes. The individual cells are then filled with the elevation for the given latitude/longitude combination. The elevation most often is given in meters above (or below) sea level.
Data Sources
There are many sources on the Internet which provide elevation data for many regions. However, many of them either lack accuracy, or don't provide the desired resolution, or are restricted to very small areas, or all of these.
Fortunately, since the advent of the Shuttle Radar Topography Mission (SRTM) things have changed. The aim of that project was set high: it encomprised measuring the elevation of a vast part of the Earth (more accurately the area between 56° Southern latitude and 60° Northern latitude), at a resolution of one arc second (corresponding to approx. 30 meters horizontally and vertically). The project was carried out by the Space Shuttle Endeavour during its mission in February of 2000, in which all measurements occured with the help of two radar antennas placed 60 meters apart, employing a technique called interferometric Synthetic Aperture Radar.
For us non-US-government people, it is absolutely a highlight that the National Geospatial-Intelligence Agency (NGA) and the National Aeronautics and Space Administration (NASA), the leaders of this international project, decided that this material shall be available for free, for everyone having access to the Internet [1].
It comes ever better, though: much of the US territory data was improved to an accuracy of 1/3 arc seconds, meaning that the resolution was brought to approx. 10 meters between two measurement points. The data was made available in a very comfortable manner by the U.S. Geological Survey (USGS). The section after the next one is dealing with their web interface.
[1] Note, that the availability of the data is somewhat restricted: the resolution of one arc second is available for US territories only; for all other areas the data made available to the public has a resolution of 3 arc seconds, corresponding to approx. 90 meters between two measurement points.
Quantities
Data Volume
Before retrieval of actual data, let's consider for a moment about what data volumes we are talking.
Let's assume, that we have available data for the whole planet, at a resolution of 10 meters. The Earth has a surface of approx. 510 million km2, of which approx. 149 million km2 (29.2%) is land. At an accuracy of 10 meters, the data would amount to 100 x 100 = 10,000 elevation measurements per square kilometer, and thus (if we restrict ourselves to the land areas) to a total of 1.49 trillion data points (short scale terminology, i.e. 1 trillion = 1 million millions).
Depending on the desired resolution, storing a single data point requires more than just a byte, though: a byte at 8 bits can only hold 28=256 distinct values. If we aim for an elevation expressed in meters above sea level, then we need at least 14 bits (for values between 0 and 16,383 meters) to accommodate such heights as the Mount Everest's peak at 8,848 meters. And if we also want to express negative values, another bit is needed to act as a sign bit. Since 15 bits are quite unhandy to deal with in modern computer systems, two bytes with a total of 16 bits will be used instead (and this is in fact one of the data formats which is provided by the USGS).
From above it became obvious, that the total storage requirements amounts to approx. 3 Terabytes, or at the time of writing this article about thrice the size of the largest hard disk drives available for PCs [1].
[1] Hitachi Deskstar 7K1000.
Evaluating a Download
For the reasons given above, we are forcing ourselves to some modesty: what we are looking for is high-resolution data (1/3 arc second, or 10 meters) not exceeding hard disk space of about 100 MB. Given a comparably modern PC, this quantity has the following advantages:
- the download via a DSL or cable connection should be feasible (a few minutes at max.);
- it should be possible to store it onto your hard disk drive; and
- it should be possible to process it within your local RAM.
100 MB imply, that we are looking for approx. 50 million data points (because each elevation point requires 16 bits = 2 bytes to be encoded). These 50 million points cover a rectangular area of whatever dimensions. To get an estimation of what we are looking for, let's assume the area is approximately square. The square root of 50 million will tell us the length of the square's edges then (in data points): the result is approx. 7,000 data points, and since we are after a resolution of 10 meters, we can go for a square with an edge length of approx. 70,000 meters = 70 km, covering an area of approx. 5,000 km2.
Because a goal of this article is the verification of once downloaded data, it would be ideal if we could tell quite instantly if we got it right. For this reason, it would be an advantage if we could find an island of roughly the area outlined above: assuming that the sea level is given as 0 (meters above sea level), we would be able to immediately identify the isolated landmass within a body of water, simply by distinguishing 0 and non-0 values and coloring the two entities differently.
Bearing in mind, that such high-resolution data was made available only for US territories, prospective area candidates were searched within the islands forming the US state of Hawaii. The "big island" Hawai'i is too big for our purpose (its landmass alone is larger than 10,000 km2), but the second-largest island Maui, having a landmass of less than 2,000 km2, with the inclusion of it's small Southwestern neighbor Kaho'olawe (approx. 116 km2) would suit us well: about half of the area consists of water surrounding the islands and the other half of land.
A quick check on Google Maps (Fig. 1) reveals, that the rectangle surrounding the two islands would have a dimension of roughly 80 km horizontally and 60 km vertically: exactly what we want.
Fig. 1: Finding a suitable area, in this case Maui and Kaho'olawe. © Google, Imagery © TerraMetrics, NASA, Map data © NAVTEQ ™
Obtaining Elevation Data
USGS Interface
Now that we know what we want, it's about the time to deal with the USGS web interface. It can be reached via http://seamless.usgs.gov/website/seamless/viewer.php. Notice that it might take a while until the map within the page has been loaded.
Fig. 2: Starting the USGS web interface
In the elevation menu to the right deselect the default mark and select "NED shaded relief (1/3 arc second)" instead (Fig. 3). You will notice that the map's appearance changes on deselection and selection to reflect the availability within the chosen resolution.
Fig. 3: Selecting the resolution
Since we by default already are in the "Zoom" menu, we can zoom into the appropriate area simply by dragging the mouse while its left button is pressed: a red rectangle will give you feedback about the currently selected area (Fig. 4).
Fig. 4: Zooming into the desired area
Upon releasing the mouse button, a zoomed-in map as per your rectangle will be loaded. Notice that the areas for which data is available appear highlighted. Repeat this process of zooming in closer and closer until you have selected the desired area, in the case of our example the islands of Maui and Kaho'olawe. If the resolution is high enough, the names of places are blended in to facilitate orientation.
Now select the "Downloads" menu group by clicking on the icon "Define Download Area". After this, the mouse essentially will react as before, i.e. one can click and drag to define a rectangular area within the currently displayed map. Being in the "Downloads" menu gives you an additional feedback (as explained in the notes below the map): if the color of the rectangle's borders is red, then the resulting file would be too large to be made available for download. This will be the case for data volumes exceeding 250 MB. You then either need to proceed in smaller steps and assemble the download files programmatically, or you need to order the the file on media if you indeed want the data as one single file. As long as the color of the rectangle is green, though, all is fine and you can continue. Try to not to be too generous with the surrounding water area: every pixel you include in your selection results in a larger download file (Fig. 5).
Fig. 5: Selecting the download area
Now another window opens with the details of your selection. Do not just accept the data presented there. Instead, use the option "Modify Data Request". You will notice, that all of the products listed here should be deselected, with the exception of "National Elevation Dataset (NED) 1 Arc Second" near the end of the list. This is not what we want: deselect it and choose "National Elevation Dataset (NED) 1/3 Arc Second" instead (Fig. 6). Notice, that this data is not the original SRTM project data: the 1 arc second and 3 arc seconds data are in the two bottommost options.
Just after your selection three scrolldown lists appear. If you follow this example, you must select "BIL" as the data format to be delivered, not "ArcGrid" as is suggested by default. "BIL" will return 16 bit data per elevation point, which is the easiest format to deal with. Also, select "TXT" as the data format, we don't want fancy markups, just the data. Depending on your computing environment, select the compression format ("ZIP" or "TGZ"). When done, click on the button labeled "Save changes and return to summary".
Fig. 6: Settings of your request
Once back to the summary, make sure you do indeed have the space on your machine to hold the requested data. Recall, that you will be delivered a compressed file, which you eventually must decompress to make use of it.
You schedule your download by clicking on the according button. Notice, that depending on the workload of the USGS server it might take a while until your requested data is ready, and also the download itself will need it's time. Time to get yourself a coffee while this happens.
Examing your Download
Once the file is on your machine, decompress the file you received into a folder of your choice. The file with the extension ".bil" contains your raw data. That file is quite large (in my case it has 96,001,532 bytes).
If you open it in a text editor capable of dealing with hex files (e.g. UltraEdit), you will notice a lot of 00 bytes at the start (at least if you were following this example). More accurately we should not speak about single bytes here, though, because the BIL format delivers 16 bits = 2 bytes per data point. So, you should notice a lot of 00 00 sequences. This is because the left upper corner of the selection rectangle was located in the ocean, and since this is at sea level, its elevation is 0 meters.
Scroll somewhat down and you will notice different data (Fig. 7). I marked two elevation points in it, one is the hexadecimal sequence 2B 07, and the other 2A 07. These sequences are in Little Endian order, which means, that the bytes must be interpreted from right to left, i.e. 07 2B and 07 2A. The hexadecimal values of 072B and 072A stand for decimal 1835 and 1834, resp., which is the elevation in meters above sea level. These two particular values are not unreasonable, as the Wikipedia article about the island Maui tells us, that the island's highest point is 3,055 meters (the Haleakala volcano).
Fig. 7: Examining the BIL file
You might have observed, that there is no delimiter which structures the obtained data. So how do we know how the data is structured then? As per explanations on the USGS site, the raw data is delivered row after row, starting at the Northwestern corner of your selection rectangle. Within a row there are several cells occupying 2 bytes each. But how do we know, how many cells make up a row? This is encoded in a headers file with the extension ".hdr", which you will find in your zip delivery as well. In my case there are the following entries in there:
| BYTEORDER | I |
| LAYOUT | BIL |
| NROWS | 5951 |
| NCOLS | 8066 |
| NBANDS | 1 |
| NBITS | 16 |
| BANDROWBYTES | 16132 |
| TOTALROWBYTES | 16132 |
| BANDGAPBYTES | 0 |
As you might have suspected already, the entries labeled NROWS and NCOLS deliver the required numbers. In my case, there are 5,951 rows each having 8,066 cells, for a total of 5,951 * 8,066 = 48,000,766 cells. Multiply this by 2 bytes per cell and you get 96,001,532 bytes - exactly the length of my BIL file.
Verification
So far we are quite confident that we are able to interpret the received data. What's missing is some actual programming. How you do this and in what language is up to you. I for one have chosen Linoleum (a cross-platform Assembler) as the programming language, and the pseudo code is quite trivial:
Open File
FilePos := 0
For Row = 1 To 5951
For Col = 1 To 8066
Elevation := Get 2 Bytes at FilePos
Do Something with Elevation
FilePos := FilePos + 2
Next Col
Switch to next Row
Next Row
Close File
The task described with "Do Something with Elevation" was to find the minimum and maximum elevations, so that I could quickly tell if the data I received was reliable. Remember that the 16 bits are signed, though (which makes sense, as it is very well possible to have elevations below sea level). In my case I received a minimum of -4 meters (0xFFFC), and a maximum of 3,053 meters (0x0BED). This is what we expected, since Maui's highest point is given as 3,055 meters.
The second task was more fun, as it was about creating a miniaturized false-color image. My current screen resolution is 1,280×1,024 pixels, so to accommodate the obtained data on my screen I needed to omit 7 out of 8 points, both horizontally and vertically, resulting in an image of 1,008×743 pixels.
I arbitrarily decided to output all points with an elevation of 0 and less meters as a blue pixel, and introduced height levels 512 meters apart, such subdividing the islands' heights into 7 layers (0 to 512 meters, 512 to 1024 meters, ..., 3072 to 3584 meters). Each layer was given a distinct plain color (yellow, light green, ..., white). To make steep slopes better visible, individual heights were interpolated between two neighboring colors, in a resolution of 1/256 per color component. The scaled-down result can be seen in Fig. 8.
Fig. 8: False-color image of data sample
The islands' shapes don't differ from what we expected (compare with the map from Google Maps, Fig. 1), also smaller features like Maui's Kealia Pond (near the Southern coast of the isthmus) are represented accurately. To show this was the goal of this first article within the series Digital Terrain Analysis.
