Merge branch 'cliffs' of ../../stable/terragear into next

src/Prep/Cliff/CMakeLists.txt

@@ -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)

src/Prep/Cliff/cliff-decode.cxx

@@ -0,0 +1,434 @@
+// cliff-decode.cxx -- process OGR layers and extract lines
+//                   representing discontinuities, clipping
+//                   against and sorting
+//                   them into the relevant tiles.
+//
+// Written by James Hester starting 2018
+//
+// Based on 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 <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. Do not use usual clipper, as we wish to preserve the
+    // whole contour, even if it goes outside the bounding box, as height
+    // calculations along the boundary might be affected by it.
+    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;
+
+    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 );
+
+    // initialize persistant polygon counter
+    //string counter_file = work_dir + "/poly_counter";
+    //poly_index_init( counter_file );
+
+    // new chop api
+    //std::string counter_file2 = work_dir + "/poly_counter2";
+    //tgPolygon::ChopIdxInit( counter_file );
+
+    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 );
+        exit( 1 );
+    }
+
+    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 );
+                exit( 1 );
+            }
+            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 );
+    return 0;
+}
+
+