Winging a DEM for a mission using World View 1

The group I work for at NASA has a big robot that likes to drive in a quarry at speed. Doing this is risky as we could easily put the robot in a position to hurt itself or hurt others. One of things we do to mitigate the risk is by having a prior DEM of the test area. The path planning software can then use the DEM to determine where it is and what terrain is too difficult to cross.

Since ASP recently gained the ability to process Digital Globe and GeoEye imagery (more about that in a later post), I was given a request to make a DEM from some World View 1 imagery they purchased. The location was Basalt Hills, a quarry at the south end of the San Luis Reservoir. To process this imagery with any speed, it is required to map project the imagery on some prior DEM. My choices were SRTM or NED. In my runs, both DEMs have problems. SRTM has holes in the middle of it that needed to be filled so ASP would work correctly. NED had linear jumps in it that ASP couldn’t entirely reverse in its math.

I ended up using SRTM as a seed to create my final DEM of the quarry. If you haven’t seen this, the process looks like the following commands below in ASP 2.0+. What’s happening is that ASP uses an RPC map projection to overlay the imagery over SRTM. When it comes time for ASP to triangulate, it reverses math it used to map project, and then in the case of Digital Globe it will triangulate using the full camera model. Another thing worth noting is that ASP needs control over how the interpolation is performed when doing RPC map projection. This forces us not to use the GDAL utilities during this step and instead use our own custom utility.

parallel rpc_mapproject --tr 0.5 \
      --t_srs'"+proj=utm +zone=10 +datum=WGS84 +units=m +no_defs"' \
      filled_srtm_dem.tif {} {.}.XML {.}.srtm.crop.tif ::: left.TIF right.TIF
stereo left.srtm.crop.tif right.srtm.crop.tif left.XML right.XML \
      r1c1_srtm_crop/r1c1_srtm_crop filled_srtm_dem.tif

Afterwards we got a pretty spiffy result that definitely shows more detail than the prior DEM sources. Unfortunately the result was shifted from the NED DEM source that my crew had previously been using. This ideally would be fixed by bundle adjusting the World View camera locations. It was clearly needed as most of our projected rays only came within 3 meters of each other. Unfortunately ASP doesn’t have that implemented.

EDIT: If I had paid closer attention to my data I would have noticed that a large part of the differences I was seeing between my DEM and USGS’s NED was because the NED data uses a vertical datum. My ASP DEM are referenced against the WGS84 ellipsoid. NED data is referenced against WGS84 plus the NAVD88. This would account for a large part of the 30 meter error I was seeing. (11.19.12)

My “I’m-single-with-nothing-going-on-tonight” solution was the Point Cloud Library. It has promising iterative closest point (ICP) implementations inside it and will eventually have the normal distribution transform algorithm in it. It also has the benefit of having its libraries designed with some forethought compared to the hideous symbol mess that is OpenCV.

PCL's pcd_viewer looking at the quarry.

I achieved ICP with PCL by converted my PC (point cloud) file from ASP into a PCL PCD file [1]. I also converted the NED DEM into a PCD file [2]. I then subsampled my ASP point cloud file to something more manageable by PCL’s all-in-memory tactics [3]. Then I performed ICP to solve for the translation offset I had between the two clouds [4]. My offset ended up being about a 40 meter shift in the north and vertical direction. I then applied this translation back to the ASP PC file [5] so that the DEM and DRG could be re-rendered together using point2dem like normal.

I wrote this code in the middle of the night using a lot of C++ because I’m that guy. Here’s the code I used just for reference in the event that it might help someone. Likely some of the stuff I performed could have been done in Python using GDAL.


After rendering a new DEM of the shifted point cloud, I used MDenoise to clean up the DEM a bit. This tool is well documented at its own site (

I’ve also been learning some QGIS. Here are some screen shots where you can see the improved difference map between NED and my result after ICP. Generally this whole process was very easy. It leaves me to believe that with some polish this could make a nice automated way to build DEMs and register them against a trusted source. Ideally bundle adjustment would be performed, but I have a hunch that the satellite positioning for Earth targets is so good that very little shape distortion has happen in our DEM triangulations. I hope this has been of interest to some of you out there!

Difference map between the USGS NED map and ASP's WV01 result.

New CTX Layer in Google Earth

The Google community silently released a few new features for their Mars mode in Google Earth. MER-B, Opportunity, now has an updated traverse path thanks to a fellow at the Unmanned Spaceflight forum along with a new base map that Ross created. However I’m really excited about a new CTX Global Map that is available. Below is a screen shot:

From this high view you can see that that CTX hasn’t imaged all of Mars. This is expected. CTX isn’t trying to image all of Mars, it is simply the context imager for HiRISE. Meaning that CTX tends to roll tape only when HiRISE is. Never the less, the imagery is still beautiful and, in my opinion, shows more detail that the default base layer from an HRSC composite.

This mosaic was created by simply downloading all CTX data from NASA’s PDS servers and then calibrating and map projected the imagery with USGS’s ISIS software. The composite was then made with in house software from our team at NASA Ames Research Center. This is a Vision Workbench combination plus our Plate Filesystem. It is the same software we used to do a global MOC-NA and HiRISE mosaic for Microsoft’s World Wide Telescope. Unfortunately, I believe those servers have bit the dust. Regardless, if you have free time, I encourage you to check out the CTX Mosaic and Opportunity traverse in Google Earth. Click the planet, select “Mars”, and then in the bottom left select “CTX Mosaic” under “Global Maps”.

Ames Stereo Pipeline and Earth

I’ve been working on Ames Stereo Pipeline for about 4 years now. During that time the developers and I have always been focusing on what we could do with data that has come from NASA satellites. So we’ve produce 3D models of the Moon, Mars, and of the moons of Saturn. Recent work however has allowed us to produce 3D models of data much closer to home.

Earth example from ASP

In the next ASP release, I’m happy to announce that we will be providing limited support for processing images from Digital Globe satellites. Digital Globe sells images of Earth via World View and Quick Bird satellites. Their imagery is commonly seen in the base layer of Google Earth. ASP was recently awarded a grant to add support for those images to allow mapping of Earth’s polar regions.

Just last week I was able to produce ASP’s first 3D model of Earth using World View data. The image is a screenshot of the colormap DEM output draped in Google Earth. The area is just west of Denver in the Rockies and the imagery represented a huge challenge for ASP. There were lots of trees that make template correlation difficult. The mountains also made for a massive search range, which allowed outliers to make their way into the final result. Even though, we were still able to produce a result. I’m very proud of this result as it leaves me with confidence that ASP will do well in the polar regions where the terrain is not as aggressive.