--- /dev/null
+/*
+ * Copyright (c) 2018, Bertold Van den Bergh (vandenbergh@bertold.org)
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions are met:
+ * * Redistributions of source code must retain the above copyright
+ * notice, this list of conditions and the following disclaimer.
+ * * Redistributions in binary form must reproduce the above copyright
+ * notice, this list of conditions and the following disclaimer in the
+ * documentation and/or other materials provided with the distribution.
+ * * Neither the name of the author nor the
+ * names of its contributors may be used to endorse or promote products
+ * derived from this software without specific prior written permission.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+ * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+ * DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR DISTRIBUTOR BE LIABLE FOR ANY
+ * DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+ * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+ * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
+ * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+ * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+ * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ */
+
+#include <shapefil.h>
+#include <iostream>
+#include <limits>
+#include <fstream>
+#include <vector>
+#include <algorithm>
+#include <unordered_map>
+#include <functional>
+#include <math.h>
+#include <tuple>
+
+const double Inf = std::numeric_limits<float>::infinity();
+
+std::unordered_map<std::string, std::string> alpha2ToName;
+std::unordered_map<std::string, std::string> tzidToAlpha2;
+
+void errorFatal(std::string what)
+{
+ std::cerr<<what<<"\n";
+ exit(1);
+}
+
+int encodeVariableLength(std::vector<uint8_t>& output, int64_t valueIn, bool handleNeg = true)
+{
+ uint64_t value = valueIn * 2;
+ if(valueIn < 0) {
+ value = -valueIn * 2 + 1;
+ }
+
+ if(!handleNeg) {
+ value = valueIn;
+ }
+
+ int bytesUsed = 0;
+ do {
+ uint8_t byteOut = value & 0x7F;
+ if(value >= 128) {
+ byteOut |= 0x80;
+ }
+ output.push_back(byteOut);
+ bytesUsed ++;
+ value >>= 7;
+ } while(value);
+
+ return bytesUsed;
+}
+
+int64_t doubleToFixedPoint(double input, double scale, unsigned int precision = 32)
+{
+ double inputScaled = input / scale;
+ return inputScaled * pow(2, precision-1);
+
+}
+
+struct Point {
+ Point(double lat = 0, double lon = 0)
+ {
+ lat_ = lat;
+ lon_ = lon;
+ }
+
+ Point operator-(const Point& p)
+ {
+ Point result(lat_ - p.lat_, lon_ - p.lon_);
+ return result;
+ }
+
+ std::tuple<int64_t, int64_t> toFixedPoint(unsigned int precision = 32)
+ {
+ int64_t latFixedPoint = doubleToFixedPoint(lat_, 90, precision);
+ int64_t lonFixedPoint = doubleToFixedPoint(lon_, 180, precision);
+
+ return std::make_tuple(latFixedPoint, lonFixedPoint);
+ }
+
+ int encodePointBinary(std::vector<uint8_t>& output, unsigned int precision = 32)
+ {
+ int64_t latFixedPoint, lonFixedPoint;
+ std::tie(latFixedPoint, lonFixedPoint) = toFixedPoint(precision);
+
+ int bytesUsed = encodeVariableLength(output, latFixedPoint);
+ bytesUsed += encodeVariableLength(output, lonFixedPoint);
+
+ return bytesUsed;
+ }
+
+ double lat_;
+ double lon_;
+};
+
+struct PolygonData {
+ Point boundingMin;
+ Point boundingMax;
+ std::vector<Point> points_;
+ unsigned long fileIndex_ = 0;
+ unsigned long metadataId_;
+
+ void processPoint(const Point& p)
+ {
+ if(p.lat_ < boundingMin.lat_) {
+ boundingMin.lat_ = p.lat_;
+ }
+ if(p.lon_ < boundingMin.lon_) {
+ boundingMin.lon_ = p.lon_;
+ }
+ if(p.lat_ > boundingMax.lat_) {
+ boundingMax.lat_ = p.lat_;
+ }
+ if(p.lon_ > boundingMax.lon_) {
+ boundingMax.lon_ = p.lon_;
+ }
+
+ points_.push_back(p);
+ }
+
+ PolygonData(unsigned long id):
+ boundingMin(Inf, Inf),
+ boundingMax(-Inf, -Inf),
+ metadataId_(id)
+ {
+ }
+
+ long encodeBinaryData(std::vector<uint8_t>& output, unsigned int precision = 20)
+ {
+ long bytesEncoded = 0;
+ bool first = true;
+ int64_t latFixedPoint = 0, lonFixedPoint = 0;
+ int64_t latFixedPointPrev, lonFixedPointPrev;
+ uint64_t vertices = 0;
+
+ std::vector<uint8_t> tmp;
+
+ int64_t diffLatAcc = 0, diffLonAcc = 0, diffLatPrev = 0, diffLonPrev = 0;
+
+ for(Point& point: points_) {
+ /* The points should first be rounded, and then the integer value is differentiated */
+ latFixedPointPrev = latFixedPoint;
+ lonFixedPointPrev = lonFixedPoint;
+ std::tie(latFixedPoint, lonFixedPoint) = point.toFixedPoint(precision);
+
+ int64_t diffLat = latFixedPoint - latFixedPointPrev;
+ int64_t diffLon = lonFixedPoint - lonFixedPointPrev;
+
+ if(first) {
+ /* First point is always encoded */
+ vertices++;
+ encodeVariableLength(tmp, latFixedPoint);
+ encodeVariableLength(tmp, lonFixedPoint);
+ first = false;
+ } else {
+ /* Ignore points that are not different */
+ if(!diffLon && !diffLat) {
+ continue;
+ }
+
+ if(diffLat != diffLatPrev || diffLon != diffLonPrev) {
+ /* Encode accumulator */
+ vertices++;
+ encodeVariableLength(tmp, diffLatAcc);
+ encodeVariableLength(tmp, diffLonAcc);
+
+ diffLatAcc = 0;
+ diffLonAcc = 0;
+ }
+
+ diffLatAcc += diffLat;
+ diffLonAcc += diffLon;
+ }
+
+ diffLatPrev = diffLat;
+ diffLonPrev = diffLon;
+ }
+
+ /* Encode final point */
+ vertices++;
+ encodeVariableLength(tmp, diffLatAcc);
+ encodeVariableLength(tmp, diffLonAcc);
+
+ encodeVariableLength(output, vertices, false);
+ std::copy(tmp.begin(), tmp.end(), std::back_inserter(output));
+
+ return bytesEncoded;
+ }
+};
+
+void encodeStringToBinary(std::vector<uint8_t>& output, std::string& input)
+{
+ encodeVariableLength(output, input.size(), false);
+ for(unsigned int i=0; i<input.size(); i++) {
+ output.push_back(input[i] ^ 0x80);
+ }
+}
+
+
+std::unordered_map<std::string, uint64_t> usedStrings_;
+
+struct MetaData {
+ void encodeBinaryData(std::vector<uint8_t>& output)
+ {
+ for(std::string& str: data_) {
+ if(str.length() >= 256) {
+ std::cout << "Metadata string is too long\n";
+ exit(1);
+ }
+
+ if(!usedStrings_.count(str)) {
+ usedStrings_[str] = output.size();
+ encodeStringToBinary(output, str);
+ } else {
+ encodeVariableLength(output, usedStrings_[str] + 256, false);
+ }
+ }
+ }
+
+ std::vector<std::string> data_;
+
+ unsigned long fileIndex_;
+};
+
+
+std::vector<PolygonData*> polygons_;
+std::vector<MetaData> metadata_;
+std::vector<std::string> fieldNames_;
+
+
+unsigned int decodeVariableLength(uint8_t* buffer, int64_t* result, bool handleNeg = true)
+{
+ int64_t value = 0;
+ unsigned int i=0, shift = 0;
+
+ do {
+ value |= (buffer[i] & 0x7F) << shift;
+ shift += 7;
+ } while(buffer[i++] & 0x80);
+
+ if(!handleNeg) {
+ *result = value;
+ } else {
+ *result = (value & 1)?-(value/2):(value/2);
+ }
+ return i;
+}
+
+void readMetaDataTimezone(DBFHandle dataHandle)
+{
+ /* Specify field names */
+ fieldNames_.push_back("TimezoneIdPrefix");
+ fieldNames_.push_back("TimezoneId");
+ fieldNames_.push_back("CountryAlpha2");
+ fieldNames_.push_back("CountryName");
+
+ /* Parse attribute names */
+ for(int i = 0; i < DBFGetRecordCount(dataHandle); i++) {
+ metadata_[i].data_.resize(4);
+ for(int j = 0; j < DBFGetFieldCount(dataHandle); j++) {
+ char fieldTitle[12];
+ int fieldWidth, fieldDecimals;
+ DBFFieldType eType = DBFGetFieldInfo(dataHandle, j, fieldTitle, &fieldWidth, &fieldDecimals);
+
+ fieldTitle[11] = 0;
+ std::string fieldTitleStr(fieldTitle);
+
+ if( eType == FTString ) {
+ if(fieldTitleStr == "tzid") {
+ std::string data = DBFReadStringAttribute(dataHandle, i, j);
+ size_t pos = data.find('/');
+ if (pos == std::string::npos) {
+ metadata_[i].data_.at(0) = data;
+ } else {
+ metadata_[i].data_.at(0) = data.substr(0, pos) + "/";
+ metadata_[i].data_.at(1) = data.substr(pos + 1, std::string::npos);
+ }
+ if(tzidToAlpha2.count(data)) {
+ metadata_[i].data_.at(2) = tzidToAlpha2[data];
+ if(alpha2ToName.count(metadata_[i].data_.at(2))) {
+ metadata_[i].data_.at(3) = alpha2ToName[metadata_[i].data_.at(2)];
+ } else {
+ std::cout<<metadata_[i].data_.at(2)<< " not found in alpha2ToName! ("<<data<<")\n";
+ }
+ } else {
+ std::cout<<data<<" not found in zoneToAlpha2!\n";
+ }
+ }
+ }
+ }
+ }
+}
+
+void readMetaDataNaturalEarthCountry(DBFHandle dataHandle)
+{
+ /* Specify field names */
+ fieldNames_.push_back("Alpha2");
+ fieldNames_.push_back("Alpha3");
+ fieldNames_.push_back("Name");
+
+ /* Parse attribute names */
+ for(int i = 0; i < DBFGetRecordCount(dataHandle); i++) {
+ metadata_[i].data_.resize(3);
+ for(int j = 0; j < DBFGetFieldCount(dataHandle); j++) {
+ char fieldTitle[12];
+ int fieldWidth, fieldDecimals;
+ DBFFieldType eType = DBFGetFieldInfo(dataHandle, j, fieldTitle, &fieldWidth, &fieldDecimals);
+
+ fieldTitle[11] = 0;
+ std::string fieldTitleStr(fieldTitle);
+
+ if( eType == FTString ) {
+ if(fieldTitleStr == "ISO_A2" || fieldTitleStr == "WB_A2") {
+ std::string tmp = DBFReadStringAttribute(dataHandle, i, j);
+ if(tmp != "-99") {
+ metadata_[i].data_.at(0) = tmp;
+ }
+ } else if(fieldTitleStr == "ISO_A3" || fieldTitleStr == "WB_A3" || fieldTitleStr == "BRK_A3") {
+ std::string tmp = DBFReadStringAttribute(dataHandle, i, j);
+ if(tmp != "-99") {
+ metadata_[i].data_.at(1) = tmp;
+ }
+ } else if(fieldTitleStr == "NAME_LONG") {
+ metadata_[i].data_.at(2) = DBFReadStringAttribute(dataHandle, i, j);
+ }
+ }
+
+ }
+ }
+}
+
+std::unordered_map<std::string, std::string> parseAlpha2ToName(DBFHandle dataHandle)
+{
+ std::unordered_map<std::string, std::string> result;
+
+ for(int i = 0; i < DBFGetRecordCount(dataHandle); i++) {
+ std::string alpha2, name;
+ for(int j = 0; j < DBFGetFieldCount(dataHandle); j++) {
+ char fieldTitle[12];
+ int fieldWidth, fieldDecimals;
+ DBFFieldType eType = DBFGetFieldInfo(dataHandle, j, fieldTitle, &fieldWidth, &fieldDecimals);
+
+ fieldTitle[11] = 0;
+ std::string fieldTitleStr(fieldTitle);
+
+ if( eType == FTString ) {
+ if(fieldTitleStr == "ISO_A2" || fieldTitleStr == "WB_A2") {
+ std::string tmp = DBFReadStringAttribute(dataHandle, i, j);
+ if(tmp != "-99" && alpha2 == "") {
+ alpha2 = tmp;
+ }
+ } else if(fieldTitleStr == "NAME_LONG") {
+ name = DBFReadStringAttribute(dataHandle, i, j);
+ }
+ }
+ }
+ if(alpha2 != "") {
+ result[alpha2]=name;
+ }
+ }
+
+ result["GF"]="French Guiana";
+ result["GP"]="Guadeloupe";
+ result["BQ"]="Bonaire";
+ result["MQ"]="Martinique";
+ result["SJ"]="Svalbard and Jan Mayen Islands";
+ result["NO"]="Norway";
+ result["CX"]="Christmas Island";
+ result["CC"]="Cocos Islands";
+ result["YT"]="Mayotte";
+ result["RE"]="Réunion";
+ result["TK"]="Tokelau";
+
+ return result;
+}
+
+std::unordered_map<std::string, std::string> parseTimezoneToAlpha2(std::string path)
+{
+ std::unordered_map<std::string, std::string> result;
+ //TODO: Clean solution...
+#include "zoneToAlpha.h"
+
+ return result;
+}
+
+int main(int argc, char ** argv )
+{
+ if(argc != 6) {
+ std::cout << "Wrong number of parameters\n";
+ return 1;
+ }
+
+ tzidToAlpha2 = parseTimezoneToAlpha2("TODO");
+
+ char tableType = argv[1][0];
+ std::string path = argv[2];
+ std::string outPath = argv[3];
+ unsigned int precision = strtol(argv[4], NULL, 10);
+ std::string notice = argv[5];
+
+ DBFHandle dataHandle = DBFOpen("naturalearth/ne_10m_admin_0_countries_lakes", "rb" );
+ alpha2ToName = parseAlpha2ToName(dataHandle);
+ DBFClose(dataHandle);
+
+ dataHandle = DBFOpen(path.c_str(), "rb" );
+ if( dataHandle == NULL ) {
+ errorFatal("Could not open attribute file\n");
+ }
+
+ metadata_.resize(DBFGetRecordCount(dataHandle));
+ std::cout << "Reading "<<metadata_.size()<<" metadata records.\n";
+
+ if(tableType == 'C') {
+ readMetaDataNaturalEarthCountry(dataHandle);
+ } else if(tableType == 'T') {
+ readMetaDataTimezone(dataHandle);
+ } else {
+ std::cout << "Unknown table type\n";
+ return 1;
+ }
+
+ DBFClose(dataHandle);
+
+ SHPHandle shapeHandle = SHPOpen(path.c_str(), "rb");
+ if( shapeHandle == NULL ) {
+ errorFatal("Could not open shapefile\n");
+ }
+
+ int numEntities, shapeType, totalPolygons = 0;
+ SHPGetInfo(shapeHandle, &numEntities, &shapeType, NULL, NULL);
+
+ std::cout<<"Opened "<<SHPTypeName( shapeType )<< " file with "<<numEntities<<" entries.\n";
+
+ for(int i = 0; i < numEntities; i++ ) {
+ SHPObject *shapeObject;
+
+ shapeObject = SHPReadObject( shapeHandle, i );
+ if(shapeObject) {
+ if(shapeObject->nSHPType != 3 && shapeObject->nSHPType != 5 &&
+ shapeObject->nSHPType != 13 && shapeObject->nSHPType != 15) {
+ std::cout<<"Unsupported shape object ("<< SHPTypeName(shapeObject->nSHPType) <<")\n";
+ continue;
+ }
+
+ int partIndex = 0;
+
+ PolygonData* polygonData = nullptr;
+
+ for(int j = 0; j < shapeObject->nVertices; j++ ) {
+ if(j == 0 || j == shapeObject->panPartStart[partIndex]) {
+ totalPolygons++;
+
+ if(polygonData) {
+ /* Commit it */
+ polygons_.push_back(polygonData);
+ }
+ polygonData = new PolygonData(i);
+
+ if(partIndex + 1 < shapeObject->nParts) {
+ partIndex++;
+ }
+ }
+
+ Point p(shapeObject->padfY[j], shapeObject->padfX[j]);
+ polygonData->processPoint(p);
+
+ }
+
+ if(polygonData) {
+ /* Commit it */
+ polygons_.push_back(polygonData);
+ }
+
+ SHPDestroyObject(shapeObject);
+ }
+ }
+
+ SHPClose(shapeHandle);
+
+ std::cout<<"Parsed "<<totalPolygons<<" polygons.\n";
+
+ /* Sort according to bounding box */
+ std::sort(polygons_.begin(), polygons_.end(), [](PolygonData* a, PolygonData* b) {
+ return a->boundingMin.lat_ < b->boundingMin.lat_;
+ });
+
+ /* Encode data section and store pointers */
+ std::vector<uint8_t> outputData;
+ for(PolygonData* polygon: polygons_) {
+ polygon->fileIndex_ = outputData.size();
+ polygon->encodeBinaryData(outputData, precision);
+ }
+ std::cout << "Encoded data section into "<<outputData.size()<<" bytes.\n";
+
+ /* Encode metadata */
+ std::vector<uint8_t> outputMeta;
+ for(MetaData& metadata: metadata_) {
+ metadata.fileIndex_ = outputMeta.size();
+ metadata.encodeBinaryData(outputMeta);
+ }
+ std::cout << "Encoded metadata into "<<outputMeta.size()<<" bytes.\n";
+
+ /* Encode bounding boxes */
+ std::vector<uint8_t> outputBBox;
+ int64_t prevFileIndex = 0;
+ int64_t prevMetaIndex = 0;
+ for(PolygonData* polygon: polygons_) {
+ polygon->boundingMin.encodePointBinary(outputBBox, precision);
+ polygon->boundingMax.encodePointBinary(outputBBox, precision);
+
+ encodeVariableLength(outputBBox, metadata_.at(polygon->metadataId_).fileIndex_ - prevMetaIndex);
+ prevMetaIndex = metadata_[polygon->metadataId_].fileIndex_;
+
+ encodeVariableLength(outputBBox, polygon->fileIndex_ - prevFileIndex, false);
+ prevFileIndex = polygon->fileIndex_;
+ }
+ std::cout << "Encoded bounding box section into "<<outputBBox.size()<<" bytes.\n";
+
+ /* Encode header */
+ std::vector<uint8_t> outputHeader;
+ outputHeader.push_back('P');
+ outputHeader.push_back('L');
+ outputHeader.push_back('B');
+ outputHeader.push_back(tableType);
+ outputHeader.push_back(0);
+ outputHeader.push_back(precision);
+ outputHeader.push_back(fieldNames_.size());
+ for(unsigned int i=0; i<fieldNames_.size(); i++) {
+ encodeStringToBinary(outputHeader, fieldNames_[i]);
+ }
+ encodeStringToBinary(outputHeader, notice);
+ encodeVariableLength(outputHeader, outputBBox.size(), false);
+ encodeVariableLength(outputHeader, outputMeta.size(), false);
+ encodeVariableLength(outputHeader, outputData.size(), false);
+ std::cout << "Encoded header into "<<outputHeader.size()<<" bytes.\n";
+
+ FILE* outputFile = fopen(outPath.c_str(), "wb");
+ fwrite(outputHeader.data(), 1, outputHeader.size(), outputFile);
+ fwrite(outputBBox.data(), 1, outputBBox.size(), outputFile);
+ fwrite(outputMeta.data(), 1, outputMeta.size(), outputFile);
+ fwrite(outputData.data(), 1, outputData.size(), outputFile);
+ fclose(outputFile);
+
+}
projects / geeqie.git / blobdiff