more VirtualPlanetBuilder fun

Let's add some other imagery to our Mt. Washington terrain. The digital orthophoto quadrangles published by the USGS are a decent source of free aerial photographs of the U.S. These images are rectified, which means that any effects from the tilt of the camera and relief of the terrain have been removed; just the thing for draping over a DEM. Unfortunately they are in black-and-white, but the resolution is 1 meter, which isn't bad. The image files are georeferenced, very convenient for building a mosaic with VirtualPlanetBuilder because we don't need to specify any additional projection or position parameters with each image. Generally one topographic map quad is divided into 4 DOQs. This makes for a marathon download session at the USGS website; maybe someone knows a better way?

It's straightforward to build a new database using an osgdem command similar to that used to build the map overlay database. One complication is that we now have 24 image files to deal with, instead of 5, and I don't feel like typing all the names on the command line. The downloaded files from the USGS have the form O4407145.NES.837049.tgz, and their only contents seems to be another tar file. Extracting that, we get the TIFF file and a header file. I arranged these all in directories; for example, the path to one TIFF file is O4407145.NES/O4407145.NES.837049.tif. We generate the command line with a little Unix magic:

ls O*.*/*.tif | sed -e's;.*;-t &;' \
| xargs osgdem --TERRAIN --geocentric -d elev-meters.tif \
-l 8 -o wash-photo.ive


We have a nice image when we look at the scene with osgviewer, although the gray imagery would benefit from colors and more contrast. But what would be very cool is display the topographic map tiles and the orthophotos on the same scene and transition between the two. VirtualPlanetBuilder lets you specify that images belong to different layers, and the OSG example program osgmultitexturecontrol supports transitioning between layers based on elevation above terrain. In order to build the new database, we need to specify a "--layer 1" argument before each orthophoto:

ls O*.*/*.tif | sed -e's;.*;--layer 1 -t &;' |\
xargs 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-layer.ive

Phew! Unfortunately, we don't see the map layer at all with the osgmultitexturecontrol program from the OSG distribution. A look at the code reveals the problem: the elevation at which the transition occurs is hard-wired at 1000 kilometers. That's fine for viewing a whole-earth database from space, but doesn't do us any good. I hacked in an elevation transition argument, and got some better results. We start with the map view:

and as we zoom in, we transition from this:

to this:

Pretty cool. It would be nice to spruce the photos up a bit, perhaps by using the approach of merging color from lower resolution LANDSAT imagery with these high resolution DOQs.

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!