Using polygons with libspatialindex
An answer to this question on the GIS Stack Exchange.
Question
I am investigating using some C/C++ libraries for point-in-polygon lookups and libspatialindex seems popular. I've written some code that seems to do what I want but I'm not entirely sure if what I doing is correct. I have a requirement that the polygons have an arbitrary number of dimensions, they generally range from 5 to sometimes 16.
My problem is that when creating the RTree I have to specify the number of dimensions and so if I try to insert a polygon that has a higher number the insertion fails. Also I'm not quite sure if SpatialIndex::Region is the correct object to use to represent polygons, it seems to actually be a rectangle rather that what I'd call a Region as it only allows 2 dimensions.
A quick fix would be to create individual Regions for each pair of points in the polygon but that is quite inaccurate.
I've scoured the internet for sample code but I can't find anything worthwhile.
Am I doing something wrong?
Answer
This is my solution for using RTree to test polygons.
The idea is to load all of the polygons from a shape file. For simplicity, I'll only look at their exterior rings.
I then find the bounding boxes of each of the polygons and load the bounding boxes into the RTree along with a numerical id corresponding to that polygon's position in the list of loaded polygons.
Now, to see if a point is within a polygon I use the RTree to quickly determine which bounding boxes, if any, it falls into. If it falls into no box, then the point is not in a polygon. If it falls into one or more boxes, I then perform a much more expensive check to see if it actually lies within one of the polygons.
In true hacker style, I perform the more expensive check using Random Code From The Internetz. But, don't despair: I've added several links to folks explaining how the point-in-polygon check works.
I return the number of polygons the point falls into, but, obviously, you could do other things as well. Arbitrary things. Things without names.
Here's the code:
//Compile with: g++ `gdal-config --cflags` -march=native -mtune=native --std=c++11 -O3 -g -Wall -ffast-math -o test.exe test.cpp `gdal-config --libs` -lproj -lspatialindex_c -lspatialindex -lproj
#include <gdal/ogrsf_frmts.h>
#include <proj_api.h>
#include <spatialindex/capi/sidx_api.h>
#include <spatialindex/capi/sidx_impl.h>
#include <spatialindex/capi/sidx_config.h>
#include <iostream>
#include <vector>
#include <algorithm>
#include <array>
#include <iomanip>
#include <cassert>
class LineString {
public:
std::vector<double> x;
std::vector<double> y;
};
//Info: https://wrf.ecse.rpi.edu//Research/Short_Notes/pnpoly.html
//Info: https://stackoverflow.com/questions/8721406/how-to-determine-if-a-point-is-inside-a-2d-convex-polygon/23223947#23223947
//Info: http://stackoverflow.com/a/2922778/752843
int ContainsPoint(const std::vector<double> &vertx, const std::vector<double> &verty, const double testx, const double testy) {
unsigned int i, j;
int c = 0;
for (i = 0, j = vertx.size()-1; i < vertx.size(); j = i++) {
if ( ((verty[i]>testy) != (verty[j]>testy)) &&
(testx < (vertx[j]-vertx[i]) * (testy-verty[i]) / (verty[j]-verty[i]) + vertx[i]) )
c = !c;
}
return c;
}
class Polygon {
public:
LineString interior;
LineString exterior;
double minX() const {
return *std::min_element(exterior.x.begin(), exterior.x.end());
}
double maxX() const {
return *std::max_element(exterior.x.begin(), exterior.x.end());
}
double minY() const {
return *std::min_element(exterior.y.begin(), exterior.y.end());
}
double maxY() const {
return *std::max_element(exterior.y.begin(), exterior.y.end());
}
bool containsPoint(const double testx, const double testy) const {
return ContainsPoint(exterior.x, exterior.y, testx, testy);
}
};
void ReadShapefile(std::string filename, std::string layername, std::vector<Polygon> &geometries){
GDALAllRegister();
GDALDataset *poDS;
poDS = (GDALDataset*) GDALOpenEx(filename.c_str(), GDAL_OF_VECTOR, NULL, NULL, NULL);
if( poDS == NULL ){
std::cerr<<"Failed to open '"<<filename<<"'!"<<std::endl;
exit( 1 );
}
OGRLayer *poLayer;
poLayer = poDS->GetLayerByName(layername.c_str());
OGRFeature *poFeature;
poLayer->ResetReading();
while( (poFeature = poLayer->GetNextFeature()) != NULL ){
/*OGRFeatureDefn *poFDefn = poLayer->GetLayerDefn();
int iField;
for( iField = 0; iField < poFDefn->GetFieldCount(); iField++ ){
OGRFieldDefn *poFieldDefn = poFDefn->GetFieldDefn( iField );
if( poFieldDefn->GetType() == OFTInteger )
printf( "%d,", poFeature->GetFieldAsInteger( iField ) );
else if( poFieldDefn->GetType() == OFTInteger64 )
printf( CPL_FRMT_GIB ",", poFeature->GetFieldAsInteger64( iField ) );
else if( poFieldDefn->GetType() == OFTReal )
printf( "%.3f,", poFeature->GetFieldAsDouble(iField) );
else if( poFieldDefn->GetType() == OFTString )
printf( "%s,", poFeature->GetFieldAsString(iField) );
else
printf( "%s,", poFeature->GetFieldAsString(iField) );
}*/
OGRGeometry *poGeometry;
poGeometry = poFeature->GetGeometryRef();
if (poGeometry==NULL){
//Pass
} else if( wkbFlatten(poGeometry->getGeometryType()) == wkbPolygon ){
OGRPolygon *poly = (OGRPolygon *) poGeometry;
auto extring = poly->getExteriorRing();
//Ignore interior rings for now: they're probably lakes
geometries.emplace_back();
for(int i=0;i<extring->getNumPoints();i++){
geometries.back().exterior.x.emplace_back(extring->getX(i));
geometries.back().exterior.y.emplace_back(extring->getY(i));
}
} else {
std::cerr<<"Unrecognised geometry of type: "<<wkbFlatten(poGeometry->getGeometryType())<<std::endl;
}
OGRFeature::DestroyFeature( poFeature );
}
GDALClose( poDS );
}
//Saving to disk: https://stackoverflow.com/questions/33395445/libspatialindex-loading-storing-index-on-disk
class SpIndex {
private:
Index* idx;
public:
SpIndex(){
//Create a property set with default values.
//See utility.cc for all defaults http://libspatialindex.github.io/doxygen/Utility_8cc_source.html#l00031
Tools::PropertySet* ps = GetDefaults();
Tools::Variant var;
//Set index type to R*-Tree
var.m_varType = Tools::VT_ULONG;
var.m_val.ulVal = RT_RTree;
ps->setProperty("IndexType", var);
//Set index to store in memory (default is disk)
var.m_varType = Tools::VT_ULONG;
var.m_val.ulVal = RT_Memory;
ps->setProperty("IndexStorageType", var);
//Initalise index
idx = new Index(*ps);
delete ps;
//Check index is ok
if (!idx->index().isIndexValid())
throw std::runtime_error("Failed to create valid index");
}
void addBox(double left, double bottom, double right, double top, int64_t id){
//Create array with lat/lon points
double low[] = {left,bottom};
double high[] = {right,top};
//Shapes can also have an object associated with them, but not this one!
uint8_t* pData = 0;
uint32_t nDataLength = 0;
// create shape
auto shape = SpatialIndex::Region(low,high,2);
// insert into index along with the an object and an ID
idx->index().insertData(nDataLength,pData,shape,id);
}
void addPolygon(const Polygon &p, int64_t id){
double left = p.minX();
double bottom = p.minY();
double right = p.maxX();
double top = p.maxY();
addBox(left,bottom,right,top, id);
}
std::vector<int64_t> overlaps(double x, double y) const {
double coords[] = {x,y};
//Object that will show us what results we've found
ObjVisitor visitor;
//Two-dimensional query point
SpatialIndex::Point r(coords, 2);
//Find those regions that intersect with the query point
idx->index().intersectsWithQuery(r,visitor);
//Copy results
std::vector<int64_t> temp;
for(const auto &r: visitor.GetResults())
temp.push_back(r->getIdentifier()); // dynamic_cast<SpatialIndex::IData*>(results[i]->clone()));
//std::cout << "found " << visitor.GetResultCount() << " results." << std::endl;
return temp;
}
};
int CountOverlaps(const int x, const int y, const SpIndex &sp, const std::vector<Polygon> &polygons){
auto sp_overlaps = sp.overlaps(x,y);
//Point doesn't intersect any bounding boxes, so it can't overlap a polygon
if(sp_overlaps.size()==0)
return 0;
//Within a bounding box, but is it actually within the polygon?
int overlaps = 0;
for(auto &pi: sp_overlaps){
if(polygons[pi].containsPoint(x,y)){
overlaps++;
break;
}
}
return overlaps;
}
int main(int argc, char **argv){
std::vector<Polygon> geometries;
ReadShapefile("data/simplified-land-polygons-complete-3857/simplified_land_polygons.shp", "simplified_land_polygons", geometries);
//ReadShapefile("data/land-polygons-split-3857/land_polygons.shp", "land_polygons", geometries);
std::cerr<<"Read "<<geometries.size()<<" polygons."<<std::endl;
SpIndex sp;
std::cerr<<"Building index..."<<std::endl;
for(int64_t i=0;(unsigned)i<geometries.size();i++) //TODO: Fix this unsigned int64 mess
sp.addPolygon(geometries[i], i);
std::cerr<<"Index built."<<std::endl;
std::cerr<<"Number of overlaps = "<<CountOverlaps(-11501094,3771603,sp,geometries)<<std::endl;
}