dem_geoid: Support EGM2008

NEWS

@@ -1,11 +1,14 @@
 *** RELEASE 2.4.2, Upcoming! ***
 
+ - dem_geoid
+   * Support the EGM2008 geoid.
+
  - stereo and parallel_stereo
    * Can handle multiple input images (multi-view stereo). The first
      image is set as reference, disparities are computed from it 
      to the other ones, and joint triangulation is performed.   
 
- - wv_correct:
+ - wv_correct
     * Corrects TDI of 16, 48, 56, and 64 (forward and reverse scan
       directions) for WV01, TDI of 8 (forward only) for WV01, and TDI
       of 16, 48, 64 (forward and reverse scan directions) for
@@ -43,7 +46,7 @@
     * Bugfix, longitude could be off by 360 degrees.
     * Robustness to large jumps in point cloud values.
 
-  - pc_align:
+  - pc_align
     * Ability to read and write CSV files having UTM data (easting,
       northing, height above datum).
     * Read DEMs in the ISIS cube format.

docs/book/tools.tex

@@ -805,8 +805,12 @@ \section{dem\_geoid}
 (correcting, if need be, for differences in dimensions between the DEM
 and geoid datum ellipsoids).
 
-Two geoids and one areoid are supported. The Earth geoids are: EGM96, referenced to the WGS84 datum ellipsoid (\url{http://earth-info.nga.mil/GandG/wgs84/gravitymod/egm96/egm96.html})
-and NAVD88, referenced to the NAD83 datum ellipsoid (\url{http://www.ngs.noaa.gov/GEOID/GEOID09/}).
+Three geoids and one areoid are supported. The Earth geoids are: EGM96
+and EGM2008, relative to the WGS84 datum ellipsoid
+(\url{http://earth-info.nga.mil/GandG/wgs84/gravitymod/egm96/egm96.html},
+\url{http://earth-info.nga.mil/GandG/wgs84/gravitymod/egm2008/egm08_wgs84.html})
+and NAVD88, relative to the NAD83 datum ellipsoid
+(\url{http://www.ngs.noaa.gov/GEOID/GEOID09/}).
 
 The Mars areoid is MOLA MEGDR
 (\url{http://geo.pds.nasa.gov/missions/mgs/megdr.html}). When importing
@@ -824,11 +828,13 @@ \section{dem\_geoid}
 \endlastfoot
 \hline
 Options & Description \\ \hline \hline
-\texttt{-\/-help|-h} & Display the help message\\ \hline
-\texttt{-\/-nodata-value \textit{integer(=-32768)}} & The value of no-data pixels, unless specified in the DEM \\ \hline
-\texttt{-\/-output-prefix|-o \textit{filename}} & Specify the output file prefix \\ \hline
-\texttt{-\/-double} & Output using double precision (64 bit) instead of float (32 bit)\\ \hline
-\texttt{-\/-reverse-adjustment} & Go from DEM relative to the geoid/areoid to DEM relative to the datum ellipsoid\\ \hline
+\texttt{-\/-help|-h} & Display the help message.\\ \hline
+\texttt{-\/-nodata-value \textit{integer(=-32768)}} & The value of no-data pixels, unless specified in the DEM. \\ \hline
+\texttt{-\/-geoid \textit{string}} & Specify the geoid to use for Earth WGS84 DEMs. 
+                               Options: EGM96, EGM2008. Default: EGM96. \\ \hline
+\texttt{-\/-output-prefix|-o \textit{filename}} & Specify the output file prefix. \\ \hline
+\texttt{-\/-double} & Output using double precision (64 bit) instead of float (32 bit).\\ \hline
+\texttt{-\/-reverse-adjustment} & Go from DEM relative to the geoid/areoid to DEM relative to the datum ellipsoid.\\ \hline
 \end{longtable}
 
 \section{dg\_mosaic}

src/asp/Tools/Makefile.am

@@ -83,6 +83,7 @@ endif
 if MAKE_APP_DEM_GEOID
   bin_PROGRAMS += dem_geoid
   dem_geoid_SOURCES = dem_geoid.cc
+  dem_geoid_LDFLAGS = $(PKG_GEOID_LDFLAGS) $(AM_LDFLAGS)
   dem_geoid_LDADD = $(APP_DEM_GEOID_LIBS)
 endif
 

src/asp/Tools/dem_geoid.cc

@@ -19,6 +19,11 @@
 /// \file dem_geoid.cc
 ///
 
+extern "C" {
+  void egm2008_call_interp_(int* nriw2, int* nciw2, double* grid,
+                            double* flon, double* flat, double* val);
+}
+
 #include <vw/FileIO.h>
 #include <vw/Image.h>
 #include <vw/Cartography.h>
@@ -30,16 +35,17 @@
 namespace po = boost::program_options;
 namespace fs = boost::filesystem;
 
-using std::endl;
-using std::string;
 using namespace vw;
 using namespace vw::cartography;
+using namespace std;
 
 template <class ImageT>
 class DemGeoidView : public ImageViewBase<DemGeoidView<ImageT> >
 {
   ImageT m_img;
   GeoReference const& m_georef;
+  bool m_is_egm2008;
+  vector<double> const& m_egm2008_grid;
   ImageViewRef<PixelMask<double> > const& m_geoid;
   GeoReference const& m_geoid_georef;
   bool m_reverse_adjustment;
@@ -53,9 +59,12 @@ class DemGeoidView : public ImageViewBase<DemGeoidView<ImageT> >
   typedef ProceduralPixelAccessor<DemGeoidView> pixel_accessor;
 
   DemGeoidView(ImageT const& img, GeoReference const& georef,
+               bool is_egm2008, vector<double> const& egm2008_grid,
                ImageViewRef<PixelMask<double> > const& geoid,
-               GeoReference const& geoid_georef, bool reverse_adjustment, double correction, double nodata_val):
+               GeoReference const& geoid_georef, bool reverse_adjustment,
+               double correction, double nodata_val):
     m_img(img), m_georef(georef),
+    m_is_egm2008(is_egm2008), m_egm2008_grid(egm2008_grid),
     m_geoid(geoid), m_geoid_georef(geoid_georef),
     m_reverse_adjustment(reverse_adjustment),
     m_correction(correction),
@@ -81,7 +90,7 @@ class DemGeoidView : public ImageViewBase<DemGeoidView<ImageT> >
     // Need to carefully wrap lonlat to the [0, 360) x [-90, 90) box.
     // Note that lon = 25, lat = 91 is the same as lon = 180 + 25, lat = 89
     // as we go through the North pole and show up on the other side.
-    while ( std::abs(lonlat[1]) > 90.0 ){
+    while ( fabs(lonlat[1]) > 90.0 ){
       if ( lonlat[1] > 90.0 ){
         lonlat[1] = 180.0 - lonlat[1];
         lonlat[0] += 180.0;
@@ -94,24 +103,35 @@ class DemGeoidView : public ImageViewBase<DemGeoidView<ImageT> >
     while( lonlat[0] <   0.0  ) lonlat[0] += 360.0;
     while( lonlat[0] >= 360.0 ) lonlat[0] -= 360.0;
 
-    result_type       geoid_height           = 0.0;
-    Vector2           pix                    = m_geoid_georef.lonlat_to_pixel(lonlat);
-    PixelMask<double> interp_val             = m_geoid(pix[0], pix[1]);
-    if (!is_valid(interp_val)) return m_nodata_val;
-    geoid_height                             = interp_val.child() + m_correction;
-    result_type       height_above_ellipsoid = m_img(col, row, p);
-    double            direction              = m_reverse_adjustment?-1:1;
+    result_type geoid_height = 0.0;
+    if (m_is_egm2008){
+      int nr = m_geoid.rows(), nc = m_geoid.cols();
+      egm2008_call_interp_(&nr, &nc, (double*)&m_egm2008_grid[0],
+                           &lonlat[0], &lonlat[1], &geoid_height);
+    }else{
+      Vector2  pix = m_geoid_georef.lonlat_to_pixel(lonlat);
+      PixelMask<double> interp_val = m_geoid(pix[0], pix[1]);
+      if (!is_valid(interp_val)) return m_nodata_val;
+      geoid_height = interp_val.child();
+    }
+
+    geoid_height += m_correction;
+
+    result_type     height_above_ellipsoid = m_img(col, row, p);
+    double          direction              = m_reverse_adjustment?-1:1;
     // See the note in the main program about the formula below
-    result_type height_above_geoid           = height_above_ellipsoid - direction*geoid_height;
-
+    result_type height_above_geoid         = height_above_ellipsoid - direction*geoid_height;
+    
     return height_above_geoid;
   }
 
   /// \cond INTERNAL
   typedef DemGeoidView<typename ImageT::prerasterize_type> prerasterize_type;
   inline prerasterize_type prerasterize( BBox2i const& bbox ) const {
     return prerasterize_type( m_img.prerasterize(bbox), m_georef,
-                              m_geoid, m_geoid_georef, m_reverse_adjustment, m_correction, m_nodata_val );
+                              m_is_egm2008, m_egm2008_grid,
+                              m_geoid, m_geoid_georef, m_reverse_adjustment,
+                              m_correction, m_nodata_val );
   }
   template <class DestT> inline void rasterize( DestT const& dest, BBox2i const& bbox ) const {
     vw::rasterize( prerasterize(bbox), dest, bbox );
@@ -122,27 +142,31 @@ class DemGeoidView : public ImageViewBase<DemGeoidView<ImageT> >
 template <class ImageT>
 DemGeoidView<ImageT>
 dem_geoid( ImageViewBase<ImageT> const& img, GeoReference const& georef,
+           bool is_egm2008, vector<double> & egm2008_grid,
            ImageViewRef<PixelMask<double> > const& geoid,
-           GeoReference const& geoid_georef, bool reverse_adjustment, double correction, double nodata_val) {
-  return DemGeoidView<ImageT>( img.impl(), georef, geoid, geoid_georef,
+           GeoReference const& geoid_georef, bool reverse_adjustment,
+           double correction, double nodata_val) {
+  return DemGeoidView<ImageT>( img.impl(), georef,
+                               is_egm2008, egm2008_grid,
+                               geoid, geoid_georef,
                                reverse_adjustment, correction, nodata_val );
 }
 
 struct Options : asp::BaseOptions {
-  string dem_name, out_prefix;
+  string dem_name, geoid, out_prefix;
   double nodata_value;
   bool use_double;
   bool reverse_adjustment;
 };
 
-std::string get_geoid_full_path(std::string geoid_file){
+string get_geoid_full_path(string geoid_file){
 
   // We try two ways of finding the path to the geoid file.
 
   // 1. Using the compile flag, from the GEOID_PATH macro.
 #define STR_EXPAND(tok) #tok
 #define STR_QUOTE(tok) STR_EXPAND(tok)
-  std::string full_path = std::string(STR_QUOTE(GEOID_PATH)) + "/" + geoid_file;
+  string full_path = string(STR_QUOTE(GEOID_PATH)) + "/" + geoid_file;
   if (fs::exists(full_path)) return full_path;
 
   // 2. Using the ASP_DATA path.
@@ -152,8 +176,7 @@ std::string get_geoid_full_path(std::string geoid_file){
               << "It should point to the 'share' directory of your ASP distribution.\n" );
   }
 
-  // Must keep this version synchronized with Packages.py!
-  full_path = std::string(asp_data) + "/geoids-1.1/" + geoid_file;
+  full_path = string(asp_data) + "/geoids/" + geoid_file;
   if (!fs::exists(full_path)){
     vw_throw( ArgumentErr() << "Could not find geoid: " << full_path << ".\n"
               << "The value of environmental variable ASP_DATA is '" << asp_data << "'.\n");
@@ -168,6 +191,7 @@ void handle_arguments( int argc, char *argv[], Options& opt ){
   general_options.add_options()
     ("nodata_value", po::value(&opt.nodata_value)->default_value(-32768),
      "The value of no-data pixels, unless specified in the DEM.")
+    ("geoid", po::value(&opt.geoid), "Specify the geoid to use for Earth WGS84 DEMs. Options: EGM96, EGM2008. Default: EGM96.")
     ("output-prefix,o", po::value(&opt.out_prefix), "Specify the output prefix.")
     ("double", po::bool_switch(&opt.use_double)->default_value(false)->implicit_value(true),
      "Output using double precision (64 bit) instead of float (32 bit).")
@@ -184,9 +208,9 @@ void handle_arguments( int argc, char *argv[], Options& opt ){
   po::positional_options_description positional_desc;
   positional_desc.add("dem", 1);
 
-  std::string usage("[options] <dem>");
+  string usage("[options] <dem>");
   bool allow_unregistered = false;
-  std::vector<std::string> unregistered;
+  vector<string> unregistered;
   po::variables_map vm =
     asp::check_command_line( argc, argv, opt, general_options, general_options,
                              positional, positional_desc, usage,
@@ -196,6 +220,8 @@ void handle_arguments( int argc, char *argv[], Options& opt ){
     vw_throw( ArgumentErr() << "Requires <dem> in order to proceed.\n\n"
               << usage << general_options );
 
+  boost::to_lower(opt.geoid);
+
   if ( opt.out_prefix.empty() )
     opt.out_prefix = fs::path(opt.dem_name).stem().string();
 
@@ -221,11 +247,12 @@ void handle_arguments( int argc, char *argv[], Options& opt ){
 // N = EGM96 Geoid height
 // Therefore: H = h - N.
 
-// We support three geoids: EGM96, NAVD88 (Earth) and MOLA MEGDR (Mars).
+// We support four geoids: EGM96, EGM2008, NAVD88 (Earth) and MOLA MEGDR (Mars).
 // Online tools for verification:
 // EGM96:  http://earth-info.nga.mil/GandG/wgs84/gravitymod/egm96/intpt.html
 // NAVD88: http://www.ngs.noaa.gov/cgi-bin/GEOID_STUFF/geoid09_prompt1.prl
 
+// EGM2008: http://earth-info.nga.mil/GandG/wgs84/gravitymod/egm2008/egm08_wgs84.html
 // MOLA MEGDR: http://geo.pds.nasa.gov/missions/mgs/megdr.html
 
 // Examples: At lat = 37 and lon = -121 (Basalt Hills, CA), we have
@@ -238,6 +265,12 @@ void handle_arguments( int argc, char *argv[], Options& opt ){
 // NAVD88 geoid height: 27.353
 // EGM96 geoid height:  26.58
 
+// For EGM2008, we use the tabulated values at 2.5', and the provided routine
+// for interpolation (interp_2p5min.f), as described at
+// http://earth-info.nga.mil/GandG/wgs84/gravitymod/egm2008/egm08_wgs84.html
+// We store the tabulated values as a tif file, and modified
+// the routine to be callable from this executable.
+
 int main( int argc, char *argv[] ) {
 
   Options opt;
@@ -261,34 +294,49 @@ int main( int argc, char *argv[] ) {
       vw_throw( ArgumentErr() << "Missing georeference for DEM: " << opt.dem_name << "\n" );
 
     // Find out the datum from the DEM. If we fail, we do an educated guess.
-    std::string datum_name = dem_georef.datum().name();
-    std::string lname = boost::to_lower_copy(datum_name);
-    std::string geoid_file;
+    string datum_name = dem_georef.datum().name();
+    string lname = boost::to_lower_copy(datum_name);
+    string geoid_file;
+    bool is_wgs84 = false;
     if ( lname == "wgs_1984" || lname == "wgs 1984" || lname == "wgs1984" ||
          lname == "wgs84" || lname == "world geodetic system 1984" ){
-      geoid_file = "egm96-5.tif";
+      is_wgs84 = true;
     }else if (lname == "north_american_datum_1983"){
       geoid_file = "navd88.tif";
     }else if (lname == "d_mars"){
       geoid_file = "mola_areoid.tif";
-    }else if ( std::abs( dem_georef.datum().semi_major_axis() - 6378137.0 ) < 500.0){
+    }else if ( fabs( dem_georef.datum().semi_major_axis() - 6378137.0 ) < 500.0){
       // Guess Earth
       vw_out(WarningMessage) << "Unknown datum: " << datum_name << ". Guessing: WGS_1984.\n";
-      geoid_file = "egm96-5.tif";
-    }else if ( std::abs( dem_georef.datum().semi_major_axis() - 3396190.0) < 500.0){
+      is_wgs84 = true;
+    }else if ( fabs( dem_georef.datum().semi_major_axis() - 3396190.0) < 500.0){
       // Guess Mars
       vw_out(WarningMessage) << "Unknown datum: " << datum_name << ". Guessing: D_MARS.\n";
       geoid_file = "mola_areoid.tif";
     }else{
       vw_throw( ArgumentErr() << "Cannot apply geoid adjustment to DEM relative to datum: "
                 << datum_name << "\n");
     }
+
+    bool is_egm2008 = false;
+    if (is_wgs84){
+      if (opt.geoid == "egm2008"){
+        is_egm2008 = true;
+        geoid_file = "egm2008.tif";
+      }else if (opt.geoid == "egm96" || opt.geoid == "")
+        geoid_file = "egm96-5.tif";
+      else
+        vw_throw( ArgumentErr() << "Unsupported value for geoid: " << opt.geoid << "\n");
+    }else if (opt.geoid != "")
+      vw_throw( ArgumentErr() << "The geoid value: " << opt.geoid
+                << " is applicable only for the WGS_1984 datum.\n");
+
     geoid_file = get_geoid_full_path(geoid_file);
     vw_out() << "Adjusting the DEM using the geoid: " << geoid_file << endl;
 
     // Read the geoid containing the adjustments. Read it in memory
     // entirely to dramatically speed up the computations.
-    double geoid_nodata_val = std::numeric_limits<float>::quiet_NaN();
+    double geoid_nodata_val = numeric_limits<float>::quiet_NaN();
     DiskImageResourceGDAL geoid_rsrc(geoid_file);
     if ( geoid_rsrc.has_nodata_read() ) {
       geoid_nodata_val = geoid_rsrc.nodata_read();
@@ -300,28 +348,47 @@ int main( int argc, char *argv[] ) {
     // Need to apply an extra correction if the datum radius of the geoid is different
     // than the datum radius of the DEM to correct. We do this only if the datum is
     // a sphere, such on mars, as otherwise a uniform correction won't work.
-    double major_correction
-      = geoid_georef.datum().semi_major_axis() - dem_georef.datum().semi_major_axis();
-    double minor_correction
-      = geoid_georef.datum().semi_minor_axis() - dem_georef.datum().semi_minor_axis();
+    double major_correction = 0.0;
+    double minor_correction = 0.0;
+    if (!is_egm2008){
+      major_correction
+        = geoid_georef.datum().semi_major_axis() - dem_georef.datum().semi_major_axis();
+      minor_correction
+        = geoid_georef.datum().semi_minor_axis() - dem_georef.datum().semi_minor_axis();
+    }
     if (major_correction != 0){
-      if (std::abs(1.0 - minor_correction/major_correction) > 1.0e-5){
+      if (fabs(1.0 - minor_correction/major_correction) > 1.0e-5){
         vw_throw( ArgumentErr() << "The input DEM and geoid datums are incompatible. "
                   << "Cannot apply geoid adjustment.\n" );
       }
       vw_out(WarningMessage) << "Will compensate for the fact that the input DEM and geoid datums "
                              << "axis lengths differ.\n";
     }
 
+    // The EGM2008 case is special. Then, we don't do bicubic interpolation into 
+    // geoid_img, rather, we invoke some Fortran routine, which gives more accurate results.
+    vector<double> egm2008_grid;
+    if (is_egm2008){
+      int nr = geoid_img.rows(), nc = geoid_img.cols();
+      egm2008_grid.resize(nr*nc);
+      for (int col = 0; col < nc; col++){
+        for (int row = 0; row < nr; row++){
+          egm2008_grid[row + col*nr] = geoid_img(col, row); // that is, egm2008_grid(row, col);
+        }
+      }
+    }
+
     ImageViewRef<PixelMask<double> > geoid
       = interpolate(create_mask( pixel_cast<double>(geoid_img), geoid_nodata_val ),
                     BicubicInterpolation(), ZeroEdgeExtension());
 
-    ImageViewRef<double> adj_dem = dem_geoid(dem_img, dem_georef, geoid, geoid_georef,
+    ImageViewRef<double> adj_dem = dem_geoid(dem_img, dem_georef,
+                                             is_egm2008, egm2008_grid,
+                                             geoid, geoid_georef,
                                              reverse_adjustment, major_correction, dem_nodata_val);
 
-    std::string adj_dem_file = opt.out_prefix + "-adj.tif";
-    vw_out() << "Writing adjusted DEM: " << adj_dem_file << std::endl;
+    string adj_dem_file = opt.out_prefix + "-adj.tif";
+    vw_out() << "Writing adjusted DEM: " << adj_dem_file << endl;
 
     if ( opt.use_double ) {
       // Output as double