SHPParser

src/Main/CMakeLists.txt

@@ -139,6 +139,7 @@ target_link_libraries(fgfs
 	${OPENSCENEGRAPH_LIBRARIES}
 	${OPENGL_LIBRARIES}
 	${PLIB_LIBRARIES}
+        ${ZLIB_LIBRARY}
 	${HLA_LIBRARIES}
 	${EVENT_INPUT_LIBRARIES}
 	${SIMGEAR_CORE_LIBRARY_DEPENDENCIES}

src/Navaids/CMakeLists.txt

@@ -18,6 +18,7 @@ set(SOURCES
     NavDataCache.cxx
     PositionedOctree.cxx
     PolyLine.cxx
+    SHPParser.cxx
 	)
 
 set(HEADERS
@@ -38,6 +39,7 @@ set(HEADERS
     NavDataCache.hxx
     PositionedOctree.hxx
     PolyLine.hxx
+    SHPParser.hxx
     CacheSchema.h
     )
 

src/Navaids/PolyLine.cxx

@@ -95,19 +95,27 @@ PolyLineList PolyLine::createChunked(Type aTy, const SGGeodVec& aRawPoints)
     return result;
 }
 
-void PolyLine::addToSpatialIndex() const
+PolyLineRef PolyLine::create(PolyLine::Type aTy, const SGGeodVec &aRawPoints)
 {
-    std::set<Octree::Leaf*> seen;
-    
-    BOOST_FOREACH(const SGGeod& g, m_data) {
-        SGVec3d cart(SGVec3d::fromGeod(g));
-        Octree::Leaf* lf = Octree::global_spatialOctree->findLeafForPos(cart);
-        if (seen.find(lf) != seen.end()) {
-            continue; // don't insert multiple times
+    return new PolyLine(aTy, aRawPoints);
 }
 
-        lf->addPolyLine(const_cast<PolyLine*>(this));
-    } // of data points iteration
+void PolyLine::addToSpatialIndex() const
+{
+    Octree::Node* node = Octree::global_spatialOctree->findNodeForBox(cartesianBox());
+    node->addPolyLine(const_cast<PolyLine*>(this));
+}
+
+SGBoxd PolyLine::cartesianBox() const
+{
+    SGBoxd result;
+    SGGeodVec::const_iterator it;
+    for (it = m_data.begin(); it != m_data.end(); ++it) {
+        SGVec3d cart = SGVec3d::fromGeod(*it);
+        result.expandBy(cart);
+    }
+
+    return result;
 }
 
 class SingleTypeFilter : public PolyLine::TypeFilter

src/Navaids/PolyLine.hxx

@@ -27,6 +27,8 @@
 #include <simgear/sg_inlines.h>
 #include <simgear/structure/SGSharedPtr.hxx>
 #include <simgear/math/SGMath.hxx>
+#include <simgear/math/SGBox.hxx>
+#include <simgear/math/SGGeometryFwd.hxx>
 
 namespace flightgear
 {
@@ -60,6 +62,8 @@ public:
         NATIONAL_BOUNDARY, /// aka a border
         REGIONAL_BOUNDARY, /// state / province / country / department
         RIVER,
+        LAKE,
+        URBAN,
         // airspace types in the future
         LAST_TYPE
     };
@@ -85,6 +89,8 @@ public:
      */
     static PolyLineList createChunked(Type aTy, const SGGeodVec& aRawPoints);
     
+    static PolyLineRef create(Type aTy, const SGGeodVec& aRawPoints);
+
     /**
      * retrieve all the lines within a range of a search point.
      * lines are returned if any point is near the search location.
@@ -98,14 +104,18 @@ public:
     };
     
     static PolyLineList linesNearPos(const SGGeod& aPos, double aRangeNm, const TypeFilter& aFilter);
-private:
+
+    SGBoxd cartesianBox() const;
+
     void addToSpatialIndex() const;
 
+private:
+    
     PolyLine(Type aTy, const SGGeodVec& aPoints);
     
     Type m_type;
     SGGeodVec m_data;
-    // cache the bounding box?
+
 };
     
 

src/Navaids/PositionedOctree.cxx

@@ -47,6 +47,27 @@ namespace Octree
   
 Node* global_spatialOctree = NULL;
 
+
+void Node::addPolyLine(const PolyLineRef& aLine)
+{
+    lines.push_back(aLine);
+}
+
+void Node::visitForLines(const SGVec3d& aPos, double aCutoff,
+                           PolyLineList& aLines,
+                           FindLinesDeque& aQ) const
+{
+    SG_UNUSED(aPos);
+    SG_UNUSED(aCutoff);
+
+    aLines.insert(aLines.end(), lines.begin(), lines.end());
+}
+
+Node *Node::findNodeForBox(const SGBoxd&) const
+{
+    return const_cast<Node*>(this);
+}
+
 Leaf::Leaf(const SGBoxd& aBox, int64_t aIdent) :
   Node(aBox, aIdent),
   childrenLoaded(false)
@@ -114,18 +135,6 @@ void Leaf::loadChildren()
   childrenLoaded = true;
 }
     
-void Leaf::addPolyLine(PolyLineRef aLine)
-{
-    lines.push_back(aLine);
-}
-
-void Leaf::visitForLines(const SGVec3d& aPos, double aCutoff,
-                           PolyLineList& aLines,
-                           FindLinesDeque& aQ) const
-{
-    aLines.insert(aLines.end(), lines.begin(), lines.end());
-}
-    
 ///////////////////////////////////////////////////////////////////////////////
     
 Branch::Branch(const SGBoxd& aBox, int64_t aIdent) :
@@ -158,6 +167,9 @@ void Branch::visitForLines(const SGVec3d& aPos, double aCutoff,
                            PolyLineList& aLines,
                            FindLinesDeque& aQ) const
 {
+    // add our own lines, easy
+    Node::visitForLines(aPos, aCutoff, aLines, aQ);
+
     for (unsigned int i=0; i<8; ++i) {
         if (!children[i]) {
             continue;
@@ -172,6 +184,35 @@ void Branch::visitForLines(const SGVec3d& aPos, double aCutoff,
     } // of child iteration
 }
 
+static bool boxContainsBox(const SGBoxd& a, const SGBoxd& b)
+{
+    const SGVec3d aMin(a.getMin()),
+            aMax(a.getMax()),
+            bMin(b.getMin()),
+            bMax(b.getMax());
+    for (int i=0; i<3; ++i) {
+        if ((bMin[i] < aMin[i]) || (bMax[i] > aMax[i])) return false;
+    }
+
+    return true;
+}
+
+Node *Branch::findNodeForBox(const SGBoxd &box) const
+{
+    // do this so childAtIndex sees consistent state of
+    // children[] and loaded flag.
+    loadChildren();
+
+    for (unsigned int i=0; i<8; ++i) {
+        const SGBoxd childBox(boxForChild(i));
+        if (boxContainsBox(childBox, box)) {
+            return childAtIndex(i)->findNodeForBox(box);
+        }
+    }
+
+    return Node::findNodeForBox(box);
+}
+
 Node* Branch::childForPos(const SGVec3d& aCart) const
 {
   assert(contains(aCart));

src/Navaids/PositionedOctree.hxx

@@ -152,8 +152,13 @@ namespace Octree
       
     virtual void visitForLines(const SGVec3d& aPos, double aCutoff,
                          PolyLineList& aLines,
-                         FindLinesDeque& aQ) const = 0;
+                         FindLinesDeque& aQ) const;
+
+    virtual Node* findNodeForBox(const SGBoxd& box) const;
+
     virtual ~Node() {}
+
+    void addPolyLine(const PolyLineRef&);
   protected:
     Node(const SGBoxd &aBox, int64_t aIdent) :
     _ident(aIdent),
@@ -163,6 +168,8 @@ namespace Octree
     
     const int64_t _ident;
     const SGBoxd _box;
+
+    PolyLineList lines;
   };
   
   class Leaf : public Node
@@ -181,19 +188,12 @@ namespace Octree
     
     void insertChild(FGPositioned::Type ty, PositionedID id);
 
-    void addPolyLine(PolyLineRef);
-    
-    virtual void visitForLines(const SGVec3d& aPos, double aCutoff,
-                               PolyLineList& aLines,
-                               FindLinesDeque& aQ) const;
   private:
     bool childrenLoaded;
     
     typedef std::multimap<FGPositioned::Type, PositionedID> ChildMap;
     ChildMap children;
             
-    PolyLineList lines;
-      
     void loadChildren();
   };
   
@@ -217,6 +217,9 @@ namespace Octree
     virtual void visitForLines(const SGVec3d& aPos, double aCutoff,
                                PolyLineList& aLines,
                                FindLinesDeque& aQ) const;
+
+    virtual Node* findNodeForBox(const SGBoxd& box) const;
+
   private:
     Node* childForPos(const SGVec3d& aCart) const;
     Node* childAtIndex(int childIndex) const;

src/Navaids/SHPParser.cxx

@@ -0,0 +1,190 @@
+/**
+ * SHPParser - parse ESRI ShapeFiles containing PolyLines */
+
+// Written by James Turner, started 2013.
+//
+// Copyright (C) 2013 James Turner <zakalawe@mac.com>
+//
+// 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 Street, Fifth Floor, Boston, MA  02110-1301, USA.
+
+#ifdef HAVE_CONFIG_H
+     #include "config.h"
+#endif
+
+#include "SHPParser.hxx"
+
+#include <simgear/debug/logstream.hxx>
+#include <simgear/structure/exception.hxx>
+#include <simgear/misc/sg_path.hxx>
+#include <simgear/io/lowlevel.hxx>
+
+// http://www.esri.com/library/whitepapers/pdfs/shapefile.pdf table 1
+const int SHP_FILE_MAGIC = 9994;
+const int SHP_FILE_VERSION = 1000;
+
+const int SHP_NULL_TYPE = 0;
+const int SHP_POLYLINE_TYPE = 3;
+const int SHP_POLYGON_TYPE = 5;
+
+namespace
+{
+
+void sgReadIntBE ( gzFile fd, int& var )
+{
+    if ( gzread ( fd, &var, sizeof(int) ) != sizeof(int) ) {
+        throw sg_io_exception("gzread failed");
+    }
+
+    if ( sgIsLittleEndian() ) {
+        sgEndianSwap( (uint32_t *) &var);
+    }
+}
+
+void sgReadIntLE ( gzFile fd, int& var )
+{
+    if ( gzread ( fd, &var, sizeof(int) ) != sizeof(int) ) {
+        throw sg_io_exception("gzread failed");
+    }
+
+    if ( sgIsBigEndian() ) {
+        sgEndianSwap( (uint32_t *) &var);
+    }
+}
+
+
+void readSHPRecordHeader(gzFile fd, int &recordNumber, int& contentLength)
+{
+    sgReadIntBE(fd, recordNumber);
+    sgReadIntBE(fd, contentLength);
+}
+
+void parseSHPPoints2D(gzFile fd, int numPoints, flightgear::SGGeodVec& aPoints)
+{
+    aPoints.reserve(numPoints);
+    std::vector<double> ds;
+    ds.resize(numPoints * 2);
+    sgReadDouble(fd, numPoints * 2, ds.data());
+
+    unsigned int index = 0;
+    for (int i=0; i<numPoints; ++i, index += 2) {
+        aPoints.push_back(SGGeod::fromDeg(ds[index], ds[index+1]));
+    }
+}
+
+} // anonymous namespace
+
+namespace flightgear
+{
+
+void SHPParser::parsePolyLines(const SGPath& aPath, PolyLine::Type aTy,
+                               PolyLineList& aResult, bool aClosed)
+{
+    gzFile file = gzopen(aPath.c_str(), "rb");
+    if (!file) {
+        throw sg_io_exception("couldn't open file:" + aPath.str());
+    }
+
+    try {
+        int header, fileLength, fileVersion, shapeType;
+        sgReadIntBE(file, header);
+        if (header != SHP_FILE_MAGIC) {
+            throw sg_io_exception("bad SHP header value", aPath);
+        }
+
+        // skip 5 ints, then read the file length
+        for (int i=0; i<6; ++i) {
+            sgReadIntBE(file, fileLength);
+        }
+
+        sgReadIntLE(file, fileVersion);
+        sgReadIntLE(file, shapeType);
+
+        if (fileVersion != SHP_FILE_VERSION) {
+            throw sg_io_exception("bad SHP file version", aPath);
+        }
+
+        if (aClosed && (shapeType != SHP_POLYGON_TYPE)) {
+            throw sg_io_exception("SHP file does not contain Polygon data", aPath);
+        }
+
+       if (!aClosed && (shapeType != SHP_POLYLINE_TYPE)) {
+            throw sg_io_exception("SHP file does not contain PolyLine data", aPath);
+        }
+
+        // we don't care about range values
+        double range;
+        for (int i=0; i<8; ++i) {
+            sgReadDouble(file, &range);
+        }
+
+        int readLength = 100; // sizeof the header
+        while (readLength < fileLength) {
+            int recordNumber, contentLengthWords;
+            readSHPRecordHeader(file, recordNumber, contentLengthWords);
+
+            int recordShapeType;
+            sgReadIntLE(file, recordShapeType);
+            if (recordShapeType == SHP_NULL_TYPE) {
+                continue; // nothing else to do
+            }
+
+            if (recordShapeType != shapeType) {
+                // vesion 1000 requires files to have homogenous shape type
+                throw sg_io_exception("SHP file shape-type mismatch", aPath);
+            }
+        // read PolyLine record from now on
+            double box[4];
+            for (int i=0; i<4; ++i) {
+                sgReadDouble(file, &box[i]);
+            }
+
+            int numParts, numPoints;
+            sgReadInt(file, &numParts);
+            sgReadInt(file, &numPoints);
+
+            std::vector<int> parts;
+            parts.resize(numParts);
+            sgReadInt(file, numParts, parts.data());
+
+            SGGeodVec points;
+            parseSHPPoints2D(file, numPoints, points);
+
+            for (int part=0; part<numParts; ++part) {
+                SGGeodVec partPoints;
+                unsigned int startIndex = parts[part];
+                unsigned int endIndex = ((part + 1) == numParts) ? numPoints : parts[part + 1];
+                partPoints.insert(partPoints.begin(), points.begin() + startIndex, points.begin() + endIndex);
+
+                if (aClosed) {
+                    aResult.push_back(PolyLine::create(aTy, partPoints));
+                } else {
+                    PolyLineList lines = PolyLine::createChunked(aTy, partPoints);
+                    aResult.insert(aResult.end(), lines.begin(), lines.end());
+                }
+            }
+
+            // total record size if contentLenght + 4 words for the two record fields
+            readLength += (contentLengthWords + 4);
+        // partition
+        } // of record parsing
+
+    } catch (sg_exception& e) {
+        aResult.clear();
+        gzclose(file);
+        throw e; // rethrow
+    }
+}
+
+} // of namespace flightgear

src/Navaids/SHPParser.hxx

@@ -0,0 +1,46 @@
+/**
+ * SHPParser - parse ESRI ShapeFiles containing PolyLines */
+
+// Written by James Turner, started 2013.
+//
+// Copyright (C) 2013 James Turner <zakalawe@mac.com>
+//
+// 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 Street, Fifth Floor, Boston, MA  02110-1301, USA.
+
+#ifndef FG_SHP_PARSER_HXX
+#define FG_SHP_PARSER_HXX
+
+#include <Navaids/PolyLine.hxx>
+
+// forward decls
+class SGPath;
+
+namespace flightgear
+{
+
+class SHPParser
+{
+public:
+    /**
+     * Parse a shape file containing PolyLine data.
+     *
+     * Throws sg_exceptions if parsing problems occur.
+     */
+    static void parsePolyLines(const SGPath&, PolyLine::Type aTy, PolyLineList& aResult, bool aClosed);
+};
+
+} // of namespace flightgear
+
+#endif // of FG_SHP_PARSER_HXX