Fun with VirtualPlanetBuilder

I've been playing with the VirtualPlanetBuilder tool for Open Scene Graph (OSG). This flexible program builds paged terrain databases with automatic level of detail that can be loaded into the standard osgviewer program or any other program that uses the standard OSG loaders. It takes heightfield data and imagery as input, and it accepts a variety of formats because it uses the GDAL geospatial data library. If the input data is georeferenced, the excercise of building a terrain database is almost too easy.

While VirtualPlanetBuilder can be used to build terrabyte-sized whole-earth databases, my own ambitions are so far more modest: build a small patch of terrain, drape some terrain over it, and fly around. Inspired by the vinyl relief maps that my parents mounted on the wall of our apartment when I was growing up, I decided that it would be fun to drape topographic maps over the terrain. I made a little geocentric database of the Mt. Washington area. I found this useful set of links to sources of digital elevation data.

For the DEM data, I went with webGIS, which provides DEMs based on USGS quadrangles. This turned out to not be terribly convenient, as there's state and county map interface but nothing that graphically shows the layout of the quadrangles. I flipped back-and-forth to the USGS web site and, armed with my knowledge of the White Mountains, downloaded the DEM data for the Mount Washington, Carter Dome, Crawford Notch, Stairs Mountain, Jackson, and Mount Dartmouth quadrangles. There's probably a more automated way to do this.

I downloaded the images of the map quadrangles using the usgsquads.com Goggle Earth interface, but this was a bit broken too. The map that contains Mt. Washington itself is 7.5 x 15 minutes, and the download link to it was broken. I went to the USGS EarthExplorer to grab that map. I would have been better off getting all the map images there. The maps at usgsquads.com seem to be pretty old and their elevations are in feet, whereas the map from USGS was in meters. Oh well, it still looks nice.

Because the .dem height data files and the TIFF map files are geolocated, it is quite straight forward, after a little preparation, to feed them to VirtualPlanetBuilder's osgdem program and get a data base. But my first effort had a problem; the maps have borders, and they show up on the terrain!

Whoops. Because the map images are projected in UTM, the latitude and longitude boundaries of the maps are not horizontal and vertical in the image, so it would be awkward to trim the maps in a program like the Gimp.

Fortunately osgdem needs its data in a geographic coordinate system and reprojects if necessary, so we can do the reprojection ourselves and trim the borders using a GDALutility.

That looks better, but the relief of the mountains seems exaggerated, just like those vinyl relief maps. A little snooping with gdalinfo shows that the elevation data is in feet, but osgdem builds a geocentric database using meters, and that info got lost somewhere along the way. Fortunately it's easy to fix this up with a GDAL tool.

Here's the complete recipe. First, many of the .dem files have the same name, so I renamed them in the directories where I unpacked them:

for i in *.zip; do
mkdir ${i%.zip}
mv $i ${i%.zip}
done


Then I combined them into one TIFF and built overviews:

gdalwarp 16*/*dem elev.tif
gdaladdo -r average elev.tif 2 4 8 16 32

The overviews are probably not necessary for such a tiny database, but it's recommended anyway. We have to fix the feet / meters problem:

gdal_translate -scale 0 1.0 0 0.3048 elev.tif elev-meters.tif


Next, fix up the map borders. The maps are encoded using index color, which doesn't doesn't work well with fancy interpolation methods, so expand them into RGB:

for i in carter-dome.tif jackson.tif stairs.tif k44071c3/k44071c3.tif ; do
gdal_translate -expand rgb $i ${i%.tif}-rgb.tif
done

and then use gdalwarp to reproject and trim. For example:

gdalwarp -t_srs "WGS84" -te -71.25 44.25 -71.125 44.375 \
-r cubicspline carter-dome.tif carter-dome-trim.tif

Build overviews:

for i in carter-dome-trim.tif crawford-trim.tif jackson-trim.tif \
stairs-trim.tif k44071c3/k44071c3-trim.tif ;do
gdaladdo -r average $i 2 4 8 16 32;
done

Finally, run osgdem:

osgdem --TERRAIN --geocentric -d elev-meters.tif -t carter-dome-trim.tif \
-t crawford-trim.tif -t jackson-trim.tif -t stairs-trim.tif \
-t k44071c3/k44071c3-trim.tif -l 8 -o wash.ive

and have fun!