From 1d71909929813d4a197cab27423a6309d8698334 Mon Sep 17 00:00:00 2001 From: curt Date: Wed, 26 Feb 2003 19:30:25 +0000 Subject: [PATCH] Initial revision of code to read SRTM "hgt" format data. --- src/Lib/HGT/Makefile.am | 14 +++ src/Lib/HGT/hgt.cxx | 252 ++++++++++++++++++++++++++++++++++++++++ src/Lib/HGT/hgt.hxx | 104 +++++++++++++++++ src/Lib/HGT/testhgt.cxx | 17 +++ 4 files changed, 387 insertions(+) create mode 100644 src/Lib/HGT/Makefile.am create mode 100644 src/Lib/HGT/hgt.cxx create mode 100644 src/Lib/HGT/hgt.hxx create mode 100644 src/Lib/HGT/testhgt.cxx diff --git a/src/Lib/HGT/Makefile.am b/src/Lib/HGT/Makefile.am new file mode 100644 index 00000000..451232a2 --- /dev/null +++ b/src/Lib/HGT/Makefile.am @@ -0,0 +1,14 @@ +noinst_LIBRARIES = libHGT.a + +libHGT_a_SOURCES = hgt.cxx hgt.hxx + +noinst_PROGRAMS = testhgt + +testhgt_SOURCES = testhgt.cxx + +testhgt_LDADD = \ + $(top_builddir)/src/Lib/HGT/libHGT.a \ + -lsgbucket -lsgmisc -lz + +INCLUDES = -I$(top_srcdir)/src + diff --git a/src/Lib/HGT/hgt.cxx b/src/Lib/HGT/hgt.cxx new file mode 100644 index 00000000..39a31fec --- /dev/null +++ b/src/Lib/HGT/hgt.cxx @@ -0,0 +1,252 @@ +// hgt.hxx -- SRTM "hgt" data management class +// +// Written by Curtis Olson, started February 2003. +// +// Copyright (C) 2003 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$ + + +#ifdef HAVE_CONFIG_H +# include +#endif + +#include + +#include // atof() + +#ifdef SG_HAVE_STD_INCLUDES +# include +#else +# include +#endif + +#include +#include + +#ifdef _MSC_VER +# include +#endif + +#include "hgt.hxx" + +SG_USING_STD(cout); +SG_USING_STD(endl); + + +TGHgt::TGHgt( int _res ) { + hgt_resolution = _res; + + data = new short int[MAX_HGT_SIZE][MAX_HGT_SIZE]; + output_data = new short int[MAX_HGT_SIZE][MAX_HGT_SIZE]; +} + + +TGHgt::TGHgt( int _res, const SGPath &file ) { + hgt_resolution = _res; + + data = new short int[MAX_HGT_SIZE][MAX_HGT_SIZE]; + output_data = new short int[MAX_HGT_SIZE][MAX_HGT_SIZE]; + + TGHgt::open( file ); +} + + +// open an HGT file +bool +TGHgt::open ( const SGPath &f ) { + SGPath file_name = f; + + // open input file (or read from stdin) + if ( file_name.str() == "-" ) { + cout << "Loading HGT data file: stdin" << endl; + if ( (fd = gzdopen(STDIN_FILENO, "r")) == NULL ) { + cout << "ERROR: opening stdin" << endl; + return false; + } + } else { + cout << "Loading HGT data file: " << file_name.str() << endl; + if ( (fd = gzopen( file_name.c_str(), "rb" )) == NULL ) { + SGPath file_name_gz = file_name; + file_name_gz.append( ".gz" ); + if ( (fd = gzopen( file_name_gz.c_str(), "rb" )) == NULL ) { + cout << "ERROR: opening " << file_name.str() << " or " + << file_name_gz.str() << "for reading!" << endl; + return false; + } + } + } + + // Determine originx/originy from file name + string name = file_name.file(); + cout << " Name = " << name << endl; + originy = atof( name.substr(1, 2).c_str() ) * 3600.0; + if ( name.substr(0, 1) == "S" ) { + originy = -originy; + } + originx = atof( name.substr(4, 3).c_str() ) * 3600.0; + if ( name.substr(3, 1) == "W" ) { + originx = -originx; + } + cout << " Origin = " << originx << ", " << originy << endl; + + return true; +} + + +// close an HGT file +bool +TGHgt::close () { + gzclose(fd); + return true; +} + + +// load an hgt file +bool +TGHgt::load( ) { + int size; + if ( hgt_resolution == 1 ) { + cols = rows = size = 3601; + col_step = row_step = 1; + } else if ( hgt_resolution == 3 ) { + cols = rows = size = 1201; + col_step = row_step = 3; + } else { + cout << "Unknown HGT resolution, only 1 and 3 arcsec formats" << endl; + cout << " are supported!" << endl; + return false; + } + + short int *var; + for ( int row = 0; row < size; ++row ) { + for ( int col = 0; col < size; ++col ) { + var = &data[row][col]; + if ( gzread ( fd, var, sizeof(short) ) != sizeof(short) ) { + return false; + } + if ( sgIsLittleEndian() ) { + sgEndianSwap( (unsigned short int*)var); + } + } + } + + return true; +} + + +// write out the area of data covered by the specified bucket. Data +// is written out column by column starting at the lower left hand +// corner. +bool +TGHgt::write_area( const string& root, SGBucket& b ) { + // calculate some boundaries + double min_x = ( b.get_center_lon() - 0.5 * b.get_width() ) * 3600.0; + double max_x = ( b.get_center_lon() + 0.5 * b.get_width() ) * 3600.0; + + double min_y = ( b.get_center_lat() - 0.5 * b.get_height() ) * 3600.0; + double max_y = ( b.get_center_lat() + 0.5 * b.get_height() ) * 3600.0; + + cout << b << endl; + cout << "width = " << b.get_width() << " height = " << b.get_height() + << endl; + cout << "min = " << min_x << "," << min_y + << " max = " << max_x << "," << max_y << endl; + int start_x = (int)((min_x - originx) / col_step); + int span_x = (int)(b.get_width() * 3600.0 / col_step); + + int start_y = (int)((min_y - originy) / row_step); + int span_y = (int)(b.get_height() * 3600.0 / row_step); + + cout << "start_x = " << start_x << " span_x = " << span_x << endl; + cout << "start_y = " << start_y << " span_y = " << span_y << endl; + + // Do a simple sanity checking. But, please, please be nice to + // this write_area() routine and feed it buckets that coincide + // well with the underlying grid structure and spacing. + + if ( ( min_x < originx ) + || ( max_x > originx + cols * col_step ) + || ( min_y < originy ) + || ( max_y > originy + rows * row_step ) ) { + cout << " ERROR: bucket at least partially outside HGT data range!" << + endl; + return false; + } + + // If the area is all ocean, skip it. + if ( !has_non_zero_elev(start_x, span_x, start_y, span_y) ) { + cout << "Tile is all zero elevation: skipping" << endl; + return false; + } + + // generate output file name + string base = b.gen_base_path(); + string path = root + "/" + base; +#ifdef _MSC_VER + fg_mkdir( path.c_str() ); +#else + string command = "mkdir -p " + path; + system( command.c_str() ); +#endif + + string array_file = path + "/" + b.gen_index_str() + ".arr.gz"; + cout << "array_file = " << array_file << endl; + + // write the file + gzFile fp; + if ( (fp = gzopen( array_file.c_str(), "wb9" )) == NULL ) { + cout << "ERROR: cannot open " << array_file << " for writing!" << endl; + exit(-1); + } + + gzprintf( fp, "%d %d\n", (int)min_x, (int)min_y ); + gzprintf( fp, "%d %d %d %d\n", span_x + 1, (int)col_step, + span_y + 1, (int)row_step ); + for ( int i = start_x; i <= start_x + span_x; ++i ) { + for ( int j = start_y; j <= start_y + span_y; ++j ) { + gzprintf( fp, "%d ", (int)data[i][j] ); + } + gzprintf( fp, "\n" ); + } + gzclose(fp); + + return true; +} + + +TGHgt::~TGHgt() { + // printf("class TGHgt DEstructor called.\n"); + delete [] data; + delete [] output_data; +} + + +bool +TGHgt::has_non_zero_elev (int start_x, int span_x, + int start_y, int span_y) const +{ + for ( int i = start_x; i < start_x + span_x; i++ ) { + for ( int j = start_y; j < start_y + span_y; j++ ) { + if ( data[i][j] != 0 ) + return true; + } + } + + return false; +} + diff --git a/src/Lib/HGT/hgt.hxx b/src/Lib/HGT/hgt.hxx new file mode 100644 index 00000000..6e316dbe --- /dev/null +++ b/src/Lib/HGT/hgt.hxx @@ -0,0 +1,104 @@ +// hgt.hxx -- SRTM "hgt" data management class +// +// Written by Curtis Olson, started February 2003. +// +// Copyright (C) 2003 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$ + + +#ifndef _HGT_HXX +#define _HGT_HXX + +#ifdef HAVE_CONFIG_H +# include +#endif + +#include + +#include +#include + + +#define MAX_HGT_SIZE 3601 + + +class TGHgt { + +private: + + // file pointer for input + gzFile fd; + + // coordinates (in arc seconds) of south west corner + double originx, originy; + + // number of columns and rows + int cols, rows; + + // Distance between column and row data points (in arc seconds) + double col_step, row_step; + + // pointers to the actual grid data allocated here + short int (*data)[MAX_HGT_SIZE]; + short int (*output_data)[MAX_HGT_SIZE]; + + int hgt_resolution; + +public: + + // Constructor, _res must be either "1" for the 1arcsec data or + // "3" for the 3arcsec data. + TGHgt( int _res ); + TGHgt( int _res, const SGPath &file ); + + // Destructor + ~TGHgt(); + + // open an HGT file (use "-" if input is coming from stdin) + bool open ( const SGPath &file ); + + // close an HGT file + bool close(); + + // load an hgt file + bool load(); + + // write out the area of data covered by the specified bucket. + // Data is written out column by column starting at the lower left + // hand corner. + bool write_area( const string& root, SGBucket& b ); + + // Informational methods + inline double get_originx() const { return originx; } + inline double get_originy() const { return originy; } + inline int get_cols() const { return cols; } + inline int get_rows() const { return rows; } + inline double get_col_step() const { return col_step; } + inline double get_row_step() const { return row_step; } + + /** + * Test whether an area contains any non-zero elevations. + */ + bool has_non_zero_elev (int start_x, int span_x, + int start_y, int span_y) const; +}; + + +#endif // _HGT_HXX + + diff --git a/src/Lib/HGT/testhgt.cxx b/src/Lib/HGT/testhgt.cxx new file mode 100644 index 00000000..4d81ba04 --- /dev/null +++ b/src/Lib/HGT/testhgt.cxx @@ -0,0 +1,17 @@ +#include + +#include "hgt.hxx" + +int main() { + SGPath path1( "/stage/fgfs03/SRTM/United_States/N33W080/1arcsec/N33W080.hgt.gz" ); + SGPath path3( "/stage/fgfs03/SRTM/United_States/N33W080/3arcsec/N33W080.hgt.gz" ); + + TGHgt hgt1( 1, path1 ); + TGHgt hgt3( 3, path3 ); + + hgt1.load(); + hgt3.load(); + + SGBucket b(-79.5, 33.5); + hgt1.write_area( ".", b ); +}