diff --git a/src/Airports/GenAirports/Makefile.am b/src/Airports/GenAirports/Makefile.am index c0f7cdc1..3930364a 100644 --- a/src/Airports/GenAirports/Makefile.am +++ b/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 \ diff --git a/src/Airports/GenAirports/apt_surface.cxx b/src/Airports/GenAirports/apt_surface.cxx index 269ec7f5..316d0d69 100644 --- a/src/Airports/GenAirports/apt_surface.cxx +++ b/src/Airports/GenAirports/apt_surface.cxx @@ -33,12 +33,14 @@ #include +#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 - double average = calc_elevations( path, elev_src, dPts ); + // Lookup the elevations of all the grid points + 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 - // << " dist = " << dist << endl; + SG_LOG( SG_GENERAL, SG_DEBUG, " (a) detected slope of " + << slope << " dist = " << dist ); - double e1 = average - p1.z(); - double e2 = average - p2.z(); + double e1 = average_elev_m - p1.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 - // << " dist = " << dist << endl; + SG_LOG( SG_GENERAL, SG_DEBUG, " (b) detected slope of " + << slope << " dist = " << dist ); - double e1 = average - p1.z(); - double e2 = average - p2.z(); + double e1 = average_elev_m - p1.z(); + double e2 = average_elev_m - p2.z(); if ( fabs(e1) > fabs(e2) ) { if ( slope > 0 ) { diff --git a/src/Airports/GenAirports/apt_surface.hxx b/src/Airports/GenAirports/apt_surface.hxx index dffd66fb..046b961c 100644 --- a/src/Airports/GenAirports/apt_surface.hxx +++ b/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(); diff --git a/src/Airports/GenAirports/build.cxx b/src/Airports/GenAirports/build.cxx index 31880929..774fe212 100644 --- a/src/Airports/GenAirports/build.cxx +++ b/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 ); - max_deg.setlon( max_deg.lon() + 0.2 * dlon ); - min_deg.setlat( min_deg.lat() - 0.2 * dlat ); - max_deg.setlat( max_deg.lat() + 0.2 * dlat ); + min_deg.setlon( min_deg.lon() - 0.3 * dlon ); + max_deg.setlon( max_deg.lon() + 0.3 * dlon ); + min_deg.setlat( min_deg.lat() - 0.3 * 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"); diff --git a/src/Airports/GenAirports/elevations.cxx b/src/Airports/GenAirports/elevations.cxx new file mode 100644 index 00000000..36990dfa --- /dev/null +++ b/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 +#include + +#include +#include +#include +#include + +#include + +#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; + } + } + } +} diff --git a/src/Airports/GenAirports/elevations.hxx b/src/Airports/GenAirports/elevations.hxx new file mode 100644 index 00000000..3b85d124 --- /dev/null +++ b/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 +#include + +#include +#include +#include +#include + +#include + +#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 ); diff --git a/src/Airports/GenAirports/main.cxx b/src/Airports/GenAirports/main.cxx index af92784f..bd0e1110 100644 --- a/src/Airports/GenAirports/main.cxx +++ b/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 = "";