Exploring RPCs

TL;DR: Generating stable RPCs is hard, and I'm not sure how to do it!

RPCs in the context of the geospatial world are equations for describing how a geometrically raw sensor image can be related to a location in the real world.  RPC can variously expand to Rational Polynomial Coefficient or Rapid Positioning Coefficients.  They consist of 80 coefficents for a ratio of two third degree polynomials (for x and y) plus 10 more values for normalizing image and geographic coordinates.

A couple background papers include one by Vincent Tao and others on Understanding the Rational Function Model (pdf), and one by Gene Dial and J. Grodecki titled RPC Replacement Camera Models (pdf).  These date from the mid-2000s, but they reference some earlier work by Gene Dial and J. Grodecki who seem to have played the seminal role in promoting the use of RPCs as generic replacements for rigerous camera models.  I believe the major satellite image vendors have now been delivering RPC models (example) for raw scenes for on the order of a decade. 

I, and others, have added support to GDAL for reading and representing RPCs as metadata over the a number of years, as well as support for evaluating RPCs in gdalwarp as a way of georeferencing images.  At it's very simplest a command like the following can be used to roughly rectify a raw image that has RPCs.

  gdalwarp -rpc raw.tif rectified.tif

In fact, at FOSS4G NA this year, a command like this was the big reveal in Paul Morin's discussion of how the Polar Geospatial Center was able to efficiently process huge numbers of satellite scenes of the polar regions.  I was rather surprised since I wasn't aware of this capability being in significant use.  However the above use example will usually do a relatively poor job.  RPCs encapsulate the perspective view of the satellite and as such you need to know the elevation of an early location in order to figure out exactly where it falls on the raw satellite image.  If you had a relatively flat image at a known elevation you can specify a fixed elevation like this:

  gdalwarp -rpc -to RPC_HEIGHT=257 raw.tif rectified.tif

However, in an area that is not so flat, you really need a reasonably could DEM (digital elevation model) of the region so that distortion effects from terrain can be taken into account.  This is particularly important if the view point is relatively low (ie. for an airphoto) or if the image is oblique rather than nadir (straight down).  You should be able to do this using a command referencing a DEM file for the area like this, though to be honest I haven't tried it recently (ever?) and I'm not sure how well it works:

  gdalwarp -rpc -to RPC_DEM=srtm.tif raw.tif rectified.tif

The above is essentially review from my point of view.  But with my recent move from Google to Planet Labs I have also become interested in how to generate RPCs given knowledge of a camera model, and the perspective view from a satellite.   This would give us the opportunity to deliver raw satellite scenes with a generic description of the geometry that can already be exploited in a variety of software packages, including GDAL, OSSIM, Orfeo Toolbox and proprietary packages like ArcGIS, Imagine and PCI Geomatica.  However, I had no idea how this would be accomplished and while I earned an Honours BMath, I have let the modest amount of math I mastered in the late 80's rot ever since.

So, using what I am good at, I started digging around the web for open source implementations for generating RPCs.  Unfortunately, RPC is a very generic term much more frequently used to mean "Remote Procedure Call" in the software development world.  Nevertheless, after digging I learned that the Orfeo Toolbox (OTB) has a program for generating RPCs from a set of GCPs (ground control points).  Ground control points are just a list of locations on an image where the real world location is known.  That is a list of information like (pixel x, pixel y, longitude, latitude, height).  That was exciting, and a bit of digging uncovered the fact that in fact Orfeo Toolbox's RPC capability is really just a wrapper for services used from OSSIM.  In particular core RPC "solver" is found in the source files ossimRpcSolver.cc, and ossimRpcSolver.h.  They appear to have been authored by the very productive Garrett Potts.

Based on the example of the Orfeo Toolbox code, I wrote a gcps2rpc C++ program that called out to libossim to transform GCPs read from a text file into RPCs written to a GDAL VRT file.

For my first attempt I created 49 GCPs in a 7x7 pattern.  The elevations were all zero (ie. they were on the WGS84 ellipsoids) and they were generated using custom code mapping from satellite image pixels to ground locations on the ellipsoid.  I included simple code in gcps2rpc to do a forward and reverse evaluation of the ground control points to see how well they fit the RPC.  In the case of this initial experiment the forward transformation, (long, lat, height) to (pixel, line) worked well with errors in the order of 1/1000th of a pixel.  Errors in the inverse transformation, (pixel, line, height) to (long, lat) was less good with worst case errors on the order of 1/10th of a pixel.

At first I wasn't sure what to make of the errors in the inverse direction.  1/10th of a pixel wasn't terrible from my perspective, but it was still a fairly dramatic difference from the tiny errors in the forward direction.  I tried comparing the results from evaluating the RPCs with OSSIM and with GDAL.  To test with GDAL I just ran gdaltransform something like this.

  gdaltransform -rpc -i simulated_rgb.vrt
  -85.0495014028 39.7243133317          (entered in terminal)
  2016.00011619623 1792.00003975369 0   (output in terminal)

and:
  gdaltransform -rpc simulated_rgb.vrt
  2016 1792                             (entered in terminal)
  -85.0495014329239 39.7243132735733 0  (output in terminal)

I found that GDAL and OSSIM agreed well in the forward transform, but while they both had roughly similar sized errors in the inverse they were still fairly distinct results.  Reading through the code the reasoning became fairly obvious.  They both use iterative solutions for the inverse transform and start with very different initial guesses, and they have fairly weak requirements for covergence.  In the case of GDAL the convergence threshold was 1/10th of a pixel.  OK, so that leaves me thinking there may be room to improve the inverse functions but accepting that the errors I'm seeing are really all about the iterative convergence. 

Next I tried using the RPCs with gdalwarp to actually "rectify" the image.  Initially I just left a default elevation of zero for gdalwarp, so I used a command like:

  gdalwarp -rpc simulated_rgb.vrt ortho_0.tif

The results seemed plausible but I didn't have a particularly good way to validate it.  My test scene is in Indiana, so I used a lat/long to height web service to determine that the region of my scene had a rough elevation of 254m above sea level (we will treat that as the same as the ellipsoid for now).  So I tried:

  gdalwarp -rpc -to RPC_HEIGHT=254 simulated_rgb.vrt ortho_254.tif

On review of the resulting image, it exactly overlayed the original.  I had been hoping the result would have been slightly different based on the elevation.  After a moments reflection and a review of the RPC which had many zero values in the terms related to elevation I realized I had't provided sufficient inputs to actually model the perspective transform.  I had't given any GCPs at different heights.  Also, the example code made it clear that an elevation sensitive RPC could not be generated with less than 80 GCPs (since the equations have 80 unknowns).

My next step was to produce more than 80 GCPs and to ensure that there was more than one elevation represented.  It seemed like two layers of 49 GCPs - one at elevation zero, and one at an elevation of 300m should be enough to do the trick, so I updated my GCP generating script to produce an extra set with an ellipsoid that was 300m "bigger".  I then ran that through gcps2rpc.

I was pleased to discover that I still got average and maximum errors at the GCPs very similar to the all-zero case.  I also confirmed that the GCPs at 300m were actually at a non-trivially different location.  So it seemed the RPCs were modelling the perspective view well, at least at the defining positions.  I thought I was home free!  So I ran a gcpwarp at a height of 0m, 300m and 254m.  The results at 0m and 300m looked fine, but somehow the results at 254 meters was radically wrong - greatly inflated!  A quick test with gdaltransform gave a sense of what was going wrong:

gdaltransform -rpc  simulated_rgb.vrt
0 0 0
-85.036383862713 39.6431112403042 0
0 0 300
-85.0363998100256 39.6431575414329 300
0 0 254
-84.9815303476089 39.2011322215184 254

Essentially this showed that pixel/line location 0,0 was transforming reasonably at the elevations 0 and 300, but at 254 the results were far from the region expected.  Clearly the RPCs were not generally modelling the perspective view well.  This is the point I am at now.  Is it not clear to me what I will need to do to produce stable RPCs for my perspective view from orbit.  Perhaps there is a more suitable pattern or number of GCPs to produce.  Perhaps the OSSIM RPC generation code needs to be revisited.  Perhaps I need to manually construct RPC coefficients that model the perspective view instead of just depending on a least squares best fit to produce reasonable RPCs from a set of GCPs. 

Anyways, I hope that this post will provide a bit of background on RPCs to others walking this road, and I hope to have an actually solution to post about in the future.

Getting Started with GeoTools

TL;DR: Use Eclipse->Import->Maven to build an eclipse project for GeoTools.

The project I work on at Google (Earth Engine) uses the GeoTools Java library, primarily for coordinate system transformations.  We have been using it for a few years for this purpose and it mostly works ok though we have worked around a few quirks.  I think it is desirable for our team to be able to contribute at least some small fixes back upstream to GeoTools and so I've been struggling to get to the point where that can be done smoothly.

Contribution Agreements

The first challenge was that we needed permission to contribute.  Generally this is not a problem as Google is pretty open source friendly, but in order to get the GeoTools corporate contribution agreement signed we needed to clear it with our legal folks.  They balked at signing the original agreement which was ... ahem ... hand crafted ... during the GeoTools incubation process.  At this point I was prepared to back off on the effort as I didn't want be the guy who raises a stink because his companies legal people were being difficult, as if that should make it the rest of the worlds problems.  However, +Jody Garnett was concerned that others might run into similar concerns and pushed through a GeoTools change request adopting a slight variation on the more widely accepted Apache contribution agreement.

The next step was internal, where my coworker +Eric Engle, who maintains our internal GeoTool tree, adjusted our use of GeoTools so we build from source and are in a position to patch our local version as  we improve things.  I had tried to do this once last fall and failed to overcome some challenges related to special requirements of manifest documents for GeoTools.  Something to do with finding factories for things like EPSG databases and such.

MapProjection Tolerances

This put things back in my court.  I had a small internally reproducable problem with the spherical mercator projection not staying within the error bounds expected by GeoTools at some points on the globe.  After some debugging Eclipse I came to the conclusion that:

1. It is evil that MapProjection.transform() actually uses Java asserts to confirm that point transformation are invertable within a tolerance.  This code only triggers when assertions are enabled, and makes it impossible to run code with transformations outside a reasonable realm without triggering this error when assertions are enabled. 
2. That the default getToleranceForAssertions() result is awfully tight - 1E-5 meters or 1/10th of millimeter!
3. For spherical projections it checks orthodromicDistance() which after stepping through the code I determined is very sensitive to rounding error for small distances.  Most locations near where I ran into a problem rounded to zero for distance in meters, while the slightly greater distance in radians at the problem point evaluated to something more like one centimeter, exceeding the error threshold. 

Well, those are the details, and I hope eventually to prepare a proper GeoTools bug report with a proposed patch.  In order to get to that I decided I needed to setup a "regular" GeoTools build that I could run tests against, and use to prepare a patch.  The GeoTools manual has great getting started information which I tried to follow.

Git Hate

GeoTools is hosted on github and I have a love/hate relationship with git.  Hmm, lets amend that.  I have a hate/hate relationship with git.  As I somewhat jestingly tell people I haven't yet got a full decade out of subversion yet after learning it in the middle of last decade so I start with a grudge about being forced to learn a new version control technology while the old one continues to work fine from my perspective.  But I also find git way more complicated than necessary.  Ultimately I determined that I needed to clone geotools on github so I would have a public repository that I could stage changes on and send pull-requests from since I can't issue pull requests for repositories living on my workstations at Google or at home.  Then I cloned that repository locally at work and eventually at home.  OK, I struggled a bit before concluding that is what was needed but I'm ok with that. 

Next I came to the conclusion that there is no easy way from the github interface to refresh my github clone from the official git instance, and I had to setup "upstream" on my workstation clone with the command:

  git remote add upstream http://github.com/geotools/geotools.git

And now every time I want my geotools clone on github updated I need to issue the following commands locally in my local github clone - basically pulling the changes down to my local repository and then pushing it back up to my own clone:
  git fetch upstream
  git merge upstream/master
  git push origin master

Blech.

I think things would be a bit different if I had write permission on the GeoTools git repository as then I'd be able to make my changes in a branch in it and then send pull requests against that but that would mean lots of branching which I'm also not too fond of.

Maven

GeoTools uses Maven for building.  I installed the Maven Ubuntu package and seemed good to go.  All I should need to do is "mvn install" in the GeoTools root directory to build.

The default Java on my workstation is a (possibly googlized) OpenJDK 1.7.  I get the impression from Jody that GeoTools is targetted at JDK 1.6 and that most folks use the Oracle JDK, not OpenJDK.  I was able to able to revert to using OpenJDK 1.6 but was still getting lots of things missing (I wish I had kept better notes on this part).  The net result after searching around was that Maven has "profiles" instructing it how to find some classes in different environments that seems to mostly affect the generation of Java docs.  After some messing around we found that because my OpenJDK returned "Google Inc." as the vendor something special was needed and Jody was kind enough to incorporate this upstream.

After this change I was able to compile a Maven build of GeoTools successfully, but it took hours and hours for reasons that were never clear and likely somehow Google specific.  Today I decided to try this again at home on my fairly stock Ubuntu 10.04 machine.  While not exactly fast, it was much faster and I could compile a "mvn install" on an already built tree without testing in a couple minutes.  I will note this time I used a downloaded Maven instead of the one provided by Ubuntu, but trying that at work did not seem to help.

Eclipse

Next, I wanted to get to the point where I had an Eclipse project I could use to debug and improve GeoTools with.  I tried to following the Eclipse Quick Start docs and learned a few things.

1. Stock Eclipse 3.8 does not seem to include any of the Maven support even though there is apparently some in Eclipse 3.7.  I did find downloading and installing the latest Eclipse (Juno - 4.something) did have good Maven support.
2. The described process produces a project where you can use GeoTools as jars, and even debug in, but not modify the GeoTools source.  This is ok for a use-only GeoTools user but not for someone hoping to fix things or even make temporary modifications during debugging.

I stumbled around quite a bit trying to get past this.  What I discovered is that the Eclipse File->Import->Maven->Existing Maven Projects path can be used to setup a project for one or more GeoTools modules.  Select the root directory of the GeoTools source tree and Eclipse will walk the pom.xml files (Maven "makefiles) to figure out the dependencies and where the source is.   You will need to select which ones to include.

I ran into complaints about lacking Eclipse understanding of  git-commit-id-plugin and something about jjtree.  I was able to work around this by not including the "cql" module which needed the jjtree stuff, and by hacking out the entire git-commit-id-plugin from the GeoTools root pom.xml file.  After that I also restricted myself to selecting the main "library" and I think "ogc" modules.  

The result was an Eclipse project in my workspace for each GeoTools module, including the referencing module I was interested in.  I was able to go into the src/test tree and run the junit tests fairly easily and all succeeded except three calls to *.class.desiredAssertionStatus() which I commented out.  They did not seem terribly relevant. 

Conclusion

Well, I can't say that I have accomplished my first meaningful contribution to GeoTools yet, but I have at least gotten to the point where I can build and work on things.  Despite some great getting start docs (mostly by Jody I think) and direct support from Jody in IRC (#geotools on freenode) it was still a painful process.  This is partly due to my relative lack of familiarity with Maven, Github, and even Eclipse outside of the very narrow way I use it at Google.  But it also seems that GeoTools is a big project and there are lots of traps one can fall into with versions of the various tools (Eclipse, Maven, Java runtime, etc).  But hopefully people will see that if I can hack it, you likely can too.