SRTM-3 data has proved to be one of the most popular and utilized GIS data sets of all time. Unfortunately, the quality of the data is less than ideal, the main problem being null data areas. For most files, these areas are relatively small and can be handled successfully by a straightforward interpolation of the base file. However, when the null data areas are more pronounced, another approach is required. Ideally, data from an alternative data set can be used to "patch" the null data areas. Unfortunately, there are not many sources for patch data. ( If there were, SRTM would not be so popular in the first place.)
Two possible sources of such data include NIMA DTED0 and SRTM30 data from the NASA SRTM data set itself. However, both these data sets are of very low resolution compared to SRTM-3 data, so any use of them for patching purposes must be done with great care in order to get acceptable results. Even with the best care the results will not be very precise, but they can at least be accurate.
BLACKART version 3.81 features greatly improved ability to handle DTED0 and SRTM30 data. This tutorial walks through two techniques for patching SRTM-3 files to demonstrate how this is done, and also to illustrate the differences in the results using the two sets. DTED0 data has been discussed at length in previous articles and I will not repeat it here. SRTM30 is a near-global digital elevation model (DEM) comprising a combination of data from the Shuttle Radar Topography Mission and the U.S. Geological Survey's GTOPO30 data set. It can be considered to be either an SRTM data set enhanced with GTOPO30, or as an upgrade to GTOPO30. SRTM30 data represents the best 30 Arc-second posting data currently available, surpassing DTED0 and GTOPO30 in data quality. This is the best of the three to use for patching.
My example uses SRTM-3 file N38W084 and the corresponding DTED0 file N38.dt0 and SRTM30 file w100n40.hgt. Before I get started, a comment about file naming conventions. The SRTM-3 file name includes the latitude and longitude of the Southwest corner of the tile. DTED0 tiles are stored in folders that are named by their longitude coordinate. The folders then contain files that are named and listed sequentially by the latitude coordinate. DTED0 tiles are the same size (1 degree by 1 degree) as SRTM-3 tiles. This makes selecting and matching SRTM-3 and DTED0 pretty easy. Just choose the single tile in the longitude folder of interest with the same lower left coordinate as your SRTM-3 tile.
The SRTM30 tiles, on the other hand use a different naming convention. They are formatted and organized in a fashion that mimics the GTOPO30 convention. Unfortunately, this makes the file naming different from the SRTM standard and care must be taken in its use. The most salient differences are that the latitude and longitude coordinates are reversed as compared to SRTM-3 files, i.e. with longitude first and latitude second. The file name coordinate identifies the upper left (Northwest) corner of the tile, NOT the lower left (Southwest) as in SRTM tiles. The 27 tiles that individually cover 50 degrees of latitude and 40 degrees of longitude each have 6,000 rows and 4,800 columns.
This makes selecting the proper SRTM30 coordinate a little trickier. In order to find the minimum latitude coordinate subtract (or add, if the NW coordinate is in the southern hemisphere and is considered negative) 50 degrees from the latitude. To get the other longitude coordinate, add (or subtract if east of the prime meridian) 40 degrees. It is necessary to map out the SRTM30 tile corners in order to make sure that your SRTM-3 tile lies within its boundaries.
For example, we wish to determine which SRTM-30 tile contains the SRTM-3 tile. First, determine the coordinates of your SRTM-3 tile. They are 38 degrees north latitude and 84 degrees west longitude for my example. The SRTM30 tile that contains the same part of the globe as my SRTM-3 tile is w100n40.hgt. I know this is true because if I subtract 50 degrees from the latitude I get s10 W 100. Since the latitude coverage is -10 to +40, I know that my SRTM-3 latitude of +38 is within parameters. Recall that west longitudes decrease to the east, so I know that the longitude spans -100 to -60. Since my target longitude of -84 is between -100 and -60, I know that it is inside the tile as well.
First Technique - Decoupled Interpolation and Substitution
OK, so lets patch a tile. Load BLACKART and open the SRTM-3 file first by selecting File|Open from the main menu (don't forget to select type .hgt from the drop down box). This tile will display in the graphics window. Let's patch with a DTED0 file first, because it is a little easier. Select File|Open DTED0 Merge DEM. Open n38.hgt that has been previously down loaded from the NIMA site. This file will display in a small (121 by 121 pixel) graphic window when it loads.
We are going to do the same thing for both types of patching data. First, we will expand the file to match the SRTM-3 file size. Then we will interpolate it to compute the missing data. Finally, we will substitute the interpolated data for the null data in the target file. (Note: we will not interpolate the target SRTM-3 file at all using this procedure.)
When both files are open, select File|Input Data from the graphics menu. You will be presented with a dialog screen. Enter your desired Laplacian and LSQR iterations as before. Now select Image|Array Operations|Expand DTED0 DEM. BLACKART will expand the 121 by 121 DTED0 file to a sparse 1201 by 1201 file. A graphics window will show you the result of the expansion. It will appear as a grid of dots that form the general outline of the source DTED0 file.
What has happened is that BLACKART has just spread out the data from the DTED0 file evenly across a 1201 by 1201 array. It does this by starting at the upper left corner of the array, and puts the first data point from the DTED0 file at coordinate [0,0] of the larger array. Instead of putting the next point at [0,1], it skips nine spaces and puts the next point at [0,10]. It proceeds in this fashion, skipping rows and columns until it creates a sparse 1201 by 1201 array using the data from the DTED0 file. In between it places zeros.
Now it is time to interpolate the file. Select Run|Run Blackart as before. BLACKART will interpolate your expanded DTED0 file (the "dot" file) to a smooth elevation surface. You may inspect the result of the interpolation by selecting View|Display Image from the main menu. You can check the interpolated elevations by dragging the mouse cursor over the file and looking at the elevations reported at the bottom of the main window. You may also save this interpolated image in SRTM or any of the formats supported by BLACKART.
If the interpolated DEM is satisfactory, you can merge the DEMs by selecting Image|Array Operations|Assign Interpolated Elevations BLACKART will inspect the SRTM file, replacing any null elevations with the corresponding elevation from the iterated DEM file. When this occurs, any null data values, which appear bright orange in the target file will generally turn blue or black. Of course the rest of the data values in the target file are left untouched by the substitution. You can check the results visually, or by querying using the mouse cursor.
Now we can try the same interpolation using an SRTM30 file as the source map. Open the SRTM-3 file as before. Now open the SRTM-30 file. This is a large, 58MB file and it will take a moment to open. BLACKART will automatically find and extract the small 120 by 120 subset of data that it needs for the interpolation and then it will close the file and deallocate the internal array. When both files are open, select File|Input Data from the graphics menu. You will be presented with a dialog screen. Enter your desired Laplacian and LSQR iterations as before. Now select Image|Array Operations|Expand SRTM30 DEM. BLACKART will expand the 120 by 120 SRTM30 file to a sparse 1201 by 1201 file. A graphics window will show you the result of the expansion.
Select Run|Run Blackart to interpolate the file as before. BLACKART will interpolate your expanded SRTM30 file to a smooth elevation surface. You may inspect the result of the interpolation by selecting View|Display Image from the main menu. You can check the interpolated elevations by dragging the mouse cursor over the file and looking at the elevations reported at the bottom of the main window. You may also save this interpolated image in SRTM or any of the formats supported by BLACKART.
If the interpolated DEM is satisfactory, you can merge the DEMs by selecting Image|Array Operations|Assign Interpolated Elevations BLACKART will inspect the SRTM file, replacing any null elevations with the corresponding elevation from the interpolated SRTM30 file.
Let's look at the results of our interpolation. The first image at the upper right shows the target N34W88.hgt target image. The null data areas are clearly visible in bright orange in the second image below. The third image shows the interpolated DTED0 data. The next image shows the SRTM30 source file. The next image shows a Global Mapper rendering of the DTED0 patched file.
The next image shows the target with the SRTM30 patch data substituted in. Comparison of the two indicates that the SRTM30 patch is somewhat superior to the DTED0 patch, but both seemed to work acceptably well.
Trouble in GIS-ville
That seemed to work OK. Let's try our technique on a more difficult map. SRTM S78W075 is a file that has been used as an example on the AVSIM Scenery and Mesh Design forum. It is a particularly tough target characterized by large regions of missing data. There are a lot of phase unwrapping anomalies due to the large expanses of water. Moreover, it is an archipelago that has a lot of steep cliffs dropping off into the water. Previous experience has shown that this is one of the toughest interpolation challenges. As we will see, we will not be disappointed (in the challenge, that is).
Let's try the technique that we experimented with above, using a DTED0 file to patch the target. I loaded S48W075.hgt and S48.dt0. I then expanded, interpolated and merged as before. The result is shown to the right.
What happened? The result is a bad fit between the interpolated patches and the target. There is substantial error between the interpolated surface and the SRTM-3 surface and the edges of the patched areas do not blend smoothly with the target. This example highlights some problems that can accompany interpolations of low-resolution data sets. The first problem is a result of the patch interpolation not being of high enough quality. More LSQR interpolations would probably have helped in this case. However, the edge matching problem highlights a key weakness of the interpolative approach taken up until now, which prompted me to add a different patching strategy to BLACKART.
Second Technique - Coupled Interpolation and Substitution
The edge mismatch problem is related to the decoupled nature of the interpolative patch and the target image. In the procedure outlined above, the DTED0 data was interpolated and then patches were cut out of the interpolated surface and then just inserted into the target image wherever there were null data areas. There was no opportunity for information to flow from the target map to the interpolated map. In addition to edge mismatches, this does not allow "good" information from the SRTM-3 null data edges to flow into the patch interpolation. The result was abrupt and unsightly inconsistencies between the two.
A better approach might be as follows: expand the DTED0 image as before. Rather than interpolate this image, let's just cut the expanded grid and insert the sparse grid into the target map. Now make a copy of the target with the grid points inserted and interpolate it, subject to certain constraints on the known data points. This allows information from the DTED0 map to combine with information from the target map to produce a better interpolated surface. Now, cut sections from this improved interpolated surface and insert them back into the target map wherever null data areas exist. We would expect a much better edge match using this technique and in fact this appears to be the case. I modified BLACKART to do just that (as well as the other technique).
In order to execute this method, load S48W075.hgt as before. Then load DTED0 file S48.dt0 also as before. Next, expand the DTED0 file by selecting Image|Array Operations| Expand DTED0 DEM. Now, insert the DTED0 data points into the target file by selecting Image|Array Operations| Insert Expanded DTED0 or SRTM30 into SRTM. If you look closely at the image, you should see an array of dots across the null data areas wherever the DTED0 file was not zero. Now select File|Input Data from the graphics window menu as before. Let's try 100 LSQR and 1000 Laplacian iterations. Select Run|Run Blackart from the main menu. After processing, BLACKART will update the null data areas with the computed patches.
Although the results look promising, there are still a few problems with the file left to solve. The first is the spurious noise on the water. Notice also (if you are following along at home) that there is new noise that has been introduced in some of the valleys and fjords where null data formerly existed. As far as I can determine, this is an artifact of the LSQR solver. Apparently the cliffs dropping off sharply into the water or flat valleys causes ringing at the discontinuities. This ringing induces negative elevations, since the flat surfaces are at or very close to zero elevation and each sinusoidal valley must therefore be negative by definition. These negative values can then propagate after each iterative pass. (For a more detailed explanation of this phenomenon, with a salient example, see my Master's Project Report (Crater Lake Interpolation section) listed elsewhere on this site for download.)
There are a couple of solutions to this problem. The easiest is to remove the unwanted negative artifacts by selecting Image|Array Operations|Set Sea Level Clipping Parameters. Set the values to -200 to 0 for this example. Now select Image|Array Operations|Fix Sea Level Elevations. All elevation values between 0 and -200 are set to zero, which accounts for all of the remaining bad data. The result of our interpolation is shown at the right. It is far from perfect, but looks a lot better than before and may be serviceable for some purposes.
Further tests not shown in this article demonstrated that this technique yielded better results than either independent interpolation and patching or straight interpolation of the SRTM-3 file. In the latter case, an entire island was missing from the map! This island existed entirely in a null data area. If not for the DTED0 data, BLACKART would have no way of knowing that it was even there! The final image shows the results from the straight interpolation of the SRTM-3 file. (Note the missing mountain inside the red box on the full size image.)
The same technique can be used for SRTM30 data, except that you must expand the SRTM30 file rather than the DTED0 file by clicking the appropriate menu selection.
Although this tutorial examined some more complicated but not uncommon situations, let me repeat the point that in most cases, straightforward interpolation will provide the easiest and perhaps the best result. Only when the null data areas get sufficiently large (like for example greater than the RMS of the surface variation) are more advanced techniques required. Interpolation is a very fascinating black art indeed and I could probably write several more pages of interesting examples on the subject, but I do not want to take away all your fun. It should be clear that no one interpolation technique works best under all conditions. It is my objective to provide as many tools in the BLARKART utility to address as may patching and interpolation applications as possible.
Editor's Note: Since writing this tutorial, I added a couple of additional features to the program. It is now possible to save your result as a gray scale TGA file or in POV-Ray heightfield compatible TGA format. User input Refining iteration capability was also added. Refining takes place after the computation of the interpolated DEM by the Laplacian solver. It relaxes the known elevations and allows them to be interpolated along with the unknown elevations. This has the effect of removing many unwanted artifacts such as pimples and dimples produced by the known DTED0 or SRTM30 grid elevations. There is no single correct number of Refining iterations to run. The greater the number, the smoother the surface. However, this comes at the expense of fidelity, because with every iteration the known elevations are collapsing from their "correct" values to other values. The default ten is a good compromise.
BLACKART can also read, interpolate and save DTED1 files. Although this data is not commonly available to the civilian public, it is sometimes available to government contractors. It is also possible to use DTED1 data to patch SRTM-1 files. The program expands the mesh and then inserts the expanded data where needed into the SRTM-1 file. Then the large SRTM-1 file can be interpolated as usual.