fly/terragear
1
0
Fork 0
Code Activity

Merge /u/jreh/flightgear/ branch cliff-updates into next

Browse source 
https://sourceforge.net/p/flightgear/terragear/merge-requests/15/
This commit is contained in:
xDraconian 2019-01-14 05:09:55 +00:00
commit 810a88b754
30 changed files with 1860 additions and 102 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

54
README.howto
View file

@ -1,5 +1,6 @@
Original version by Alexei Novikov <anovikov@heron.itep.ru> Sep. 3, 1999
Updates by Curtis Olson <http://www.flightgear.org/~curt>
Cliff description by James Hester 2018
Overview
@ -120,3 +121,56 @@ tg-construct-client in rude mode. After a while scenery is ready.
That's all folks,
Alexei.
Addendum: Adding cliffs
=======================
Cliffs represent discontinuities in the height array described above. This
is problematic because:
(1) Most height grids will give an average height in the area around the
grid
(2) Interpolation between grid points to get the local height will convert
a cliff into a much more gentle slope.
To add cliffs into the scenery using OpenStreetMap data, the following
approach can be used:
(1) Extract cliffs from the OSM datafile using your favourite
tool, for example for all cliffs near Sydney, Australia:
osmosis --read-pbf australia-latest.osm.pbf --bounding-box top=-33 left=150
bottom=-35 right=152 completeWays=yes --lp --tf accept-ways 'natural=cliff'
--tf reject-relations --used-node --write-xml file='my_cliffs.osm'
(2) Convert the XML above to Shapefiles:
ogr2ogr -where "OGR_GEOMETRY='LineString'" -f 'ESRI shapefile' -lco SHPT=ARC
cs_cliffs my_cliffs.osm
(3) Run cliff-decode to place cliff information into the height grid working
directory created previously (<data>/cs_cliffs in the example below):
cliff-decode <work_base>/SRTM-3 <data>/cs_cliffs
(4) Optionally run rectify_height to fix grid heights. This is recommended
as otherwise cliffs will have uneven edges and bases as the grid points
vary in distance from the cliff:
rectify_height --work-dir=<work_base> --height-dir=SRTM-3 --min-lon=151.0
--max-lon=152.0 --min-lat=-34.0 --max-lat=-33.0 --min-dist=100
The min-dist parameter sets the distance from the cliff beyond which points
are trustworthy. 100m is the default.
(5) Place the cliffs into the scenery. This is identical to the process for
rivers and roads and uses the same cliff information from step (2):
ogr-decode --line-width 5 --area-type Rock <work_base>/Cliffs <data>/cs_cliffs"
Note that an area-type of Cliffs will allow custom cliff texturing, if this
effect is available. This command places a 5m wide strip of rock at the
position of the cliffs. One side of this strip will have a calculated elevation
at the top of the cliff, and the other will be at the bottom.
(6) Run tg-construct as usual, including 'Rock' or 'Cliffs' polygon types to
make sure that cliffs appear.

1
src/Airports/GenAirports850/CMakeLists.txt
View file

@ -1,5 +1,6 @@
include_directories(${GDAL_INCLUDE_DIR})
include_directories(${PROJECT_SOURCE_DIR}/src/Lib)
include_directories(${PROJECT_SOURCE_DIR}/src/Lib/terragear)
add_executable(genapts850
airport.hxx airport.cxx

2
src/BuildTiles/CMakeLists.txt
View file

@ -1,5 +1,7 @@
include_directories(${PROJECT_SOURCE_DIR}/src/Lib)
include_directories(${PROJECT_SOURCE_DIR}/src/BuildTiles)
include_directories(${PROJECT_SOURCE_DIR}/src/Lib/Array)
include_directories(${PROJECT_SOURCE_DIR}/src/Lib/terragear)
#add_subdirectory(Parallel)
add_subdirectory(Main)

1
src/BuildTiles/Main/CMakeLists.txt
View file

@ -1,4 +1,5 @@
include_directories(${GDAL_INCLUDE_DIR})
include_directories(${PROJECT_SOURCE_DIR}/src/Lib/terragear)
add_executable(tg-construct
tgconstruct.hxx

3
src/BuildTiles/Main/default_priorities.txt
View file

@ -30,6 +30,7 @@ River stream
IntermittentStream stream
Watercourse stream
Canal stream
Cliffs cliff # A cliff face
Glacier other # Solid ice/snow
PackIce other # Water with ice packs
PolarIce other
@ -85,7 +86,9 @@ ShrubGrassCover other # Mixed Shrubland/Grassland
SavannaCover other # Savanna
Orchard other # Orchard
DeciduousForest other # Deciduous Forest
DeciduousBroadCover other # Deciduous Forest
EvergreenForest other # Evergreen Forest
EvergreenBroadCover other # Evergreen Forest
MixedForest other # Mixed Forest
RainForest other # Rain Forest
BarrenCover other # Barren or Sparsely Vegetated

10
src/BuildTiles/Main/priorities.hxx
View file

@ -140,6 +140,14 @@ public:
}
}
bool is_cliff_area( unsigned int p ) const {
if (area_defs[p].GetCategory() == "cliff" ) {
return true;
} else {
return false;
}
}
std::string const& get_area_name( unsigned int p ) const {
return area_defs[p].GetName();
}
@ -172,4 +180,4 @@ private:
unsigned int sliver_area_priority;
};
#endif // _PRIORITIES_HXX
#endif // _PRIORITIES_HXX

2
src/BuildTiles/Main/tgconstruct_elevation.cxx
View file

@ -222,4 +222,4 @@ void TGConstruct::CalcElevations( void )
}
}
}
}
}

4
src/BuildTiles/Main/tgconstruct_poly.cxx
View file

@ -65,8 +65,10 @@ int TGConstruct::LoadLandclassPolys( void ) {
}
string lext = p.complete_lower_extension();
string lastext = p.extension();
if ((lext == "arr") || (lext == "arr.gz") || (lext == "btg.gz") ||
(lext == "fit") || (lext == "fit.gz") || (lext == "ind"))
(lext == "fit") || (lext == "fit.gz") || (lext == "ind") ||
(lastext == "cliffs") || (lext == "arr.rectified.gz") || (lext == "arr.new.gz"))
{
// skipped!
} else {

15
src/Lib/Array/CMakeLists.txt
View file

@ -3,6 +3,9 @@
add_library(Array STATIC
array.cxx array.hxx
)
include_directories(${PROJECT_SOURCE_DIR}/src/Lib/terragear)
install( TARGETS Array
RUNTIME DESTINATION bin
LIBRARY DESTINATION lib
@ -12,13 +15,23 @@ install( FILES array.hxx DESTINATION include/tg )
add_executable(test_array testarray.cxx)
add_executable(rectify_height rectify_height.cxx)
target_link_libraries(test_array
Array
terragear
${ZLIB_LIBRARY}
${SIMGEAR_CORE_LIBRARIES}
${SIMGEAR_CORE_LIBRARY_DEPENDENCIES})
install(TARGETS test_array RUNTIME DESTINATION bin)
target_link_libraries(rectify_height
Array
terragear
${ZLIB_LIBRARY}
${SIMGEAR_CORE_LIBRARIES}
${SIMGEAR_CORE_LIBRARY_DEPENDENCIES})
install(TARGETS test_array rectify_height RUNTIME DESTINATION bin)
if (MSVC)
set_target_properties( test_array PROPERTIES DEBUG_POSTFIX d )
endif ()

577
src/Lib/Array/array.cxx
View file

@ -24,6 +24,7 @@
#endif
#include <cstring>
#include <iomanip> //for setprecision
#include <simgear/compiler.h>
#include <simgear/io/iostreams/sgstream.hxx>
@ -56,15 +57,27 @@ TGArray::TGArray( const string &file ):
// open an Array file (and fitted file if it exists)
// Also open cliffs file if it exists. By default a
// rectified file is searched for, if it doesn't exist
// we load the unrectified version.
bool TGArray::open( const string& file_base ) {
// open array data file
string array_name = file_base + ".arr.gz";
string array_name = file_base + ".arr.rectified.gz";
rectified = true;
array_in = gzopen( array_name.c_str(), "rb" );
if (array_in == NULL) {
return false;
}
// try unrectified
array_name = file_base + ".arr.gz";
array_in = gzopen(array_name.c_str(), "rb");
if (array_in == NULL) {
return false;
} else {
rectified = false;
}
}
SG_LOG(SG_GENERAL,SG_DEBUG,"Loaded height array " << array_name);
// open fitted data file
string fitted_name = file_base + ".fit.gz";
fitted_in = new sg_gzifstream( fitted_name );
@ -79,7 +92,8 @@ bool TGArray::open( const string& file_base ) {
} else {
SG_LOG(SG_GENERAL, SG_DEBUG, " Opening fitted data file: " << fitted_name );
}
// open any cliffs data file
load_cliffs(file_base);
return (array_in != NULL) ? true : false;
}
@ -101,6 +115,43 @@ TGArray::close() {
return true;
}
//This code adapted from tgconstruct::LoadLandclassPolys
//All polys in the bucket should be contours which we load
//into our contour list.
void TGArray::load_cliffs(const string & height_base)
{
//Get the directory so we can list the children
tgPolygon poly; //actually a contour but whatever...
int total_contours_read = 0;
SGPath b(height_base);
simgear::Dir d(b.dir());
simgear::PathList files = d.children(simgear::Dir::TYPE_FILE);
for (const SGPath& p: files) {
if (p.file_base() != b.file_base()) {
continue;
}
string lext = p.lower_extension();
if (lext == "cliffs") {
gzFile fp = gzopen( p.c_str(), "rb" );
unsigned int count;
sgReadUInt( fp, &count );
SG_LOG( SG_GENERAL, SG_DEBUG, " Load " << count << " contours from " << p.realpath() );
for ( unsigned int i=0; i<count; i++ ) {
poly.LoadFromGzFile( fp );
if ( poly.Contours()==1 ) { //should always have one contour
cliffs_list.push_back(poly.GetContour(0));
} else {
SG_LOG( SG_GENERAL, SG_WARN, " Found " << poly.Contours() << " contours in " << p.realpath() );
}
}
}
}
}
void
TGArray::unload( void ) {
if (array_in) {
@ -166,7 +217,7 @@ TGArray::parse( SGBucket& b ) {
}
}
return true;
return true;
}
void TGArray::parse_bin()
@ -197,6 +248,40 @@ void TGArray::parse_bin()
sgReadShort(array_in, cols * rows, in_data);
}
// Write out an array. If rectified is true, the heights have been adjusted
// for discontinuities.
void TGArray::write_bin(const string root_dir, bool rectified, SGBucket& b) {
// generate output file name
string base = b.gen_base_path();
string path = root_dir + "/" + base;
string extension = ".arr.new.gz";
if (rectified) extension = ".arr.rectified.gz";
SGPath sgp( path );
sgp.append( "dummy" );
sgp.create_dir( 0755 );
string array_file = path + "/" + b.gen_index_str() + extension;
SG_LOG(SG_GENERAL, SG_DEBUG, "array_file = " << array_file );
// write the file
gzFile fp;
if ( (fp = gzopen( array_file.c_str(), "wb9" )) == NULL ) {
SG_LOG(SG_GENERAL, SG_ALERT, "ERROR: cannot open " << array_file << " for writing!" );
return;
}
int32_t header = 0x54474152; //'TGAR'
sgWriteLong(fp,header);
sgWriteInt(fp,originx);
sgWriteInt(fp,originy);
sgWriteInt(fp,cols);
sgWriteInt(fp,col_step);
sgWriteInt(fp,rows);
sgWriteInt(fp,row_step);
sgWriteShort(fp, rows*cols, in_data);
gzclose(fp);
}
// write an Array file
bool TGArray::write( const string root_dir, SGBucket& b ) {
// generate output file name
@ -230,7 +315,6 @@ bool TGArray::write( const string root_dir, SGBucket& b ) {
return true;
}
// do our best to remove voids by picking data from the nearest neighbor.
void TGArray::remove_voids( ) {
// need two passes to ensure that all voids are removed (unless entire
@ -316,14 +400,15 @@ void TGArray::remove_voids( ) {
// Return the elevation of the closest non-void grid point to lon, lat
// Lon, lat in arcsec
double TGArray::closest_nonvoid_elev( double lon, double lat ) const {
double mindist = 99999999999.9;
double minelev = -9999.0;
SGGeod p0 = SGGeod::fromDeg( lon, lat );
SGGeod p0 = SGGeod::fromDeg( lon/3600.0, lat/3600.0 );
for ( int row = 0; row < rows; row++ ) {
for ( int col = 0; col < cols; col++ ) {
SGGeod p1 = SGGeod::fromDeg( originx + col * col_step, originy + row * row_step );
SGGeod p1 = SGGeod::fromDeg( (originx + col * col_step)/3600.0, (originy + row * row_step)/3600.0 );
double dist = SGGeodesy::distanceM( p0, p1 );
double elev = get_array_elev(col, row);
if ( dist < mindist && elev > -9000 ) {
@ -340,7 +425,162 @@ double TGArray::closest_nonvoid_elev( double lon, double lat ) const {
}
}
//Find and remember all points that are bad because they are
//too close to a cliff
std::vector<int> TGArray::collect_bad_points(const double bad_zone) {
std::vector<int> bad_points; //local to avoid multi-thread issues
for( int horiz=0;horiz<cols;horiz++ ) {
double lon = (originx + col_step*horiz)/3600;
for( int vert=0;vert<rows;vert++ ) {
double lat = (originy + row_step*vert)/3600;
if( is_near_cliff(lon,lat,bad_zone) ) {
bad_points.push_back(horiz+vert*cols);
}
}
}
return bad_points;
}
// Check to see if the specified grid point is bad
bool TGArray::is_bad_point(const int xgrid, const int ygrid, const std::vector<int> bad_points) const {
int grididx;
grididx = xgrid+ygrid*cols;
auto result = std::find( std::begin(bad_points),std::end(bad_points),grididx );
if ( result != std::end(bad_points) ) return true;
return false;
}
//This may collide with other threads, but as they will both be writing
//the correct height, this is harmless.
void TGArray::rectify_heights( const double bad_zone ) {
double new_ht;
std::vector<int> rectified,bad_points;
int total_rectified;
bad_points = collect_bad_points( bad_zone );
do {
for ( auto pt : bad_points ) {
int ygrid = pt/cols;
int xgrid = pt - ygrid*cols;
new_ht = rectify_point( xgrid,ygrid,bad_points );
if (new_ht > -9999) {
rectified.push_back(pt);
set_array_elev( xgrid,ygrid,(int) new_ht );
}
}
total_rectified = rectified.size();
SG_LOG(SG_GENERAL, SG_DEBUG, "Rectified " << total_rectified << " points ");
if( total_rectified > 0 ) {
for( auto r : rectified ) {
bad_points.erase( std::remove( std::begin(bad_points), std::end(bad_points),r) );
}
rectified.clear();
}
} while ( total_rectified > 0 );
if( bad_points.size() > 0 ) {
SG_LOG(SG_GENERAL, SG_DEBUG, "Failed to rectify " << bad_points.size() << " points");
}
}
/* If we have cliffs, it is possible that a grid point will be too close
to the cliff. In this case, the SRTM-3 data appears to average the height
in the region of the point, which makes the height unreliable. This routine
searches for three neighbouring points that are reliable, and form a rectangle
with the target point, and calculates the height from the plane passing
through the three known points.
* * *
* x *
* * *
TODO: Handle points on the boundaries. */
double TGArray::rectify_point(const int xgrid, const int ygrid, const std::vector<int> bad_points) const {
//xgrid: grid units horizontally
//ygrid: grid units vertically
//Loop over corner points, if no points available, give up
int corners[4][2]; //possible corners
int final_pts[3][2]; // rectangle corners
int pt_cnt = 0;
double centre_long, centre_lat;
double cliff_error = col_step; //Assume row step, col step the same
int original_height = get_array_elev(xgrid,ygrid);
centre_long = (originx + col_step*xgrid)/3600;
centre_lat = (originy + row_step*ygrid)/3600;
for ( int horiz = -1; horiz <= 1; horiz+=2 ) {
if (xgrid + horiz >= cols || xgrid + horiz < 0) continue; //edge of bucket
double test_long = centre_long + (col_step*horiz)/3600;
for ( int vert = -1; vert <= 1; vert+=2 ) {
if (ygrid + vert >= rows || ygrid + vert < 0) continue; //edge of bucket
double test_lat = centre_lat + (row_step*vert)/3600;
if ( !is_bad_point( xgrid+horiz,ygrid+vert,bad_points ) && //can trust height
check_points( test_long,test_lat,centre_long,centre_lat ) ) { //same side
corners[pt_cnt][0] = horiz;
corners[pt_cnt][1] = vert;
pt_cnt++;
}
}
} // end of search for corners
if (pt_cnt == 0) return -9999; // no corners found
// Find two points that form a rectangle with a corner
int pt;
double height = 0;
for ( pt = 0; pt < pt_cnt; pt++ ) {
if ( !is_bad_point( xgrid+corners[pt][0],ygrid,bad_points ) &&
!is_bad_point( xgrid, ygrid+corners[pt][1],bad_points ) ) {
double test_horiz = centre_long + corners[pt][0]*col_step/3600;
double test_vert = centre_lat + corners[pt][1]*row_step/3600;
if ( check_points( test_horiz,centre_lat,centre_long,centre_lat ) &&
check_points( centre_long,test_vert,centre_long,centre_lat ) ) break;
}
}
if (pt == pt_cnt) {
// perhaps we have a concave cliff, just take the
// average of the known points
double totht = 0;
for( int pti = 0; pti <pt_cnt; pti++ ) {
totht = totht + get_array_elev( xgrid+corners[pti][0],ygrid+corners[pti][1] );
}
height = totht/pt_cnt;
} else {
// We have three points, calculate the height
// Set anything very negative to zero
double corner = get_array_elev( xgrid+corners[pt][0],ygrid+corners[pt][1] );
double horiz = get_array_elev( xgrid,ygrid+corners[pt][1] );
double vert = get_array_elev( xgrid+corners[pt][0],ygrid );
if ( corner < -9000 ) corner = 0;
if ( horiz < -9000 ) horiz = 0;
if ( vert < -9000 ) vert = 0;
height = horiz + ( vert - corner );
}
SG_LOG(SG_GENERAL, SG_DEBUG, xgrid << "," << ygrid << ": was " << original_height << " , now " << height);
return height;
}
// return the current altitude based on grid data.
// TODO: We should rewrite this to interpolate exact values, but for now this is good enough
double TGArray::altitude_from_grid( double lon, double lat ) const {
@ -360,16 +600,17 @@ double TGArray::altitude_from_grid( double lon, double lat ) const {
------
then calculate our end points
*/
*/
// Store in degrees for later
double londeg = lon/3600;
double latdeg = lat/3600;
xlocal = (lon - originx) / col_step;
ylocal = (lat - originy) / row_step;
xindex = (int)(xlocal);
yindex = (int)(ylocal);
// printf("xindex = %d yindex = %d\n", xindex, yindex);
if ( xindex + 1 == cols ) {
xindex--;
}
@ -380,75 +621,274 @@ double TGArray::altitude_from_grid( double lon, double lat ) const {
if ( (xindex < 0) || (xindex + 1 >= cols) ||
(yindex < 0) || (yindex + 1 >= rows) ) {
SG_LOG(SG_GENERAL, SG_DEBUG, "WARNING: Attempt to interpolate value outside of array!!!" );
return -9999;
}
dx = xlocal - xindex;
dy = ylocal - yindex;
// Now check if we are on the same side of any cliffs
if ( dx > dy ) {
// lower triangle
// Collect lat,long at corners of area
// remember the missing corner if three found
// Go around the rectangle clockwise from SW corner
int corners[4][2];
int ccnt = 0;
int missing = -1; //the missing point when 3 available
double lon1 = (originx+(xindex*col_step))/3600;
double lat1 = (originy+(yindex*row_step))/3600;
double lon2 = lon1 + col_step/3600;
double lat2 = lat1 + row_step/3600;
if ( check_points(lon1,lat1,londeg,latdeg) ) {
corners[ccnt][0] = xindex;
corners[ccnt][1] = yindex;
ccnt++;
} else missing = 0;
if ( check_points(lon1,lat2,londeg,latdeg) ) {
corners[ccnt][0] = xindex;
corners[ccnt][1] = yindex+1;
ccnt++;
} else missing = 1;
if ( check_points(lon2,lat2,londeg,latdeg) ) {
corners[ccnt][0] = xindex+1;
corners[ccnt][1] = yindex+1;
ccnt++;
} else missing = 2;
if ( check_points(lon2,lat1,londeg,latdeg) ) {
corners[ccnt][0] = xindex+1;
corners[ccnt][1] = yindex;
ccnt++;
} else missing = 3;
switch (ccnt) {
case 3: //3 points are corners of a rectangle
// choose the points so that x2 is the right angle
// and x1-x2 is the x arm of the triangle
// dx,dy are the (positive) distances from the x1 corner
SG_LOG(SG_GENERAL, SG_DEBUG, "3 points, missing #" << missing);
dx = xlocal -xindex;
dy = ylocal -yindex;
switch ( missing ) {
case 0: //SW corner missing
x1 = corners[0][0];
y1 = corners[0][1];
x1 = xindex;
y1 = yindex;
z1 = get_array_elev(x1, y1);
x2 = corners[1][0];
y2 = corners[1][1];
x2 = xindex + 1;
y2 = yindex;
z2 = get_array_elev(x2, y2);
x3 = corners[2][0];
y3 = corners[2][1];
x3 = xindex + 1;
y3 = yindex + 1;
z3 = get_array_elev(x3, y3);
dy = 1 - dy;
break;
case 1: //NW corner missing
x1 = corners[0][0];
y1 = corners[0][1];
x2 = corners[2][0];
y2 = corners[2][1];
x3 = corners[1][0];
y3 = corners[1][1];
break;
case 2: //NE corner missing
x1 = corners[2][0];
y1 = corners[2][1];
x2 = corners[0][0];
y2 = corners[0][1];
x3 = corners[1][0];
y3 = corners[1][1];
dx = 1 - dx; //x1 is SE corner
break;
case 3: //SE corner missing
x1 = corners[2][0];
y1 = corners[2][1];
x2 = corners[1][0];
y2 = corners[1][1];
x3 = corners[0][0];
y3 = corners[0][1];
dx = 1 - dx; //x1 is NE corner
dy = 1 - dy;
break;
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;
if ( dx > SG_EPSILON ) {
elev = dy * (zB - zA) / dx + zA;
} else {
elev = zA;
}
} else {
// upper triangle
x1 = xindex;
y1 = yindex;
z1 = get_array_elev(x1, y1);
x2 = xindex;
y2 = yindex + 1;
z2 = get_array_elev(x2, y2);
x3 = xindex + 1;
y3 = yindex + 1;
z3 = get_array_elev(x3, y3);
if ( z1 < -9000 || z2 < -9000 || z3 < -9000 ) {
// don't interpolate off a void
return closest_nonvoid_elev( lon, lat );
// Now do the calcs on the triangle
// We interpolate on height along x1-x2 and
// x1 - x3. Then interpolate between these
// two points along y.
z1 = get_array_elev(x1,y1);
z2 = get_array_elev(x2,y2);
z3 = get_array_elev(x3,y3);
zA = dx * (z2 - z1) + z1;
zB = dx * (z3 - z1) + z1;
if ( dx > SG_EPSILON ) {
elev = dy * (zB - zA) / dx + zA;
} else {
elev = zA;
}
break;
case 2: //project onto line connecting two points
x1 = corners[0][0];
y1 = corners[0][1];
z1 = get_array_elev(x1,y1);
zA = dy * (z2 - z1) + z1;
zB = dy * (z3 - z1) + z1;
x2 = corners[1][0];
y2 = corners[1][1];
z2 = get_array_elev(x2,y2);
if ( dy > SG_EPSILON ) {
elev = dx * (zB - zA) / dy + zA;
} else {
elev = zA;
}
//two points are either a side of the rectangle, or
//else the diagonal
dx = xlocal - x1;
dy = ylocal - y1;
if (x1==x2) {
elev = z1+dy*(z2-z1)/(y2-y1);
}
else if (y1==y2) {
elev = z1+dx*(z2-z1)/(x2-x1);
}
else { //diagonal: project onto 45 degree line
int comp1 = x2-x1;
int comp2 = y2-y1;
double dotprod = (dx*comp1 + dy*comp2)/sqrt(2);
double projlen = sqrt(dx*dx+dy*dy)*dotprod;
elev = (z2-z1)*projlen/sqrt(2);
}
break;
case 1: //only one point found
elev = get_array_elev( corners[0][0],corners[0][1] );
break;
case 0: // all points on wrong side, fall through to normal calc
SG_LOG(SG_GENERAL, SG_WARN, "All elevation grid points on wrong side of cliff for " << std::setprecision(10) << londeg << "," << latdeg );
SG_LOG(SG_GENERAL, SG_WARN, "Grid points ("<< std::setprecision(9) << lon1 << "," << lat1 << "),("<<lon2<<","<<lat2<<")");
default: // all corners
dx = xlocal - xindex;
dy = ylocal - yindex;
if ( dx > dy ) {
// lower triangle
x1 = xindex;
y1 = yindex;
z1 = get_array_elev(x1, y1);
x2 = xindex + 1;
y2 = yindex;
z2 = get_array_elev(x2, y2);
x3 = xindex + 1;
y3 = yindex + 1;
z3 = get_array_elev(x3, y3);
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;
if ( dx > SG_EPSILON ) {
elev = dy * (zB - zA) / dx + zA;
} else {
elev = zA;
}
} else {
// upper triangle
x1 = xindex;
y1 = yindex;
z1 = get_array_elev(x1, y1);
x2 = xindex;
y2 = yindex + 1;
z2 = get_array_elev(x2, y2);
x3 = xindex + 1;
y3 = yindex + 1;
z3 = get_array_elev(x3, y3);
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;
if ( dy > SG_EPSILON ) {
elev = dx * (zB - zA) / dy + zA;
} else {
elev = zA;
}
}
}
return elev;
}
// Check that two points are on the same side of all cliff contours
// Could speed up by checking bounding box first
bool TGArray::check_points( const double lon1, const double lat1, const double lon2, const double lat2 ) const {
if ( cliffs_list.size()==0 ) return true;
if ( fabs(lon1-lon2)<SG_EPSILON && fabs(lat1-lat2)<SG_EPSILON ) return true;
SGGeod pt1 = SGGeod::fromDeg( lon1,lat1 );
SGGeod pt2 = SGGeod::fromDeg( lon2,lat2 );
bool same_side = true;
for ( int i=0;i<cliffs_list.size();i++ ) {
bool check_result = cliffs_list[i].AreSameSide( pt1,pt2 );
if(!check_result) {
SG_LOG(SG_GENERAL, SG_DEBUG, "Cliff " << i <<":" <<pt1 << " and " << pt2 << " on opposite sides");
same_side = false;
break;
}
}
return same_side;
}
//Check that a point is more than given distance from any cliff
//Could speed up by checking bounding box
bool TGArray::is_near_cliff( const double lon1, const double lat1, const double bad_zone ) const {
if (cliffs_list.size()==0) return false;
SGGeod pt1 = SGGeod::fromDeg(lon1,lat1);
for ( int i=0;i<cliffs_list.size();i++ ) {
double dist = cliffs_list[i].MinDist(pt1);
if (dist < bad_zone) return true;
}
return false;
}
TGArray::~TGArray( void )
{
if (in_data) {
@ -480,9 +920,10 @@ void TGArray::set_array_elev( int col, int row, int val )
bool TGArray::is_open() const
{
if ( array_in != NULL ) {
return true;
} else {
return false;
}
if ( array_in != NULL ) {
return true;
} else {
return false;
}
}

26
src/Lib/Array/array.hxx
View file

@ -27,6 +27,11 @@
#include <simgear/bucket/newbucket.hxx>
#include <simgear/math/sg_types.hxx>
#include <simgear/io/iostreams/sgstream.hxx>
#include <simgear/misc/sg_dir.hxx>
#include <boost/foreach.hpp>
#include "tg_contour.hxx"
#include "tg_polygon.hxx"
class TGArray {
@ -42,6 +47,8 @@ private:
// number of columns and rows
int cols, rows;
// Whether or not the input data have been rectified
bool rectified;
// Distance between column and row data points (in arc seconds)
double col_step, row_step;
@ -52,7 +59,15 @@ private:
std::vector<SGGeod> corner_list;
std::vector<SGGeod> fitted_list;
// list of cliff contours
tgcontour_list cliffs_list;
void parse_bin();
// Routines for height rectification
std::vector<int> collect_bad_points(const double bad_zone);
bool is_bad_point(const int xgrid, const int ygrid, const std::vector<int> bad_points) const;
double rectify_point(const int xgrid, const int ygrid, const std::vector<int> bad_points) const;
bool is_near_cliff(const double lon1,const double lon2, const double bad_zone) const;
public:
// Constructor
@ -65,6 +80,9 @@ public:
// open an Array file (use "-" if input is coming from stdin)
bool open ( const std::string& file_base );
// Load contours from polygon files delineating height discontinuities
void load_cliffs(const std::string & height_base);
// return if array was successfully opened or not
bool is_open() const;
@ -77,10 +95,16 @@ public:
// write an Array file
bool write( const std::string root_dir, SGBucket& b );
// write an Array file in binary format. If ht_rect is true,
// the file will have extension 'arr.rectified.gz'
void write_bin(const std::string root_dir, bool ht_rect, SGBucket& b);
// do our best to remove voids by picking data from the nearest
// neighbor.
void remove_voids();
void rectify_heights(const double bad_zone);
// Return the elevation of the closest non-void grid point to lon, lat
double closest_nonvoid_elev( double lon, double lat ) const;
@ -103,6 +127,8 @@ public:
int get_array_elev( int col, int row ) const;
void set_array_elev( int col, int row, int val );
// Check whether or not two points are on the same side of contour
bool check_points (const double a,const double b, const double c, const double d) const;
// reset Array to initial state - ready to load another elevation file
void unload( void );
};

153
src/Lib/Array/rectify_height.cxx Normal file
View file

@ -0,0 +1,153 @@
// rectify_height.cxx -- Correct height grid for incorrect heights near
// discontinuities included in cliff files
//
// Written by James Hester 2018
//
// Copyright (C) 2018 James Hester
//
// 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., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
//
#ifdef HAVE_CONFIG_H
# include <config.h>
#endif
#include <simgear/compiler.h>
#include <simgear/bucket/newbucket.hxx>
#include <simgear/misc/sg_path.hxx>
#include <simgear/debug/logstream.hxx>
#include <sys/types.h>
#include <iostream>
#include <string.h>
#include <stdlib.h>
#include "array.hxx"
/* This program will find all height files in the directory provided,
and if they have an associated cliff file, will adjust and then output
the heights. */
// display usage and exit
static void usage( const std::string name ) {
SG_LOG(SG_GENERAL, SG_ALERT, "Usage: " << name);
SG_LOG(SG_GENERAL, SG_ALERT, " --work-dir=<directory>");
SG_LOG(SG_GENERAL, SG_ALERT, " --height-dir=<directory>");
SG_LOG(SG_GENERAL, SG_ALERT, "[ --tile-id=<id>]");
SG_LOG(SG_GENERAL, SG_ALERT, " --min-lon=<degrees>");
SG_LOG(SG_GENERAL, SG_ALERT, " --max-lon=<degrees>");
SG_LOG(SG_GENERAL, SG_ALERT, " --min-lat=<degrees>");
SG_LOG(SG_GENERAL, SG_ALERT, " --max-lat=<degrees>");
SG_LOG(SG_GENERAL, SG_ALERT, "[ --min-dist=<float>]");
exit(-1);
}
int main( int argc, char **argv) {
std::string output_dir = ".";
std::string work_dir = ".";
std::string height_dir = "";
SGGeod min,max;
long tile_id = -1;
double bad_zone = 100; //distance in m from cliff to rectify
sglog().setLogLevels(SG_ALL,SG_INFO);
sglog().set_log_priority(SG_DEBUG);
//
// Parse the command-line arguments.
//
int arg_pos;
for (arg_pos = 1; arg_pos < argc; arg_pos++) {
std::string arg = argv[arg_pos];
if (arg.find("--work-dir=") == 0) {
work_dir = arg.substr(11);
} else if (arg.find("--height-dir=") == 0) {
height_dir = arg.substr(13);
} else if (arg.find("--tile-id=") == 0) {
tile_id = atol(arg.substr(10).c_str());
} else if ( arg.find("--min-lon=") == 0 ) {
min.setLongitudeDeg(atof( arg.substr(10).c_str() ));
} else if ( arg.find("--max-lon=") == 0 ) {
max.setLongitudeDeg(atof( arg.substr(10).c_str() ));
} else if ( arg.find("--min-lat=") == 0 ) {
min.setLatitudeDeg(atof( arg.substr(10).c_str() ));
} else if ( arg.find("--max-lat=") == 0 ) {
max.setLatitudeDeg(atof( arg.substr(10).c_str() ));
} else if ( arg.find("--min-dist=") == 0) {
bad_zone = atof(arg.substr(11).c_str());
} else if (arg.find("--") == 0) {
usage(argv[0]);
} else {
break;
}
}
SG_LOG(SG_GENERAL, SG_ALERT, "Working directory is " << work_dir);
SG_LOG(SG_GENERAL, SG_ALERT, "Heights are in " << height_dir);
SG_LOG(SG_GENERAL, SG_ALERT, "Rectification zone within " << bad_zone << "m of cliffs");
if ( tile_id > 0 ) {
SG_LOG(SG_GENERAL, SG_ALERT, "Tile id is " << tile_id);
} else {
if (min.isValid() && max.isValid() && (min != max))
{
SG_LOG(SG_GENERAL, SG_ALERT, "Longitude = " << min.getLongitudeDeg() << ':' << max.getLongitudeDeg());
SG_LOG(SG_GENERAL, SG_ALERT, "Latitude = " << min.getLatitudeDeg() << ':' << max.getLatitudeDeg());
} else
{
SG_LOG(SG_GENERAL, SG_ALERT, "Lon/Lat unset or wrong");
exit(1);
}
}
// Now generate the buckets
std::vector<SGBucket> bucketList;
if (tile_id == -1) {
// build all the tiles in an area
SG_LOG(SG_GENERAL, SG_ALERT, "Fixing heights within given bounding box");
SGBucket b_min( min );
SGBucket b_max( max );
if ( b_min == b_max ) {
bucketList.push_back( b_min );
} else {
SG_LOG(SG_GENERAL, SG_ALERT, " adjustment area spans tile boundaries");
sgGetBuckets( min, max, bucketList );
}
} else {
// adjust the specified tile
SG_LOG(SG_GENERAL, SG_ALERT, "Adjusting tile " << tile_id);
bucketList.push_back( SGBucket( tile_id ) );
}
// Finally, loop over the buckets adjusting heights
for (unsigned int i=0; i < bucketList.size(); i++) {
SGBucket bucket = bucketList[i];
TGArray array;
std::string base_path = bucket.gen_base_path();
std::string array_path = work_dir + "/" + height_dir + "/" + base_path + "/" + bucket.gen_index_str();
if (!array.open(array_path)) {
SG_LOG(SG_GENERAL,SG_DEBUG, "Failed to open array file " << array_path);
continue;
}
SG_LOG(SG_GENERAL,SG_DEBUG, "Successfully opened array file");
array.parse(bucket);
array.close();
array.remove_voids();
array.rectify_heights(bad_zone);
array.write_bin(work_dir + "/" + height_dir,true,bucket);
}
}

4
src/Lib/Array/testarray.cxx
View file

@ -7,6 +7,7 @@
#include <simgear/compiler.h>
#include <simgear/bucket/newbucket.hxx>
#include <simgear/misc/sg_path.hxx>
#include <simgear/debug/logstream.hxx>
#include <sys/types.h>
#include <iostream>
@ -42,6 +43,9 @@ int main( int argc, char **argv ) {
double lon, lat;
check_for_help(argc, argv);
sglog().setLogLevels(SG_ALL,SG_INFO);
sglog().set_log_priority(SG_DEBUG);
if ( argc != 4 ) {
cout << "ERROR: Needs 3 arguments!" << endl;

1
src/Lib/HGT/CMakeLists.txt
View file

@ -1,4 +1,5 @@
add_library(HGT STATIC
dted.cxx dted.hxx
hgt.cxx hgt.hxx
srtmbase.cxx srtmbase.hxx
)

232
src/Lib/HGT/dted.cxx Normal file
View file

@ -0,0 +1,232 @@
// dted.cxx -- SRTM "dted" data management class
//
// Written by James Hester based on hgt code of
// Curtis Olson, started February 2003.
//
// Copyright (C) 2003 Curtis L. Olson - http://www.flightgear.org/~curt
// Copyright (C) 2018 James Hester
//
// 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., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
//
// $Id: hgt.cxx,v 1.7 2005-12-19 16:06:45 curt Exp $
#ifdef HAVE_CONFIG_H
# include <config.h>
#endif
#include <simgear/compiler.h>
#include <stdlib.h> // atof()
#include <iostream>
#ifdef SG_HAVE_STD_INCLUDES
# include <cerrno>
#else
# include <errno.h>
#endif
#ifdef _MSC_VER
# include <direct.h>
#endif
#include <simgear/constants.h>
#include <simgear/io/lowlevel.hxx>
#include <simgear/misc/sg_dir.hxx>
#include <simgear/debug/logstream.hxx>
#include "dted.hxx"
using std::cout;
using std::endl;
using std::string;
TGDted::TGDted( int _res )
{
dted_resolution = _res;
data = new short int[MAX_DTED_SIZE][MAX_DTED_SIZE];
output_data = new short int[MAX_DTED_SIZE][MAX_DTED_SIZE];
}
TGDted::TGDted( int _res, const SGPath &file ):TGDted(_res)
{
TGDted::open( file );
}
// open an DTED file
bool
TGDted::open ( const SGPath &f ) {
SGPath file_name = f;
// open input file (or read from stdin)
if ( file_name.str() == "-" ) {
cout << "Loading DTED data file: stdin" << endl;
if ( (fd = gzdopen(0, "r")) == NULL ) { // 0 == STDIN_FILENO
cout << "ERROR: opening stdin" << endl;
return false;
}
} else {
if ( file_name.extension() == "zip" ) {
// extract the .zip file to /tmp and point the file name
// to the extracted file
tmp_dir = simgear::Dir::tempDir("dted");
cout << "Extracting " << file_name.str() << " to " << tmp_dir.path().str() << endl;
string command = "unzip -d \"" + tmp_dir.path().str() + "\" " + file_name.base();
if ( system( command.c_str() ) != -1 )
{
simgear::PathList files = tmp_dir.children(simgear::Dir::TYPE_FILE | simgear::Dir::NO_DOT_OR_DOTDOT);
for (const SGPath& file : files) {
string ext = file.lower_extension();
if ( ext == "dted" ) {
file_name = file;
break;
}
}
remove_tmp_file = true;
cout << "Proceeding with " << file_name.str() << endl;
} else {
SG_LOG(SG_GENERAL, SG_ALERT, "Failed to issue system call " << command );
return EXIT_FAILURE;
}
}
cout << "Loading DTED 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 contents
// User Header Label
char header[3];
char location[7];
int degrees,minutes,seconds;
char hemisphere;
// Check header
gzread(fd,header,3);
if (strncmp(header,"UHL",3) != 0) {
cout << "UHL User Header Label not found" << endl;
return false;
}
gzread(fd,header,1); //dummy
gzread(fd,location,8);//longitude
sscanf(location,"%3d%2d%2d%c",&degrees,&minutes,&seconds,&hemisphere);
originx = degrees *3600 + minutes*60 + seconds;
if(hemisphere == 'W') {
originx = - originx;
}
gzread(fd,location,8);//latitude
sscanf(location,"%3d%2d%2d%c",&degrees,&minutes,&seconds,&hemisphere);
originy = degrees *3600 + minutes*60 + seconds;
if(hemisphere == 'S') {
originy = - originy;
}
cout << " Origin = " << originx << ", " << originy << endl;
return true;
}
// close an DTED file
bool
TGDted::close () {
gzclose(fd);
return true;
}
// load an DTED file
bool
TGDted::load( ) {
int size;
bool little_endian = sgIsLittleEndian();
if ( dted_resolution == 1 ) {
cols = rows = size = 3601;
col_step = row_step = 1;
} else if ( dted_resolution == 3 ) {
cols = rows = size = 1201;
col_step = row_step = 3;
} else {
cout << "Unknown DTED resolution, only 1 and 3 arcsec formats" << endl;
cout << " are supported!" << endl;
return false;
}
if (little_endian) {
cout << "Little Endian: swapping input values" << endl;
}
//Skip headers
gzseek(fd,3428,SEEK_SET);
unsigned short int latct, longct;
short int *var;
int dummy; //to read column header
for ( int col = 0; col < size; ++col ) {
dummy = 0; // zero out all bytes
longct = 0;
latct = 0;
gzread(fd,&dummy,1); //sentinel
if(dummy != 170) {
cout << "Failed to find sentinel at col " << col << endl;
return false;
}
gzread(fd,&dummy,3); //block count
gzread(fd,&longct,2); //Longitude count
gzread(fd,&latct,2); //Latitude count
if ( little_endian ) {
sgEndianSwap(&longct);
}
// cout << "Longitude count " << longct << endl;
for ( int row = 0; row < size; ++row ) {
var = &data[col][row];
if ( gzread ( fd, var, 2 ) != sizeof(short) ) {
return false;
}
if ( little_endian ) {
sgEndianSwap( (unsigned short int*)var);
}
}
gzread(fd,&dummy,4); //Checksum
// Check values are right
if (col == 0) {
cout << data[col][0] << endl;
cout << data[col][1] << endl;
cout << data[col][2] << endl;
}
}
return true;
}
TGDted::~TGDted() {
// printf("class TGSrtmBase DEstructor called.\n");
delete [] data;
delete [] output_data;
}

84
src/Lib/HGT/dted.hxx Normal file
View file

@ -0,0 +1,84 @@
// dted.hxx -- SRTM "dted" data management class
//
// Written by James Hester based on hgt.hxx, which was
// written by Curtis Olson who started February 2003.
//
// Copyright (C) 2003 Curtis L. Olson - http://www.flightgear.org/~curt
// Copyright (c) 2018 James Hester
//
// 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., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
//
#ifndef _DTED_HXX
#define _DTED_HXX
#ifdef HAVE_CONFIG_H
# include <config.h>
#endif
#include <simgear/compiler.h>
#include "srtmbase.hxx"
#include <zlib.h>
#include <string>
#include <simgear/bucket/newbucket.hxx>
#include <simgear/misc/sg_path.hxx>
#define MAX_DTED_SIZE 3601
class TGDted : public TGSrtmBase {
private:
// file pointer for input
gzFile fd;
int dted_resolution;
// pointers to the actual grid data allocated here
short int (*data)[MAX_DTED_SIZE];
short int (*output_data)[MAX_DTED_SIZE];
public:
// Constructor, _res must be either "1" for the 1arcsec data or
// "3" for the 3arcsec data.
TGDted( int _res );
TGDted( int _res, const SGPath &file );
// Destructor
~TGDted();
// open an DTED file (use "-" if input is coming from stdin)
bool open ( const SGPath &file );
// close an DTED file
bool close();
// load an dted file
bool load();
virtual short height( int x, int y ) const { return data[x][y]; }
};
#endif // _DTED_HXX

3
src/Lib/terragear/CMakeLists.txt
View file

@ -1,4 +1,5 @@
include_directories(${GDAL_INCLUDE_DIR})
include_directories(${PROJECT_SOURCE_DIR}/src/Lib/terragear)
include( ${CGAL_USE_FILE} )
@ -32,4 +33,4 @@ add_library(terragear STATIC
tg_unique_vec2f.hxx
tg_unique_vec3d.hxx
tg_unique_vec3f.hxx
)
)

12
src/Lib/terragear/tg_chopper.cxx
View file

@ -11,8 +11,7 @@
#include "tg_misc.hxx"
tgPolygon tgChopper::Clip( const tgPolygon& subject,
const std::string& type,
SGBucket& b )
const std::string& type, SGBucket& b )
{
tgPolygon base, result;
@ -22,8 +21,9 @@ tgPolygon tgChopper::Clip( const tgPolygon& subject,
base.AddNode( 0, b.get_corner( SG_BUCKET_NE ) );
base.AddNode( 0, b.get_corner( SG_BUCKET_NW ) );
result = tgPolygon::Intersect( subject, base );
result = tgPolygon::Intersect( subject, base );
if ( result.Contours() > 0 ) {
if ( subject.GetPreserve3D() ) {
result.InheritElevations( subject );
result.SetPreserve3D( true );
@ -34,6 +34,9 @@ tgPolygon tgChopper::Clip( const tgPolygon& subject,
result.SetTexMethod( TG_TEX_BY_GEODE, b.get_center_lat() );
}
result.SetFlag(type);
if (!subject.IsClosed()) {
result.SetOpen();
}
lock.lock();
bp_map[b.gen_index()].push_back( result );
@ -63,6 +66,7 @@ void tgChopper::ClipRow( const tgPolygon& subject, const double& center_lat, con
}
}
void tgChopper::Add( const tgPolygon& subject, const std::string& type )
{
// bail out immediately if polygon is empty
@ -205,7 +209,7 @@ void tgChopper::Save( bool DebugShapefiles )
long int poly_index = GenerateIndex( path );
sprintf( poly_ext, "%ld", poly_index );
polyfile = polyfile + "." + poly_ext;
polyfile = polyfile + "." + poly_ext + "." + extra_extension;
gzFile fp;
if ( (fp = gzopen( polyfile.c_str(), "wb9" )) == NULL ) {

5
src/Lib/terragear/tg_chopper.hxx
View file

@ -11,10 +11,14 @@ class tgChopper
public:
tgChopper( const std::string& path ) {
root_path = path;
extra_extension = "";
}
void Add( const tgPolygon& poly, const std::string& type );
void Save( bool DebugShapes );
void Add_Extension( const std::string& extension) {
extra_extension = extension;
}
private:
long int GenerateIndex( std::string path );
@ -25,4 +29,5 @@ private:
std::string root_path;
bucket_polys_map bp_map;
SGMutex lock;
std::string extra_extension; //add at end of file name
};

124
src/Lib/terragear/tg_contour.cxx
View file

@ -1,4 +1,5 @@
#include <simgear/math/sg_geodesy.hxx>
#include <simgear/math/SGGeometry.hxx>
#include <simgear/io/lowlevel.hxx>
#include <simgear/debug/logstream.hxx>
@ -65,7 +66,6 @@ double tgContour::GetArea( void ) const
for (i=0; i<node_list.size(); i++) {
a = SGGeod_ToSGVec2d( node_list[i] );
b = SGGeod_ToSGVec2d( node_list[j] );
area += (b.x() + a.x()) * (b.y() - a.y());
j=i;
}
@ -74,6 +74,93 @@ double tgContour::GetArea( void ) const
return fabs(area * 0.5);
}
// Check that the two supplied points are on the same side of the contour
bool tgContour::AreSameSide( const SGGeod& firstpt, const SGGeod& secondpt) const
{
//Find equation of line segment joining the points
double x1 = firstpt.getLatitudeDeg();
double x2 = secondpt.getLatitudeDeg();
double y1 = firstpt.getLongitudeDeg();
double y2 = secondpt.getLongitudeDeg();
//Store differences for later
double xdif = x2-x1;
double ydif = y2-y1;
/*We describe a line parametrically:
x1 (x2-x1)
L = + t
y1 (y2-y1)
with u the parametric coefficient for the second line.
Then the line segments intersect if 0 <= t,u <= 1.
To determine t and u we use the approach of Goldman ("Graphics
Gems" as described in Stack Overflow question 563198).
if r x s = r_x * s_y - r_y * s_x, then
t = (q - p) x s / (r x s)
and
u = (q - p) x r / (r x s)
for line 1 = p + t r, line 2 = q + u s
*/
//Now cycle over all nodes and count how many times we intersect
int intersect_ct = 0;
if (node_list.size()) {
int j = node_list.size() - 1;
for (int i=0;i<node_list.size()-1;i++) {
double nx1 = node_list[i].getLatitudeDeg();
double ny1 = node_list[i].getLongitudeDeg();
double nx2 = node_list[i+1].getLatitudeDeg();
double ny2 = node_list[i+1].getLongitudeDeg();
double nydif = ny2-ny1;
double nxdif = nx2-nx1;
double denom = xdif*nydif - ydif*nxdif;
if (denom != 0) { //Not parallel
double crossx = nx1-x1; double crossy = ny1-y1;
double t = (crossx*nydif - crossy*nxdif)/denom;
double u = -1*(xdif*crossy - ydif*crossx)/denom;
// We consider that an intersection at the edge of the line has
// crossed
// over, that is, they lie on opposite sides. This way we capture
// places where the chopper has clipped a cliff on the tile edge
if (t > -0.0001 && t < 1.0001 && u > -0.0001 && u < 1.0001) intersect_ct++;
}
}
}
bool isinter = (intersect_ct%2 == 0);
return isinter;
}
double tgContour::MinDist(const SGGeod& probe) const {
SGVec3d probexyz;
SGGeodesy::SGGeodToCart( probe,probexyz );
double mindist = 100000.0;
double dist;
if ( node_list.size() ) {
int j = node_list.size() - 1;
for (int i=0;i<j;i++) {
SGVec3d start,end;
SGGeodesy::SGGeodToCart( node_list[i],start );
SGGeodesy::SGGeodToCart( node_list[i+1],end );
SGLineSegment<double> piece = SGLineSegment<double>(start,end);
dist = distSqr( piece,probexyz );
if (dist < mindist) mindist = dist;
}
}
return sqrt(mindist);
}
bool tgContour::IsInside( const tgContour& inside, const tgContour& outside )
{
// first contour is inside second if the intersection of first with second is == first
@ -424,6 +511,7 @@ tgContour tgContour::FromClipper( const ClipperLib::Path& subject )
return result;
}
tgRectangle tgContour::GetBoundingBox( void ) const
{
SGGeod min, max;
@ -803,14 +891,16 @@ tgContour tgContour::AddColinearNodes( const tgContour& subject, std::vector<SGG
p0 = subject.GetNode( subject.GetSize() - 1 );
p1 = subject.GetNode( 0 );
// add start of segment
result.AddNode( p0 );
// add intermediate points
AddIntermediateNodes( p0, p1, nodes, result, SG_EPSILON*10, SG_EPSILON*4 );
// maintain original hole flag setting
if (!subject.GetOpen()) {
// add start of segment
result.AddNode( p0 );
// add intermediate points
AddIntermediateNodes( p0, p1, nodes, result, SG_EPSILON*10, SG_EPSILON*4 );
}
// maintain original hole and openness flag setting
result.SetHole( subject.GetHole() );
result.SetOpen( subject.GetOpen() );
return result;
}
@ -833,15 +923,17 @@ tgContour tgContour::AddColinearNodes( const tgContour& subject, bool preserve3d
p0 = subject.GetNode( subject.GetSize() - 1 );
p1 = subject.GetNode( 0 );
// add start of segment
result.AddNode( p0 );
// add intermediate points
AddIntermediateNodes( p0, p1, preserve3d, nodes, result, SG_EPSILON*10, SG_EPSILON*4 );
// maintain original hole flag setting
if(!subject.GetOpen()) {
// add start of segment
result.AddNode( p0 );
// add intermediate points
AddIntermediateNodes( p0, p1, preserve3d, nodes, result, SG_EPSILON*10, SG_EPSILON*4 );
}
// maintain original hole and open flag settings
result.SetHole( subject.GetHole() );
result.SetOpen( subject.GetOpen() );
return result;
}

19
src/Lib/terragear/tg_contour.hxx
View file

@ -7,6 +7,8 @@
#include <simgear/compiler.h>
#include <simgear/math/sg_types.hxx>
#include <simgear/math/SGLineSegment.hxx>
#include <simgear/math/SGGeodesy.hxx>
#include <boost/concept_check.hpp>
#include "tg_unique_geod.hxx"
@ -29,6 +31,7 @@ class tgContour
public:
tgContour() {
hole = false;
keep_open = false;
}
void Erase() {
@ -38,10 +41,19 @@ public:
void SetHole( bool h ) {
hole = h;
}
bool GetHole( void ) const {
return hole;
}
bool GetOpen(void) const {
return keep_open;
}
void SetOpen(bool o) {
keep_open = o;
}
unsigned int GetSize( void ) const {
return node_list.size();
}
@ -101,6 +113,10 @@ public:
}
// Return true if the two points are on the same side of the contour
bool AreSameSide(const SGGeod& firstpt, const SGGeod& secondpt) const;
// Return minimum distance of point from contour
double MinDist(const SGGeod& probe) const;
static tgContour Snap( const tgContour& subject, double snap );
static tgContour RemoveDups( const tgContour& subject );
static tgContour SplitLongEdges( const tgContour& subject, double dist );
@ -135,10 +151,11 @@ public:
private:
std::vector<SGGeod> node_list;
bool hole;
bool keep_open; //If non-closed contour, keep open
};
typedef std::vector <tgContour> tgcontour_list;
typedef tgcontour_list::iterator tgcontour_list_iterator;
typedef tgcontour_list::const_iterator const_tgcontour_list_iterator;
#endif // _TGCONTOUR_HXX
#endif // _TGCONTOUR_HXX

15
src/Lib/terragear/tg_polygon.hxx
View file

@ -243,6 +243,7 @@ class tgPolygon
public:
tgPolygon() {
preserve3d = false;
closed = true;
}
~tgPolygon() {
contours.clear();
@ -351,6 +352,18 @@ public:
flag = f;
}
void SetOpen( void ) { //Do not try to close the contours
closed = false;
}
void SetClosed (void) {
closed = true;
}
bool IsClosed( void ) const {
return closed;
}
bool GetPreserve3D( void ) const {
return preserve3d;
}
@ -436,6 +449,7 @@ public:
// Conversions
static ClipperLib::Paths ToClipper( const tgPolygon& subject );
static tgPolygon FromClipper( const ClipperLib::Paths& subject );
static tgPolygon FromClipper( const ClipperLib::PolyTree& subject );
static void ToClipperFile( const tgPolygon& subject, const std::string& path, const std::string& filename );
// T-Junctions and segment search
@ -460,6 +474,7 @@ private:
bool preserve3d;
unsigned int id; // unique polygon id for debug
tgTexParams tp;
bool closed; // if we treat it as a closed shape
};
#endif // _POLYGON_HXX

35
src/Lib/terragear/tg_polygon_clip.cxx
View file

@ -136,6 +136,7 @@ tgPolygon tgPolygon::Diff( const tgPolygon& subject, tgPolygon& clip )
return result;
}
//Intersect will keep open paths open
tgPolygon tgPolygon::Intersect( const tgPolygon& subject, const tgPolygon& clip )
{
tgPolygon result;
@ -156,22 +157,31 @@ tgPolygon tgPolygon::Intersect( const tgPolygon& subject, const tgPolygon& clip
ClipperLib::Paths clipper_subject = tgPolygon::ToClipper( subject );
ClipperLib::Paths clipper_clip = tgPolygon::ToClipper( clip );
ClipperLib::Paths clipper_result;
ClipperLib::Paths clipper_result; // for closed polygons
ClipperLib::PolyTree clipper_tree_result; //for open polygons
ClipperLib::Clipper c;
c.Clear();
c.AddPaths(clipper_subject, ClipperLib::PolyType::Subject, true);
c.AddPaths(clipper_subject, ClipperLib::PolyType::Subject, subject.IsClosed() );
c.AddPaths(clipper_clip, ClipperLib::PolyType::Clip, true);
c.Execute(ClipperLib::ClipType::Intersection, clipper_result, ClipperLib::PolyFillType::EvenOdd, ClipperLib::PolyFillType::EvenOdd);
result = tgPolygon::FromClipper( clipper_result );
if(subject.IsClosed()) {
c.Execute(ClipperLib::ClipType::Intersection, clipper_result, ClipperLib::PolyFillType::EvenOdd, ClipperLib::PolyFillType::EvenOdd );
result = tgPolygon::FromClipper( clipper_result );
}
else {
c.Execute(ClipperLib::ClipType::Intersection, clipper_tree_result, ClipperLib::PolyFillType::EvenOdd, ClipperLib::PolyFillType::EvenOdd );
result = tgPolygon::FromClipper( clipper_tree_result );
}
result = tgPolygon::AddColinearNodes( result, all_nodes );
result.SetMaterial( subject.GetMaterial() );
result.SetTexParams( subject.GetTexParams() );
result.SetId( subject.GetId() );
result.SetPreserve3D( subject.GetPreserve3D() );
if(subject.IsClosed()) {
result.SetClosed();
} else {
result.SetOpen();
}
return result;
}
@ -199,6 +209,17 @@ tgPolygon tgPolygon::FromClipper( const ClipperLib::Paths& subject )
return result;
}
// A polytree is returned when an open polygon has been clipped
tgPolygon tgPolygon::FromClipper( const ClipperLib::PolyTree& subject )
{
ClipperLib::Paths path_result;
tgPolygon result;
ClipperLib::PolyTreeToPaths(subject,path_result);
result = FromClipper(path_result);
return result;
}
ClipperLib::Path tgTriangle::ToClipper( const tgTriangle& subject )
{
ClipperLib::Path contour;

1
src/Prep/CMakeLists.txt
View file

@ -5,3 +5,4 @@ add_subdirectory(DemChop)
add_subdirectory(Terra)
add_subdirectory(TerraFit)
add_subdirectory(OGRDecode)
add_subdirectory(Cliff)

14
src/Prep/Cliff/CMakeLists.txt Normal file
View file

@ -0,0 +1,14 @@
include_directories(${GDAL_INCLUDE_DIR})
add_executable(cliff-decode cliff-decode.cxx)
target_link_libraries(cliff-decode
${GDAL_LIBRARY}
terragear
${ZLIB_LIBRARY}
${CMAKE_THREAD_LIBS_INIT}
${SIMGEAR_CORE_LIBRARIES}
${SIMGEAR_CORE_LIBRARY_DEPENDENCIES}
${Boost_LIBRARIES}
)
install(TARGETS cliff-decode RUNTIME DESTINATION bin)

435
src/Prep/Cliff/cliff-decode.cxx Normal file
View file

@ -0,0 +1,435 @@
// cliff-decode.cxx -- process OGR layers and extract lines
// representing discontinuities, clipping
// against and sorting
// them into the relevant tiles.
//
// Written by James Hester 2018
//
// Largely copied from ogr-decode by Ralf Gerlich.
//
// 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., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
//
#include <chrono>
#include <string>
#include <map>
#include <boost/thread.hpp>
#include <ogrsf_frmts.h>
#include <gdal_priv.h>
#include <simgear/compiler.h>
#include <simgear/threads/SGThread.hxx>
#include <simgear/threads/SGQueue.hxx>
#include <simgear/debug/logstream.hxx>
#include <simgear/math/sg_geodesy.hxx>
#include <simgear/misc/sg_path.hxx>
#include <Include/version.h>
#include <terragear/tg_chopper.hxx>
#include <terragear/tg_shapefile.hxx>
using std::string;
int continue_on_errors=0;
int num_threads=1;
bool use_attribute_query=false;
string attribute_query;
int start_record=0;
bool use_spatial_query=false;
int seperate_segments = 0;
double spat_min_x, spat_min_y, spat_max_x, spat_max_y;
bool save_shapefiles=false;
SGLockedQueue<OGRFeature *> global_workQueue;
class Decoder : public SGThread
{
public:
Decoder( OGRCoordinateTransformation *poct, tgChopper& c ) : chopper(c) {
poCT = poct;
}
private:
virtual void run();
void processLineString(OGRLineString* poGeometry);
private:
// The transformation for each geometry object
OGRCoordinateTransformation *poCT;
// Store the reults per tile
tgChopper& chopper;
};
void Decoder::processLineString(OGRLineString* poGeometry)
{
tgContour line;
SGGeod p0, p1;
double heading, dist, az2;
int i, j, numPoints, numSegs;
double max_dist;
numPoints = poGeometry->getNumPoints();
if (numPoints < 2) {
SG_LOG( SG_GENERAL, SG_WARN, "Skipping line with less than two points" );
return;
}
heading = SGGeodesy::courseDeg( p1, p0 );
// now add the middle points : if they are too far apart, add intermediate nodes
for ( i=0;i<numPoints;i++) {
p0 = SGGeod::fromDeg( poGeometry->getX(i), poGeometry->getY(i) );
line.AddNode( p0 );
}
// Do not close this contour
line.SetOpen(true);
// clip the contour.
// TODO: don't bother clipping, as the contour is informational
// only, and simply output to all relevant buckets instead.
tgPolygon open_poly;
open_poly.SetOpen(); //do not try to close this one up
open_poly.AddContour(line);
chopper.Add( open_poly, "Cliffs" );
}
void Decoder::run()
{
// as long as we have geometry to parse, do so
while (!global_workQueue.empty()) {
OGRFeature *poFeature = global_workQueue.pop();
SG_LOG( SG_GENERAL, SG_BULK, " remaining features is " << global_workQueue.size() );
if ( poFeature ) {
OGRGeometry *poGeometry = poFeature->GetGeometryRef();
if (poGeometry==NULL) {
SG_LOG( SG_GENERAL, SG_INFO, "Found feature without geometry!" );
if (!continue_on_errors) {
SG_LOG( SG_GENERAL, SG_ALERT, "Aborting!" );
exit( 1 );
} else {
continue;
}
}
OGRwkbGeometryType geoType=wkbFlatten(poGeometry->getGeometryType());
if (geoType!=wkbLineString && geoType!=wkbMultiLineString) {
SG_LOG( SG_GENERAL, SG_INFO, "Non-line feature" );
continue;
}
poGeometry->transform( poCT );
switch (geoType) {
case wkbLineString: {
SG_LOG( SG_GENERAL, SG_DEBUG, "LineString feature" );
processLineString((OGRLineString*)poGeometry);
break;
}
case wkbMultiLineString: {
SG_LOG( SG_GENERAL, SG_DEBUG, "MultiLineString feature" );
OGRMultiLineString* multilines=(OGRMultiLineString*)poGeometry;
for (int i=0;i<multilines->getNumGeometries();i++) {
processLineString((OGRLineString*)(multilines->getGeometryRef(i)));
}
break;
}
default:
/* Ignore unhandled objects */
break;
}
OGRFeature::DestroyFeature( poFeature );
}
}
}
// Main Thread
void processLayer(OGRLayer* poLayer, tgChopper& results )
{
int feature_count=poLayer->GetFeatureCount();
if (feature_count!=-1 && start_record>0 && start_record>=feature_count) {
SG_LOG( SG_GENERAL, SG_ALERT, "Layer has only " << feature_count << " records, but start record is set to " << start_record );
exit( 1 );
}
/* determine the indices of the required columns */
OGRFeatureDefn *poFDefn = poLayer->GetLayerDefn();
string layername=poFDefn->GetName();
/* setup a transformation to WGS84 */
OGRSpatialReference *oSourceSRS, oTargetSRS;
oSourceSRS=poLayer->GetSpatialRef();
if (oSourceSRS == NULL) {
SG_LOG( SG_GENERAL, SG_ALERT, "Layer " << layername << " has no defined spatial reference system" );
exit( 1 );
}
char* srsWkt;
oSourceSRS->exportToWkt(&srsWkt);
SG_LOG( SG_GENERAL, SG_DEBUG, "Source spatial reference system: " << srsWkt );
OGRFree(srsWkt);
oTargetSRS.SetWellKnownGeogCS( "WGS84" );
OGRCoordinateTransformation *poCT = OGRCreateCoordinateTransformation(oSourceSRS, &oTargetSRS);
/* setup attribute and spatial queries */
if (use_spatial_query) {
double trans_min_x,trans_min_y,trans_max_x,trans_max_y;
/* do a simple reprojection of the source SRS */
OGRCoordinateTransformation *poCTinverse;
poCTinverse = OGRCreateCoordinateTransformation(&oTargetSRS, oSourceSRS);
trans_min_x=spat_min_x;
trans_min_y=spat_min_y;
trans_max_x=spat_max_x;
trans_max_y=spat_max_y;
poCTinverse->Transform(1,&trans_min_x,&trans_min_y);
poCTinverse->Transform(1,&trans_max_x,&trans_max_y);
poLayer->SetSpatialFilterRect(trans_min_x, trans_min_y,
trans_max_x, trans_max_y);
}
if (use_attribute_query) {
if (poLayer->SetAttributeFilter(attribute_query.c_str()) != OGRERR_NONE) {
SG_LOG( SG_GENERAL, SG_ALERT, "Error in query expression '" << attribute_query << "'" );
exit( 1 );
}
}
// Generate the work queue for this layer
OGRFeature *poFeature;
poLayer->SetNextByIndex(start_record);
while ( ( poFeature = poLayer->GetNextFeature()) != NULL )
{
global_workQueue.push( poFeature );
}
// Now process the workqueue with threads
std::vector<Decoder *> decoders;
for (int i=0; i<num_threads; i++) {
Decoder* decoder = new Decoder( poCT, results );
decoder->start();
decoders.push_back( decoder );
}
// Then wait until they are finished
for (unsigned int i=0; i<decoders.size(); i++) {
decoders[i]->join();
}
OCTDestroyCoordinateTransformation ( poCT );
}
void usage(char* progname) {
SG_LOG( SG_GENERAL, SG_ALERT, "Usage: " <<
progname << " [options...] <work_dir> <datasource> [<layername>...]" );
SG_LOG( SG_GENERAL, SG_ALERT, "Options:" );
SG_LOG( SG_GENERAL, SG_ALERT, "--log-level priority" );
SG_LOG( SG_GENERAL, SG_ALERT, " Width in priority being bulk|debug|info|warn|alert" );
SG_LOG( SG_GENERAL, SG_ALERT, "--continue-on-errors" );
SG_LOG( SG_GENERAL, SG_ALERT, " Continue even if the file seems fishy" );
SG_LOG( SG_GENERAL, SG_ALERT, "--max-segment max_segment_length" );
SG_LOG( SG_GENERAL, SG_ALERT, " Maximum segment length in meters" );
SG_LOG( SG_GENERAL, SG_ALERT, "--start-record record_no" );
SG_LOG( SG_GENERAL, SG_ALERT, " Start processing at the specified record number (first record num=0)" );
SG_LOG( SG_GENERAL, SG_ALERT, "--where attrib_query" );
SG_LOG( SG_GENERAL, SG_ALERT, " Use an attribute query (like SQL WHERE)" );
SG_LOG( SG_GENERAL, SG_ALERT, "--spat xmin ymin xmax ymax" );
SG_LOG( SG_GENERAL, SG_ALERT, " spatial query extents" );
SG_LOG( SG_GENERAL, SG_ALERT, "--threads" );
SG_LOG( SG_GENERAL, SG_ALERT, " Enable multithreading with user specified number of threads" );
SG_LOG( SG_GENERAL, SG_ALERT, "--all-threads" );
SG_LOG( SG_GENERAL, SG_ALERT, " Enable multithreading with all available cpu cores" );
SG_LOG( SG_GENERAL, SG_ALERT, "" );
SG_LOG( SG_GENERAL, SG_ALERT, "<work_dir>" );
SG_LOG( SG_GENERAL, SG_ALERT, " Directory to put the cliff files in" );
SG_LOG( SG_GENERAL, SG_ALERT, "<datasource>" );
SG_LOG( SG_GENERAL, SG_ALERT, " The datasource from which to fetch the data" );
SG_LOG( SG_GENERAL, SG_ALERT, "<layername>..." );
SG_LOG( SG_GENERAL, SG_ALERT, " The layers to process." );
SG_LOG( SG_GENERAL, SG_ALERT, " If no layer is given, all layers in the datasource are used" );
exit(-1);
}
void
setLoggingPriority (const char * p)
{
if (p == 0)
return;
string priority = p;
if (priority == "bulk") {
sglog().set_log_priority(SG_BULK);
} else if (priority == "debug") {
sglog().set_log_priority(SG_DEBUG);
} else if (priority.empty() || priority == "info") { // default
sglog().set_log_priority(SG_INFO);
} else if (priority == "warn") {
sglog().set_log_priority(SG_WARN);
} else if (priority == "alert") {
sglog().set_log_priority(SG_ALERT);
} else {
SG_LOG(SG_GENERAL, SG_WARN, "Unknown logging priority " << priority);
}
}
int main( int argc, char **argv ) {
char* progname=argv[0];
string datasource,work_dir;
auto start_time = std::chrono::high_resolution_clock::now();
sglog().setLogLevels( SG_ALL, SG_INFO );
while (argc>1) {
if (!strcmp(argv[1],"--log-level")) {
if (argc<3) {
usage(progname);
}
setLoggingPriority(argv[2]);
argv+=2;
argc-=2;
} else if (!strcmp(argv[1],"--seperate-segments")) {
argv++;
argc--;
seperate_segments=1;
} else if (!strcmp(argv[1],"--continue-on-errors")) {
argv++;
argc--;
continue_on_errors=1;
} else if (!strcmp(argv[1],"--start-record")) {
if (argc<3) {
usage(progname);
}
start_record=atoi(argv[2]);
argv+=2;
argc-=2;
} else if (!strcmp(argv[1],"--where")) {
if (argc<3) {
usage(progname);
}
use_attribute_query=true;
attribute_query=argv[2];
argv+=2;
argc-=2;
} else if (!strcmp(argv[1],"--spat")) {
if (argc<6) {
usage(progname);
}
use_spatial_query=true;
spat_min_x=atof(argv[2]);
spat_min_y=atof(argv[3]);
spat_max_x=atof(argv[4]);
spat_max_y=atof(argv[5]);
argv+=5;
argc-=5;
} else if (!strcmp(argv[1],"--threads")) {
if (argc<3) {
usage(progname);
}
num_threads=atoi(argv[2]);
argv+=2;
argc-=2;
} else if (!strcmp(argv[1],"--all-threads")) {
num_threads=boost::thread::hardware_concurrency();
argv+=1;
argc-=1;
} else if (!strcmp(argv[1],"--debug")) {
argv++;
argc--;
save_shapefiles=true;
} else if (!strcmp(argv[1],"--help")) {
usage(progname);
} else {
break;
}
}
SG_LOG( SG_GENERAL, SG_ALERT, "\ncliff-decode version " << getTGVersion() );
if (argc<3) {
usage(progname);
}
work_dir=argv[1];
datasource=argv[2];
SGPath sgp( work_dir );
sgp.append( "dummy" );
sgp.create_dir( 0755 );
tgChopper results( work_dir );
SG_LOG( SG_GENERAL, SG_DEBUG, "Opening datasource " << datasource << " for reading." );
GDALAllRegister();
GDALDataset *poDS;
poDS = (GDALDataset*) GDALOpenEx( datasource.c_str(), GDAL_OF_VECTOR, NULL, NULL, NULL );
if( poDS == NULL )
{
SG_LOG( SG_GENERAL, SG_ALERT, "Failed opening datasource " << datasource );
return EXIT_FAILURE;
}
SG_LOG( SG_GENERAL, SG_ALERT, "Processing datasource " << datasource );
OGRLayer *poLayer;
if (argc>3) {
for (int i=3;i<argc;i++) {
poLayer = poDS->GetLayerByName( argv[i] );
if (poLayer == NULL )
{
SG_LOG( SG_GENERAL, SG_ALERT, "Failed opening layer " << argv[i] << " from datasource " << datasource );
return EXIT_FAILURE;
}
processLayer(poLayer, results );
}
} else {
for (int i=0;i<poDS->GetLayerCount();i++) {
poLayer = poDS->GetLayer(i);
assert(poLayer != NULL);
processLayer(poLayer, results );
}
}
GDALClose(poDS);
SG_LOG(SG_GENERAL, SG_ALERT, "Saving to buckets");
results.Add_Extension("cliffs");
results.Save( save_shapefiles );
auto finish_time = std::chrono::high_resolution_clock::now();
std::chrono::duration<double> elapsed = finish_time - start_time;
std::cout << std::endl << "Elapsed time: " << elapsed.count() << " seconds" << std::endl << std::endl;
return 0;
}

14
src/Prep/DemChop/CMakeLists.txt
View file

@ -1,3 +1,4 @@
include_directories(${PROJECT_SOURCE_DIR}/src/Lib/terragear)
add_executable(demchop demchop.cxx)
@ -19,6 +20,17 @@ target_link_libraries(hgtchop
install(TARGETS hgtchop RUNTIME DESTINATION bin)
add_executable(dtedchop dtedchop.cxx)
target_link_libraries(dtedchop
HGT
${ZLIB_LIBRARY}
${SIMGEAR_CORE_LIBRARIES}
${SIMGEAR_CORE_LIBRARY_DEPENDENCIES})
install(TARGETS dtedchop RUNTIME DESTINATION bin)
if(TIFF_FOUND)
if(MSVC AND CMAKE_CL_64)
set( SRTMCHOP_LIBRARIES ${JPEG_LIBRARY} )
@ -38,6 +50,7 @@ endif(TIFF_FOUND)
add_executable(fillvoids fillvoids.cxx)
target_link_libraries(fillvoids
Array
terragear
${ZLIB_LIBRARY}
${SIMGEAR_CORE_LIBRARIES}
${SIMGEAR_CORE_LIBRARY_DEPENDENCIES})
@ -47,6 +60,7 @@ install(TARGETS fillvoids RUNTIME DESTINATION bin)
add_executable(testassem testassem.cxx)
target_link_libraries(testassem
Array
terragear
${ZLIB_LIBRARY}
${SIMGEAR_CORE_LIBRARIES}
${SIMGEAR_CORE_LIBRARY_DEPENDENCIES})

110
src/Prep/DemChop/dtedchop.cxx Normal file
View file

@ -0,0 +1,110 @@
// dtedchop.cxx -- chop up a dted file into it's corresponding pieces and stuff
// them into the workspace directory
//
// Adapted by James Hester from hgtchop
// Written by Curtis Olson, started March 1999.
//
// Copyright (C) 1997 Curtis L. Olson - http://www.flightgear.org/~curt
// Copyright (C) 2018 James Hester
//
// 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., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
//
#ifdef HAVE_CONFIG_H
# include <config.h>
#endif
#include <simgear/compiler.h>
#include <cstdlib>
#include <string>
#include <iostream>
#include <simgear/bucket/newbucket.hxx>
#include <simgear/debug/logstream.hxx>
#include <simgear/misc/sg_path.hxx>
#include <simgear/misc/sg_dir.hxx>
#include <Include/version.h>
#include <HGT/dted.hxx>
using std::cout;
using std::endl;
using std::string;
int main(int argc, char **argv) {
sglog().setLogLevels( SG_ALL, SG_WARN );
SG_LOG( SG_GENERAL, SG_ALERT, "dtedchop version " << getTGVersion() << "\n" );
if ( argc != 4 ) {
cout << "Usage " << argv[0] << " <resolution> <dted_file> <work_dir>" << endl;
cout << endl;
cout << "\tresolution must be either 1 or 3 (1-arc-sec or 3-arc-sec)" << endl;
return EXIT_FAILURE;
}
int resolution = std::stoi(string(argv[1]));
string dted_name = string(argv[2]);
string work_dir = string(argv[3]);
// determine if file is 1arc-sec or 3arc-sec variety
if ( resolution != 1 && resolution != 3 ) {
cout << "ERROR: resolution must be 1 or 3." << endl;
return EXIT_FAILURE;
}
SGPath sgp( work_dir );
simgear::Dir workDir(sgp);
workDir.create(0755);
TGDted dted(resolution, dted_name);
dted.load();
dted.close();
SGGeod min = SGGeod::fromDeg( dted.get_originx() / 3600.0 + SG_HALF_BUCKET_SPAN,
dted.get_originy() / 3600.0 + SG_HALF_BUCKET_SPAN );
SGGeod max = SGGeod::fromDeg( (dted.get_originx() + dted.get_cols() * dted.get_col_step()) / 3600.0 - SG_HALF_BUCKET_SPAN,
(dted.get_originy() + dted.get_rows() * dted.get_row_step()) / 3600.0 - SG_HALF_BUCKET_SPAN );
SGBucket b_min( min );
SGBucket b_max( max );
if ( b_min == b_max ) {
dted.write_area( work_dir, b_min );
} else {
SGBucket b_cur;
int dx, dy;
sgBucketDiff(b_min, b_max, &dx, &dy);
cout << "DTED file spans tile boundaries (ok)" << endl;
cout << " dx = " << dx << " dy = " << dy << endl;
if ( (dx > 20) || (dy > 20) ) {
cout << "somethings really wrong!!!!" << endl;
return EXIT_FAILURE;
}
for ( int j = 0; j <= dy; ++j ) {
for ( int i = 0; i <= dx; ++i ) {
b_cur = b_min.sibling(i, j);
dted.write_area( work_dir, b_cur );
}
}
}
return EXIT_SUCCESS;
}

3
src/Prep/TerraFit/CMakeLists.txt
View file

@ -1,8 +1,11 @@
add_executable(terrafit terrafit.cc)
include_directories(${PROJECT_SOURCE_DIR}/src/Lib/terragear)
target_link_libraries(terrafit
Array Terra
terragear
${ZLIB_LIBRARY}
${SIMGEAR_CORE_LIBRARIES}
${SIMGEAR_CORE_LIBRARY_DEPENDENCIES})

3
src/Prep/TerraFit/terrafit.cc
View file

@ -239,7 +239,8 @@ void walk_path(const SGPath& path) {
return;
}
if ((path.lower_extension() == "arr") || (path.complete_lower_extension() == "arr.gz")) {
if ((path.lower_extension() == "arr") || (path.complete_lower_extension() == "arr.rectified.gz") ||
(path.complete_lower_extension() == "arr.gz")) {
SG_LOG(SG_GENERAL, SG_DEBUG, "will queue " << path);
queue_fit_file(path);
} else if (path.isDir()) {