cuebiq-subset.cpp

// Compile with: g++ ./cuebiq-subset.cpp -std=c++11 -Wall -O3 -funroll-loops -o cuebiq-subset
#include <algorithm>
#include <fstream>
#include <iostream>
#include <string>

#include <boost/geometry.hpp>
#include <boost/geometry/geometries/point_xy.hpp>
#include <boost/geometry/geometries/polygon.hpp>

typedef boost::geometry::model::d2::point_xy<double> point_type;
const int LATITUDE_COLUMN = 4;
const char COLUMN_SEPARATOR = '\t';

int main(int argc, char* argv[]) {

    if (argc != 2) {
        std::cout << "usage: " << argv[0] << " path/to/bounds.wkt" << std::endl;
        return 1;
    }

    // Used for extracting latitudes and longitudes from cuebiq data
    size_t latitude_start = 0;
    size_t longitude_start = 0;
    size_t longitude_end = 0;
    std::string latitude_str;
    std::string longitude_str;

    // Used for reading WKT from file
    std::ifstream bounds_wkt_file;
    std::stringstream bounds_buffer;
    std::string bounds_wkt_str;

    // Geometry types
    boost::geometry::model::polygon<point_type> bounds;
    point_type point;

    // Read WKT from file into a boost::geometry::model::polygon object
    bounds_wkt_file.open(argv[1]);
    bounds_buffer << bounds_wkt_file.rdbuf();
    bounds_wkt_str = bounds_buffer.str();
    bounds_wkt_str.erase(
        std::remove(bounds_wkt_str.begin(), bounds_wkt_str.end(), '\n'),
        bounds_wkt_str.end()
    );
    boost::geometry::read_wkt(bounds_wkt_str, bounds);

    for (std::string line; std::getline(std::cin, line);) {

        // Find tab breaks for latitude and longitude columns
        latitude_start = 0;
        for (int i = 0; i < LATITUDE_COLUMN - 1; i++) {
            latitude_start = line.find(COLUMN_SEPARATOR, latitude_start + 1);
        }
        longitude_start = line.find(COLUMN_SEPARATOR, latitude_start + 1);
        longitude_end = line.find(COLUMN_SEPARATOR, longitude_start + 1);

        // Small optimization: only process the longitude if the latitude is valid
        latitude_str = line.substr(latitude_start + 1, longitude_start - latitude_start - 1);
        if (latitude_str.length() > 0) {
            longitude_str = line.substr(longitude_start + 1, longitude_end - longitude_start - 1);

            // Overwrite the point object and do a point-in-polygon
            point.set<0>(std::stod(longitude_str));
            point.set<1>(std::stod(latitude_str));
            if (boost::geometry::within(point, bounds)) {
                std::cout << line << std::endl;
            }
        }

    }

    return 0;
}