diff --git a/src/Lib/Array/array.cxx b/src/Lib/Array/array.cxx index 822be173..21b8b730 100644 --- a/src/Lib/Array/array.cxx +++ b/src/Lib/Array/array.cxx @@ -80,33 +80,29 @@ TGArray::TGArray( const string &file ) { // open an Array file (and fitted file if it exists) bool TGArray::open( const string& file_base ) { - // open input file (or read from stdin) - if ( file_base == "-" ) { - cout << " Opening array data pipe from stdin" << endl; - // fd = stdin; - // fd = gzdopen(STDIN_FILENO, "r"); - cout << " Not yet ported ..." << endl; - return false; + bool success = true; + + // open array data file + string array_name = file_base + ".arr.gz"; + array_in = new sg_gzifstream( array_name ); + if ( ! array_in->is_open() ) { + cout << " Cannot open " << array_name << endl; + success = false; } else { - // open array data file - string array_name = file_base + ".arr.gz"; - array_in = new sg_gzifstream( array_name ); - // open fitted data file - string fitted_name = file_base + ".fit.gz"; - fitted_in = new sg_gzifstream( fitted_name ); - if ( ! array_in->is_open() ) { - cout << " Cannot open " << array_name << endl; - return false; - } - cout << " Opening array data file: " << array_name << endl; - if ( ! fitted_in->is_open() ) { - cout << " Cannot open " << fitted_name << endl; - return false; - } - cout << " Opening fitted data file: " << fitted_name << endl; + cout << " Opening array data file: " << array_name << endl; } - return true; + // open fitted data file + string fitted_name = file_base + ".fit.gz"; + fitted_in = new sg_gzifstream( fitted_name ); + if ( ! fitted_in->is_open() ) { + cout << " Cannot open " << fitted_name << endl; + success = false; + } else { + cout << " Opening fitted data file: " << fitted_name << endl; + } + + return success; } @@ -146,6 +142,7 @@ TGArray::parse( SGBucket& b ) { cout << " Done parsing\n"; +#if 0 // need two passes to ensure that all voids are removed (unless entire // array is a void.) bool have_void = true; @@ -221,6 +218,7 @@ TGArray::parse( SGBucket& b ) { } } } +#endif } else { // file not open (not found?), fill with zero'd data @@ -250,7 +248,7 @@ TGArray::parse( SGBucket& b ) { } // Parse/load the array data file - if ( fitted_in && fitted_in->is_open() ) { + if ( fitted_in->is_open() ) { fit_on_the_fly = false; int fitted_size; double x, y, z, error; @@ -272,6 +270,8 @@ TGArray::parse( SGBucket& b ) { } + + // add a node to the output corner node list void TGArray::add_corner_node( int i, int j, double val ) { @@ -456,9 +456,37 @@ int TGArray::fit( double error ) { } +// Return the elevation of the closest non-void grid point to lon, lat +double TGArray::closest_nonvoid_elev( double lon, double lat ) const { + double mindist = 99999999999.9; + double minelev = -9999.0; + Point3D p0( lon, lat, 0.0 ); + + for ( int row = 0; row < rows; row++ ) { + for ( int col = 0; col < cols; col++ ) { + Point3D p1(originx + col * col_step, originy + row * row_step, 0.0); + double dist = p0.distance3D( p1 ); + double elev = in_data[col][row]; + if ( dist < mindist && elev > -9000 ) { + mindist = dist; + minelev = elev; + cout << "dist = " << mindist; + cout << " elev = " << elev << endl; + } + } + } + + if ( minelev > -9999.0 ) { + return minelev; + } else { + return 0.0; + } +} + + // return the current altitude based on grid data. We should rewrite // this to interpolate exact values, but for now this is good enough -double TGArray::interpolate_altitude( double lon, double lat ) const { +double TGArray::altitude_from_grid( double lon, double lat ) const { // we expect incoming (lon,lat) to be in arcsec for now double xlocal, ylocal, dx, dy, zA, zB, elev; @@ -524,6 +552,11 @@ double TGArray::interpolate_altitude( double lon, double lat ) const { // printf(" (x2,y2,z2) = (%d,%d,%d)\n", x2, y2, z2); // printf(" (x3,y3,z3) = (%d,%d,%d)\n", x3, y3, z3); + if ( z1 < -9000 || z2 < -9000 || z3 < -9000 ) { + // don't interpolate off a void + return closest_nonvoid_elev( lon, lat ); + } + zA = dx * (z2 - z1) + z1; zB = dx * (z3 - z1) + z1; @@ -555,6 +588,11 @@ double TGArray::interpolate_altitude( double lon, double lat ) const { // printf(" (x2,y2,z2) = (%d,%d,%d)\n", x2, y2, z2); // printf(" (x3,y3,z3) = (%d,%d,%d)\n", x3, y3, z3); + if ( z1 < -9000 || z2 < -9000 || z3 < -9000 ) { + // don't interpolate off a void + return closest_nonvoid_elev( lon, lat ); + } + zA = dy * (z2 - z1) + z1; zB = dy * (z3 - z1) + z1; @@ -568,7 +606,7 @@ double TGArray::interpolate_altitude( double lon, double lat ) const { } } - return(elev); + return elev; } diff --git a/src/Lib/Array/array.hxx b/src/Lib/Array/array.hxx index ad0c9b6d..0f28780a 100644 --- a/src/Lib/Array/array.hxx +++ b/src/Lib/Array/array.hxx @@ -111,10 +111,13 @@ public: // add a node to the output fitted node list void add_fit_node( int i, int j, double val ); + // Return the elevation of the closest non-void grid point to lon, lat + double closest_nonvoid_elev( double lon, double lat ) const; + // return the current altitude based on grid data. We should // rewrite this to interpolate exact values, but for now this is // good enough - double interpolate_altitude( double lon, double lat ) const; + double altitude_from_grid( double lon, double lat ) const; // Informational methods inline double get_originx() const { return originx; } diff --git a/src/Lib/Array/testarray.cxx b/src/Lib/Array/testarray.cxx index 973f5fc4..13e34d15 100644 --- a/src/Lib/Array/testarray.cxx +++ b/src/Lib/Array/testarray.cxx @@ -30,7 +30,7 @@ int main( int argc, char **argv ) { lon *= 3600; lat *= 3600; - cout << " " << a.interpolate_altitude(lon, lat) << endl; + cout << " " << a.altitude_from_grid(lon, lat) << endl; a.fit( 100 ); diff --git a/src/Lib/Array/testgts.cxx b/src/Lib/Array/testgts.cxx index c7ef92ae..92172f6d 100644 --- a/src/Lib/Array/testgts.cxx +++ b/src/Lib/Array/testgts.cxx @@ -131,25 +131,25 @@ int main( int argc, char **argv ) { // Make the corner vertices (enclosing exactly the DEM coverage area) x = basex; y = basey; - z = a.interpolate_altitude( x, y ); + z = a.altitude_from_grid( x, y ); cout << "adding = " << Point3D( x, y, z) << endl; GtsVertex *v1 = gts_vertex_new( gts_vertex_class(), x, y, z ); x = basex + dx * (a.get_cols() - 1); y = basey; - z = a.interpolate_altitude( x, y ); + z = a.altitude_from_grid( x, y ); cout << "adding = " << Point3D( x, y, z) << endl; GtsVertex *v2 = gts_vertex_new( gts_vertex_class(), x, y, z ); x = basex + dx * (a.get_cols() - 1); y = basey + dy * (a.get_rows() - 1); - z = a.interpolate_altitude( x, y ); + z = a.altitude_from_grid( x, y ); cout << "adding = " << Point3D( x, y, z) << endl; GtsVertex *v3 = gts_vertex_new( gts_vertex_class(), x, y, z ); x = basex; y = basey + dy * (a.get_rows() - 1); - z = a.interpolate_altitude( x, y ); + z = a.altitude_from_grid( x, y ); cout << "adding = " << Point3D( x, y, z) << endl; GtsVertex *v4 = gts_vertex_new( gts_vertex_class(), x, y, z );