src/Navaids/SHPParser.cxx

   1 /**
   2  * SHPParser - parse ESRI ShapeFiles containing PolyLines */
   3
   4 // Written by James Turner, started 2013.
   5 //
   6 // Copyright (C) 2013 James Turner <zakalawe@mac.com>
   7 //
   8 // This program is free software; you can redistribute it and/or
   9 // modify it under the terms of the GNU General Public License as
  10 // published by the Free Software Foundation; either version 2 of the
  11 // License, or (at your option) any later version.
  12 //
  13 // This program is distributed in the hope that it will be useful, but
  14 // WITHOUT ANY WARRANTY; without even the implied warranty of
  15 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  16 // General Public License for more details.
  17 //
  18 // You should have received a copy of the GNU General Public License
  19 // along with this program; if not, write to the Free Software
  20 // Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
  21
  22 #ifdef HAVE_CONFIG_H
  23      #include "config.h"
  24 #endif
  25
  26 #include "SHPParser.hxx"
  27
  28 #include <simgear/debug/logstream.hxx>
  29 #include <simgear/structure/exception.hxx>
  30 #include <simgear/misc/sg_path.hxx>
  31 #include <simgear/io/lowlevel.hxx>
  32
  33 // http://www.esri.com/library/whitepapers/pdfs/shapefile.pdf table 1
  34 const int SHP_FILE_MAGIC = 9994;
  35 const int SHP_FILE_VERSION = 1000;
  36
  37 const int SHP_NULL_TYPE = 0;
  38 const int SHP_POLYLINE_TYPE = 3;
  39 const int SHP_POLYGON_TYPE = 5;
  40
  41 namespace
  42 {
  43
  44 void sgReadIntBE ( gzFile fd, int& var )
  45 {
  46     if ( gzread ( fd, &var, sizeof(int) ) != sizeof(int) ) {
  47         throw sg_io_exception("gzread failed");
  48     }
  49
  50     if ( sgIsLittleEndian() ) {
  51         sgEndianSwap( (uint32_t *) &var);
  52     }
  53 }
  54
  55 void sgReadIntLE ( gzFile fd, int& var )
  56 {
  57     if ( gzread ( fd, &var, sizeof(int) ) != sizeof(int) ) {
  58         throw sg_io_exception("gzread failed");
  59     }
  60
  61     if ( sgIsBigEndian() ) {
  62         sgEndianSwap( (uint32_t *) &var);
  63     }
  64 }
  65
  66
  67 void readSHPRecordHeader(gzFile fd, int &recordNumber, int& contentLength)
  68 {
  69     sgReadIntBE(fd, recordNumber);
  70     sgReadIntBE(fd, contentLength);
  71 }
  72
  73 void parseSHPPoints2D(gzFile fd, int numPoints, flightgear::SGGeodVec& aPoints)
  74 {
  75     aPoints.reserve(numPoints);
  76     std::vector<double> ds;
  77     ds.resize(numPoints * 2);
  78     sgReadDouble(fd, numPoints * 2, ds.data());
  79
  80     unsigned int index = 0;
  81     for (int i=0; i<numPoints; ++i, index += 2) {
  82         aPoints.push_back(SGGeod::fromDeg(ds[index], ds[index+1]));
  83     }
  84 }
  85
  86 } // anonymous namespace
  87
  88 namespace flightgear
  89 {
  90
  91 void SHPParser::parsePolyLines(const SGPath& aPath, PolyLine::Type aTy,
  92                                PolyLineList& aResult, bool aClosed)
  93 {
  94     gzFile file = gzopen(aPath.c_str(), "rb");
  95     if (!file) {
  96         throw sg_io_exception("couldn't open file:" + aPath.str());
  97     }
  98
  99     try {
 100         int header, fileLength, fileVersion, shapeType;
 101         sgReadIntBE(file, header);
 102         if (header != SHP_FILE_MAGIC) {
 103             throw sg_io_exception("bad SHP header value", aPath);
 104         }
 105
 106         // skip 5 ints, then read the file length
 107         for (int i=0; i<6; ++i) {
 108             sgReadIntBE(file, fileLength);
 109         }
 110
 111         sgReadIntLE(file, fileVersion);
 112         sgReadIntLE(file, shapeType);
 113
 114         if (fileVersion != SHP_FILE_VERSION) {
 115             throw sg_io_exception("bad SHP file version", aPath);
 116         }
 117
 118         if (aClosed && (shapeType != SHP_POLYGON_TYPE)) {
 119             throw sg_io_exception("SHP file does not contain Polygon data", aPath);
 120         }
 121
 122        if (!aClosed && (shapeType != SHP_POLYLINE_TYPE)) {
 123             throw sg_io_exception("SHP file does not contain PolyLine data", aPath);
 124         }
 125
 126         // we don't care about range values
 127         double range;
 128         for (int i=0; i<8; ++i) {
 129             sgReadDouble(file, &range);
 130         }
 131
 132         int readLength = 100; // sizeof the header
 133         while (readLength < fileLength) {
 134             int recordNumber, contentLengthWords;
 135             readSHPRecordHeader(file, recordNumber, contentLengthWords);
 136
 137             int recordShapeType;
 138             sgReadIntLE(file, recordShapeType);
 139             if (recordShapeType == SHP_NULL_TYPE) {
 140                 continue; // nothing else to do
 141             }
 142
 143             if (recordShapeType != shapeType) {
 144                 // vesion 1000 requires files to have homogenous shape type
 145                 throw sg_io_exception("SHP file shape-type mismatch", aPath);
 146             }
 147         // read PolyLine record from now on
 148             double box[4];
 149             for (int i=0; i<4; ++i) {
 150                 sgReadDouble(file, &box[i]);
 151             }
 152
 153             int numParts, numPoints;
 154             sgReadInt(file, &numParts);
 155             sgReadInt(file, &numPoints);
 156
 157             std::vector<int> parts;
 158             parts.resize(numParts);
 159             sgReadInt(file, numParts, parts.data());
 160
 161             SGGeodVec points;
 162             parseSHPPoints2D(file, numPoints, points);
 163
 164             for (int part=0; part<numParts; ++part) {
 165                 SGGeodVec partPoints;
 166                 unsigned int startIndex = parts[part];
 167                 unsigned int endIndex = ((part + 1) == numParts) ? numPoints : parts[part + 1];
 168                 partPoints.insert(partPoints.begin(), points.begin() + startIndex, points.begin() + endIndex);
 169
 170                 if (aClosed) {
 171                     aResult.push_back(PolyLine::create(aTy, partPoints));
 172                 } else {
 173                     PolyLineList lines = PolyLine::createChunked(aTy, partPoints);
 174                     aResult.insert(aResult.end(), lines.begin(), lines.end());
 175                 }
 176             }
 177
 178             // total record size if contentLenght + 4 words for the two record fields
 179             readLength += (contentLengthWords + 4);
 180         // partition
 181         } // of record parsing
 182
 183     } catch (sg_exception& e) {
 184         aResult.clear();
 185         gzclose(file);
 186         throw e; // rethrow
 187     }
 188 }
 189
 190 } // of namespace flightgear