Ok, this change seems to help a lot.

There are quite a few difficult cases where an airport is located on top
of a hill, or in a bowl, or has a significant elevation change near by.
I take the average elevation of the area and clamp the outlyers.  However
these difficult cases "bias" the average elevation because the airport
surface may include much of the surrounding area.

This change to the code computes the airport elevation *only* based on the
actual airport geometry node and ignores all the surrounding nonsense that
might exist.  This doesn't make things perfect, but is a *big* step forward
for airports in areas with significant elevation change nearby.

src/Airports/GenAirports/Makefile.am

@@ -29,6 +29,7 @@ genapts_SOURCES = \
 	apt_surface.hxx apt_surface.cxx \
 	build.cxx build.hxx \
 	convex_hull.cxx convex_hull.hxx \
+	elevations.cxx elevations.hxx \
 	global.hxx \
 	lights.hxx lights.cxx \
 	main.cxx \

src/Airports/GenAirports/apt_surface.cxx

@@ -33,12 +33,14 @@
 
 #include <Array/array.hxx>
 
+#include "elevations.hxx"
 #include "global.hxx"
 #include "apt_surface.hxx"
 
 SG_USING_NAMESPACE( PLib );
 
 
+#if 0
 // fix node elevations.  Offset is added to the final elevation,
 // returns average of all points.
 static double calc_elevations( const string &root, const string_list elev_src,
@@ -164,17 +166,20 @@ static double calc_elevations( const string &root, const string_list elev_src,
 
     return average;
 }
+#endif
 
 
 // Constructor, specify min and max coordinates of desired area in
 // lon/lat degrees
 TGAptSurface::TGAptSurface( const string& path,
                             const string_list& elev_src,
-                            Point3D _min_deg, Point3D _max_deg )
+                            Point3D _min_deg, Point3D _max_deg,
+                            double _average_elev_m )
 {
     // Calculate desired size of grid
     min_deg = _min_deg;
     max_deg = _max_deg;
+    average_elev_m = _average_elev_m;
 
     // The following size calculations are for the purpose of
     // determining grid divisions so it's not important that they be
@@ -220,8 +225,12 @@ TGAptSurface::TGAptSurface( const string& path,
         }
     }
 
-    // Determine elevation of the grid points
+    // Lookup the elevations of all the grid points
-    double average = calc_elevations( path, elev_src, dPts );
+    tgCalcElevations( path, elev_src, dPts );
+
+    // Clamp the elevations against the externally provided average
+    // elevation.
+    tgClampElevations( dPts, average_elev_m, max_clamp );
 
     // Build the normal res input grid from the double res version
     Matrix_Point3Dd Pts(ydivs + 1, xdivs + 1);
@@ -259,7 +268,7 @@ TGAptSurface::TGAptSurface( const string& path,
 
     bool slope_error = true;
     while ( slope_error ) {
-        // cout << "start of slope processing pass" << endl;
+        SG_LOG( SG_GENERAL, SG_DEBUG, "start of slope processing pass" );
         slope_error = false;
         // Add some "slope" sanity to the resulting surface grid points
         for ( int i = 0; i < Pts.cols() - 1; ++i ) {
@@ -280,11 +289,11 @@ TGAptSurface::TGAptSurface( const string& path,
 
                     slope_error = true;
 
-                    // cout << " (a) detected slope of " << slope
+                    SG_LOG( SG_GENERAL, SG_DEBUG, " (a) detected slope of "
-                    //      << " dist = " << dist << endl;
+                            << slope << " dist = " << dist );
 
-                    double e1 = average - p1.z();
+                    double e1 = average_elev_m - p1.z();
-                    double e2 = average - p2.z();
+                    double e2 = average_elev_m - p2.z();
                     // cout << "  z1 = " << p1.z() << "  z2 = " << p2.z() << endl;
                     // cout << "  e1 = " << e1 << "  e2 = " << e2 << endl;
 
@@ -320,11 +329,11 @@ TGAptSurface::TGAptSurface( const string& path,
 
                     slope_error = true;
 
-                    // cout << " (b) detected slope of " << slope 
+                    SG_LOG( SG_GENERAL, SG_DEBUG, " (b) detected slope of "
-                    //      << " dist = " << dist << endl;
+                            << slope << " dist = " << dist );
 
-                    double e1 = average - p1.z();
+                    double e1 = average_elev_m - p1.z();
-                    double e2 = average - p2.z();
+                    double e2 = average_elev_m - p2.z();
 
                     if ( fabs(e1) > fabs(e2) ) {
                         if ( slope > 0 ) {

src/Airports/GenAirports/apt_surface.hxx

@@ -57,12 +57,16 @@ private:
 
     Point3D min_deg, max_deg;
 
+    // externally seeded average airport elevation
+    double average_elev_m;
+
 public:
 
     // Constructor, specify min and max coordinates of desired area in
-    // lon/lat degrees
+    // lon/lat degrees, also please specify an "average" airport
+    // elevations in meters.
     TGAptSurface( const string &path, const string_list& elev_src,
-                  Point3D _min_deg, Point3D _max_deg );
+                  Point3D _min_deg, Point3D _max_deg, double _average_elev_m );
 
     // Destructor
     ~TGAptSurface();

src/Airports/GenAirports/build.cxx

@@ -58,6 +58,7 @@
 
 #include "apt_surface.hxx"
 #include "convex_hull.hxx"
+#include "elevations.hxx"
 #include "lights.hxx"
 #include "point2d.hxx"
 #include "poly_extra.hxx"
@@ -961,6 +962,18 @@ void build_airport( string airport_id, float alt_m,
 	tris_tc.push_back( base_tc );
     }
 
+    // Now that we have assembled all the airport geometry nodes into
+    // a list, calculate an "average" airport elevation based on all
+    // the actual airport node points.  This is more useful than
+    // calculating an average over the entire airport surface because
+    // it avoids biases introduced from the surrounding area if the
+    // airport is located in a bowl or on a hill.
+
+    double average = tgAverageElevation( root, elev_src,
+                                         nodes.get_node_list() );
+
+    // Now build the airport nurbs surface ...
+
     // calculation min/max coordinates of airport area
     Point3D min_deg(9999.0, 9999.0, 0), max_deg(-9999.0, -9999.0, 0);
     for ( j = 0; j < (int)nodes.get_node_list().size(); ++j ) {
@@ -1005,16 +1018,17 @@ void build_airport( string airport_id, float alt_m,
     // Extend the area a bit so we don't have wierd things on the edges
     double dlon = max_deg.lon() - min_deg.lon();
     double dlat = max_deg.lat() - min_deg.lat();
-    min_deg.setlon( min_deg.lon() - 0.2 * dlon );
+    min_deg.setlon( min_deg.lon() - 0.3 * dlon );
-    max_deg.setlon( max_deg.lon() + 0.2 * dlon );
+    max_deg.setlon( max_deg.lon() + 0.3 * dlon );
-    min_deg.setlat( min_deg.lat() - 0.2 * dlat );
+    min_deg.setlat( min_deg.lat() - 0.3 * dlat );
-    max_deg.setlat( max_deg.lat() + 0.2 * dlat );
+    max_deg.setlat( max_deg.lat() + 0.3 * dlat );
 
-    TGAptSurface apt_surf( root, elev_src, min_deg, max_deg );
+    TGAptSurface apt_surf( root, elev_src, min_deg, max_deg, average );
     SG_LOG(SG_GENERAL, SG_DEBUG, "Surface created");
 
     // calculate node elevations
-    point_list geod_nodes = calc_elevations( apt_surf, nodes.get_node_list(),
+    point_list geod_nodes = calc_elevations( apt_surf,
+                                             nodes.get_node_list(),
                                              0.0 );
     divided_base = calc_elevations( apt_surf, divided_base, 0.0 );
     SG_LOG(SG_GENERAL, SG_DEBUG, "DIVIDED");

src/Airports/GenAirports/elevations.cxx

@@ -0,0 +1,297 @@
+// elevations.cxx -- routines to help calculate DEM elevations for a
+//                   set of points
+//
+// Written by Curtis Olson, started April 2004.
+//
+// Copyright (C) 2004  Curtis L. Olson  - curt@flightgear.org
+//
+// This program is free software; you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation; either version 2 of the License, or
+// (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
+//
+// $Id$
+//
+
+
+#include <nurbs++/nurbsS.h>
+#include <nurbs++/nurbsSub.h>
+
+#include <simgear/constants.h>
+#include <simgear/math/sg_geodesy.hxx>
+#include <simgear/math/sg_types.hxx>
+#include <simgear/debug/logstream.hxx>
+
+#include <Array/array.hxx>
+
+#include "global.hxx"
+#include "apt_surface.hxx"
+
+SG_USING_NAMESPACE( PLib );
+
+
+// lookup node elevations for each point in the point_list.  Returns
+// average of all points.  Doesn't modify the original list.
+
+double tgAverageElevation( const string &root, const string_list elev_src,
+                           const point_list points_source )
+{
+    bool done = false;
+    unsigned int i;
+    TGArray array;
+
+    // make a copy so our routine is non-destructive.
+    point_list points = points_source;
+
+    // just bail if no work to do
+    if ( points.size() == 0 ) {
+        return 0.0;
+    }
+
+    // set all elevations to -9999
+    for ( i = 0; i < points.size(); ++i ) {
+        points[i].setz( -9999.0 );
+    }
+
+    while ( !done ) {
+	// find first node with -9999 elevation
+        Point3D first;
+        bool found_one = false;
+        for ( i = 0; i < points.size(); ++i ) {
+            if ( points[i].z() < -9000.0 && !found_one ) {
+                first = points[i];
+                found_one = true;
+            }
+        }
+
+	if ( found_one ) {
+	    SGBucket b( first.x(), first.y() );
+	    string base = b.gen_base_path();
+
+	    // try the various elevation sources
+            i = 0;
+            bool found_file = false;
+            while ( !found_file && i < elev_src.size() ) {
+                string array_path = root + "/" + elev_src[i] + "/" + base 
+                    + "/" + b.gen_index_str();
+                if ( array.open(array_path) ) {
+                    found_file = true;
+                    SG_LOG( SG_GENERAL, SG_DEBUG, "Using array_path = "
+                            << array_path );
+                }                    
+                i++;
+            }
+            
+            // this will fill in a zero structure if no array data
+            // found/opened
+	    array.parse( b );
+
+            // this will do a hasty job of removing voids by inserting
+            // data from the nearest neighbor (sort of)
+            array.remove_voids();
+
+	    // update all the non-updated elevations that are inside
+	    // this array file
+	    double elev;
+	    done = true;
+            for ( i = 0; i < points.size(); ++i ) {
+                if ( points[i].z() < -9000.0 ) {
+                    done = false;
+                    elev = array.altitude_from_grid( points[i].lon() * 3600.0,
+                                                     points[i].lat() * 3600.0 );
+                    if ( elev > -9000 ) {
+                        points[i].setz( elev );
+                        // cout << "interpolating for " << p << endl;
+                        // cout << p.x() << " " << p.y() << " " << p.z()
+                        //      << endl;
+                    }
+                }
+            }
+
+	    array.close();
+	} else {
+	    done = true;
+	}
+    }
+
+    // now find the average height of the queried points
+    double total = 0.0;
+    int count = 0;
+    for ( i = 0; i < points.size(); ++i ) {
+        total += points[i].z();
+        count++;
+    }
+    double average = total / (double) count;
+    SG_LOG(SG_GENERAL, SG_DEBUG, "Average surface height of point list = "
+           << average);
+
+    return average;
+}
+
+
+// lookup node elevations for each point in the specified nurbs++
+// matrix.  Returns average of all points.
+
+double tgCalcElevations( const string &root, const string_list elev_src,
+                         Matrix_Point3Dd &Pts ) {
+    bool done = false;
+    int i, j;
+    TGArray array;
+
+    // just bail if no work to do
+    if ( Pts.rows() == 0 || Pts.cols() == 0 ) {
+        return 0.0;
+    }
+
+    // set all elevations to -9999
+    for ( i = 0; i < Pts.cols(); ++i ) {
+        for ( j = 0; j < Pts.rows(); ++j ) {
+            Point3Dd p = Pts(j,i);
+            p.z() = -9999.0;
+        }
+    }
+
+    while ( !done ) {
+	// find first node with -9999 elevation
+        Point3Dd first;
+        bool found_one = false;
+        for ( i = 0; i < Pts.cols(); ++i ) {
+            for ( j = 0; j < Pts.rows(); ++j ) {
+                Point3Dd p = Pts(j,i);
+                if ( p.z() < -9000.0 && !found_one ) {
+                    first = p;
+                    found_one = true;
+                }
+            }
+        }
+
+	if ( found_one ) {
+	    SGBucket b( first.x(), first.y() );
+	    string base = b.gen_base_path();
+
+	    // try the various elevation sources
+            j = 0;
+            bool found_file = false;
+            while ( !found_file && j < (int)elev_src.size() ) {
+                string array_path = root + "/" + elev_src[j] + "/" + base 
+                    + "/" + b.gen_index_str();
+                if ( array.open(array_path) ) {
+                    found_file = true;
+                    SG_LOG( SG_GENERAL, SG_DEBUG, "Using array_path = "
+                            << array_path );
+                }                    
+                j++;
+            }
+            
+            // this will fill in a zero structure if no array data
+            // found/opened
+	    array.parse( b );
+
+            // this will do a hasty job of removing voids by inserting
+            // data from the nearest neighbor (sort of)
+            array.remove_voids();
+
+	    // update all the non-updated elevations that are inside
+	    // this array file
+	    double elev;
+	    done = true;
+            for ( i = 0; i < Pts.cols(); ++i ) {
+                for ( j = 0; j < Pts.rows(); ++j ) {
+                    Point3Dd p = Pts(j,i);
+                    if ( p.z() < -9000.0 ) {
+                        done = false;
+                        elev = array.altitude_from_grid( p.x() * 3600.0,
+                                                         p.y() * 3600.0 );
+                        if ( elev > -9000 ) {
+                            p.z() = elev;
+                            Pts(j,i) = p;
+                            // cout << "interpolating for " << p << endl;
+                            // cout << p.x() << " " << p.y() << " " << p.z()
+                            //      << endl;
+                        }
+                    }
+                }
+            }
+
+	    array.close();
+	} else {
+	    done = true;
+	}
+    }
+
+    // do some post processing for sanity's sake
+
+    // find the average height of the queried points
+    double total = 0.0;
+    int count = 0;
+    for ( i = 0; i < Pts.cols(); ++i ) {
+        for ( j = 0; j < Pts.rows(); ++j ) {
+            Point3Dd p = Pts(j,i);
+            total += p.z();
+            count++;
+        }
+    }
+    double average = total / (double) count;
+    SG_LOG(SG_GENERAL, SG_DEBUG, "Average surface height of nurbs matrix = "
+           << average);
+
+    // now go through the elevations and clamp them all to within
+    // +/-50m (164') of the average.
+    for ( i = 0; i < Pts.cols(); ++i ) {
+        for ( j = 0; j < Pts.rows(); ++j ) {
+            Point3Dd p = Pts(j,i);
+            if ( p.z() < average - max_clamp ) {
+                SG_LOG(SG_GENERAL, SG_DEBUG, "   clamping " << p.z()
+                       << " to " << average - max_clamp );
+                p.z() = average - max_clamp;
+                Pts(j,i) = p;
+            }
+            if ( p.z() > average + max_clamp ) {
+                SG_LOG(SG_GENERAL, SG_DEBUG, "   clamping " << p.z()
+                       << " to " << average + max_clamp );
+                p.z() = average + max_clamp;
+                Pts(j,i) = p;
+            }
+        }
+    }
+
+    return average;
+}
+
+
+// clamp all elevations to the specified range
+
+void tgClampElevations( Matrix_Point3Dd &Pts,
+                        double center_m, double max_clamp_m )
+{
+    int i, j;
+
+    // go through the elevations and clamp all elevations to within
+    // +/-max_m of the center_m elevation.
+    for ( i = 0; i < Pts.cols(); ++i ) {
+        for ( j = 0; j < Pts.rows(); ++j ) {
+            Point3Dd p = Pts(j,i);
+            if ( p.z() < center_m - max_clamp_m ) {
+                SG_LOG(SG_GENERAL, SG_DEBUG, "   clamping " << p.z()
+                       << " to " << center_m - max_clamp_m );
+                p.z() = center_m - max_clamp_m;
+                Pts(j,i) = p;
+            }
+            if ( p.z() > center_m + max_clamp_m ) {
+                SG_LOG(SG_GENERAL, SG_DEBUG, "   clamping " << p.z()
+                       << " to " << center_m + max_clamp_m );
+                p.z() = center_m + max_clamp_m;
+                Pts(j,i) = p;
+            }
+        }
+    }
+}

src/Airports/GenAirports/elevations.hxx

@@ -0,0 +1,54 @@
+// elevations.hxx -- routines to help calculate DEM elevations for a
+//                   set of points
+//
+// Written by Curtis Olson, started April 2004.
+//
+// Copyright (C) 2004  Curtis L. Olson  - curt@flightgear.org
+//
+// This program is free software; you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation; either version 2 of the License, or
+// (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
+//
+// $Id$
+//
+
+
+#include <nurbs++/nurbsS.h>
+#include <nurbs++/nurbsSub.h>
+
+#include <simgear/constants.h>
+#include <simgear/math/sg_geodesy.hxx>
+#include <simgear/math/sg_types.hxx>
+#include <simgear/debug/logstream.hxx>
+
+#include <Array/array.hxx>
+
+#include "global.hxx"
+#include "apt_surface.hxx"
+
+SG_USING_NAMESPACE( PLib );
+
+
+// lookup node elevations for each point in the point_list.  Returns
+// average of all points.  Doesn't modify the original list.
+double tgAverageElevation( const string &root, const string_list elev_src,
+                           const point_list points_source );
+
+// lookup node elevations for each point in the specified nurbs++
+// matrix.  Returns average of all points.
+double tgCalcElevations( const string &root, const string_list elev_src,
+                         Matrix_Point3Dd &Pts );
+
+// clamp all elevations to the specified range
+void tgClampElevations( Matrix_Point3Dd &Pts,
+                        double center_m, double max_clamp_m );

src/Airports/GenAirports/main.cxx

@@ -118,7 +118,7 @@ int main( int argc, char **argv ) {
     string_list elev_src;
     elev_src.clear();
 
-    sglog().setLogLevels( SG_GENERAL, SG_DEBUG );
+    sglog().setLogLevels( SG_GENERAL, SG_INFO );
 
     // parse arguments
     string work_dir = "";