Refactor shapefile spatial queries (systemed#228)

docs/CONFIGURATION.md

@@ -197,7 +197,6 @@ Or you can find out what country(/ies) the node is within using `FindIntersectin
 
 To enable these functions, set `index` to true in your shapefile layer definition. `index_column` is not needed for `Intersects` but required for `FindIntersecting`.
 
-Note these significant provisos:
+`CoveredBy` and `FindCovering` work similarly but check if the object is covered by a shapefile layer object.
 
-* Way queries are performed on the start and end points of ways, not the full way geometry. So if your way starts and ends outside a polygon, `Intersects` will return false, even if the midpoints are within the polygon. This may be changed in a future version.
-* Spatial queries do not work where the OSM object is a multipolygon, and will return false/empty. This may be changed in a future version.
+`AreaIntersecting` returns the area of the current way's intersection with the shapefile layer. You can use this to find whether a water body is already represented in a shapefile ocean layer.

include/osm_lua_processing.h

@@ -89,13 +89,21 @@ class OsmLuaProcessing : public PbfReaderOutput {
 
 	// Find intersecting shapefile layer
 	std::vector<std::string> FindIntersecting(const std::string &layerName);
+	double AreaIntersecting(const std::string &layerName);
 	bool Intersects(const std::string &layerName);
+	template <typename GeometryT> double intersectsArea(const std::string &layerName, GeometryT &geom) const;
+	template <typename GeometryT> std::vector<uint> intersectsQuery(const std::string &layerName, bool once, GeometryT &geom) const;
 
+	std::vector<std::string> FindCovering(const std::string &layerName);
+	bool CoveredBy(const std::string &layerName);
+	template <typename GeometryT> std::vector<uint> coveredQuery(const std::string &layerName, bool once, GeometryT &geom) const;
+
 	// Returns whether it is closed polygon
 	bool IsClosed() const;
 
 	// Returns area
 	double Area();
+	double multiPolygonArea(const MultiPolygon &mp) const;
 
 	// Returns length
 	double Length();
@@ -157,9 +165,8 @@ class OsmLuaProcessing : public PbfReaderOutput {
 		multiPolygonInited = false;
 	}
 
-	// Internal: set start/end co-ordinates
-	inline void setLocation(int32_t a, int32_t b, int32_t c, int32_t d) {
-		lon1=a; latp1=b; lon2=c; latp2=d;
+	const inline Point getPoint() {
+		return Point(lon/10000000.0,latp/10000000.0);
 	}
 
 	OSMStore const *indexStore;				// global OSM for reading input
@@ -176,7 +183,7 @@ class OsmLuaProcessing : public PbfReaderOutput {
 	WayID newWayID = MAX_WAY_ID;			///< Decrementing new ID for relations
 	bool isWay, isRelation, isClosed;		///< Way, node, relation?
 
-	int32_t lon1,latp1,lon2,latp2;			///< Start/end co-ordinates of OSM object
+	int32_t lon,latp;						///< Node coordinates
 	OSMStore::handle_t nodeVecHandle;
 	OSMStore::handle_t relationHandle;
 

include/output_object.h

@@ -125,9 +125,6 @@ Geometry buildWayGeometry(OSMStore &osmStore, OutputObject const &oo, const Tile
 //\brief Build a node geometry
 LatpLon buildNodeGeometry(OSMStore &osmStore, OutputObject const &oo, const TileBbox &bbox);
 
-///\brief Check if the object intersects with the given point
-bool intersects(OSMStore &osmStore, OutputObject const &oo, Point const &p);
-
 // Comparison functions
 
 bool operator==(const OutputObjectRef &x, const OutputObjectRef &y);

include/shp_mem_tiles.h

@@ -9,9 +9,6 @@ class ShpMemTiles : public TileDataSource
 public:
 	ShpMemTiles(OSMStore &osmStore, uint baseZoom);
 
-	// Find intersecting shapefile layer
-	std::vector<std::string> FindIntersecting(const std::string &layerName, Box &box) const;
-	bool Intersects(const std::string &layerName, Box &box) const;
 	void CreateNamedLayerIndex(const std::string &layerName);
 
 	// Used in shape file loading
@@ -24,11 +21,30 @@ class ShpMemTiles : public TileDataSource
 	void AddObject(TileCoordinates const &index, OutputObjectRef const &oo) {
 		tileIndex[index].push_back(oo);
 	}
-private:
-	std::vector<uint> findIntersectingGeometries(const std::string &layerName, Box &box) const;
-	std::vector<uint> verifyIntersectResults(std::vector<IndexValue> &results, Point &p1, Point &p2) const;
-	std::vector<std::string> namesOfGeometries(std::vector<uint> &ids) const;
+	std::vector<uint> QueryMatchingGeometries(const std::string &layerName, bool once, Box &box, 
+		std::function<std::vector<IndexValue>(const RTree &rtree)> indexQuery, 
+		std::function<bool(OutputObject &oo)> checkQuery) const;
+	std::vector<std::string> namesOfGeometries(const std::vector<uint> &ids) const;
+
+	template <typename GeometryT>
+	double AreaIntersecting(const std::string &layerName, GeometryT &g) const {
+		auto f = indices.find(layerName);
+		if (f==indices.end()) { std::cerr << "Couldn't find indexed layer " << layerName << std::endl; return false;  }
+		Box box; geom::envelope(g, box);
+		std::vector <IndexValue> results;
+		f->second.query(geom::index::intersects(box), back_inserter(results));
+		MultiPolygon mp, tmp;
+		for (auto it : results) {
+			OutputObjectRef oo = cachedGeometries.at(it.second);
+			if (oo->geomType!=OutputGeometryType::POLYGON) continue;
+			geom::union_(mp, osmStore.retrieve<mmap::multi_polygon_t>(oo->handle), tmp);
+			geom::assign(mp, tmp);
+		}
+		geom::correct(mp);
+		return geom::covered_by(g, mp);
+	}
 
+private:
 	/// Add an OutputObject to all tiles between min/max lat/lon
 	void addToTileIndexByBbox(OutputObjectRef &oo, 
 		double minLon, double minLatp, double maxLon, double maxLatp);

resources/config-openmaptiles.json

@@ -22,7 +22,7 @@
 		"building":          { "minzoom": 13, "maxzoom": 14 },
 
 		"water":             { "minzoom": 6,  "maxzoom": 14, "simplify_below": 12, "simplify_length": 1, "simplify_ratio": 2},
-		"ocean":             { "minzoom": 0,  "maxzoom": 14, "source": "coastline/water_polygons.shp", "simplify_below": 13, "simplify_level": 0.0005, "write_to": "water" },
+		"ocean":             { "minzoom": 0,  "maxzoom": 14, "source": "coastline/water_polygons.shp", "simplify_below": 13, "simplify_level": 0.0005, "write_to": "water", "index": true },
 		"water_name":        { "minzoom": 14, "maxzoom": 14 },
 		"water_name_detail": { "minzoom": 14, "maxzoom": 14, "write_to": "water_name" },
 

resources/process-openmaptiles.lua

@@ -396,6 +396,7 @@ function way_function(way)
 		if way:Find("covered")=="yes" or not isClosed then return end
 		local class="lake"; if natural=="bay" then class="ocean" elseif waterway~="" then class="river" end
 		if class=="lake" and way:Find("wikidata")=="Q192770" then return end
+		if class=="ocean" and isClosed and (way:AreaIntersecting("ocean")/way:Area() > 0.98) then return end
 		way:Layer("water",true)
 		SetMinZoomByArea(way)
 		way:Attribute("class",class)

src/osm_lua_processing.cpp

@@ -42,8 +42,11 @@ OsmLuaProcessing::OsmLuaProcessing(
 		.addFunction("Find", &OsmLuaProcessing::Find)
 		.addFunction("FindIntersecting", &OsmLuaProcessing::FindIntersecting)
 		.addFunction("Intersects", &OsmLuaProcessing::Intersects)
+		.addFunction("FindCovering", &OsmLuaProcessing::FindCovering)
+		.addFunction("CoveredBy", &OsmLuaProcessing::CoveredBy)
 		.addFunction("IsClosed", &OsmLuaProcessing::IsClosed)
 		.addFunction("Area", &OsmLuaProcessing::Area)
+		.addFunction("AreaIntersecting", &OsmLuaProcessing::AreaIntersecting)
 		.addFunction("Length", &OsmLuaProcessing::Length)
 		.addFunction("Layer", &OsmLuaProcessing::Layer)
 		.addFunction("LayerAsCentroid", &OsmLuaProcessing::LayerAsCentroid)
@@ -110,23 +113,92 @@ string OsmLuaProcessing::Find(const string& key) const {
 
 // ----	Spatial queries called from Lua
 
-// Find intersecting shapefile layer
 vector<string> OsmLuaProcessing::FindIntersecting(const string &layerName) {
-	// TODO: multipolygon relations not supported, will always return empty vector
-	if(isRelation) return vector<string>();
-	Point p1(lon1/10000000.0,latp1/10000000.0);
-	Point p2(lon2/10000000.0,latp2/10000000.0);
-	Box box = Box(p1,p2);
-	return shpMemTiles.FindIntersecting(layerName, box);
+	if      (!isWay   ) { return shpMemTiles.namesOfGeometries(intersectsQuery(layerName, false, getPoint())); }
+	else if (!isClosed) { return shpMemTiles.namesOfGeometries(intersectsQuery(layerName, false, linestringCached())); }
+	else if (isRelation){ return shpMemTiles.namesOfGeometries(intersectsQuery(layerName, false, multiPolygonCached())); }
+	else                { return shpMemTiles.namesOfGeometries(intersectsQuery(layerName, false, polygonCached())); }
 }
 
 bool OsmLuaProcessing::Intersects(const string &layerName) {
-	// TODO: multipolygon relations not supported, will always return false
-	if(isRelation) return false;
-	Point p1(lon1/10000000.0,latp1/10000000.0);
-	Point p2(lon2/10000000.0,latp2/10000000.0);
-	Box box = Box(p1,p2);
-	return shpMemTiles.Intersects(layerName, box);
+	if      (!isWay   ) { return !intersectsQuery(layerName, true, getPoint()).empty(); }
+	else if (!isClosed) { return !intersectsQuery(layerName, true, linestringCached()).empty(); }
+	else if (isRelation){ return !intersectsQuery(layerName, true, multiPolygonCached()).empty(); }
+	else                { return !intersectsQuery(layerName, true, polygonCached()).empty(); }
+}
+
+vector<string> OsmLuaProcessing::FindCovering(const string &layerName) {
+	if      (!isWay   ) { return shpMemTiles.namesOfGeometries(coveredQuery(layerName, false, getPoint())); }
+	else if (!isClosed) { return shpMemTiles.namesOfGeometries(coveredQuery(layerName, false, linestringCached())); }
+	else if (isRelation){ return shpMemTiles.namesOfGeometries(coveredQuery(layerName, false, multiPolygonCached())); }
+	else                { return shpMemTiles.namesOfGeometries(coveredQuery(layerName, false, polygonCached())); }
+}
+
+bool OsmLuaProcessing::CoveredBy(const string &layerName) {
+	if      (!isWay   ) { return !coveredQuery(layerName, true, getPoint()).empty(); }
+	else if (!isClosed) { return !coveredQuery(layerName, true, linestringCached()).empty(); }
+	else if (isRelation){ return !coveredQuery(layerName, true, multiPolygonCached()).empty(); }
+	else                { return !coveredQuery(layerName, true, polygonCached()).empty(); }
+}
+
+double OsmLuaProcessing::AreaIntersecting(const string &layerName) {
+	if      (!isWay || !isClosed) { return 0.0; }
+	else if (isRelation){ return intersectsArea(layerName, multiPolygonCached()); }
+	else                { return intersectsArea(layerName, polygonCached()); }
+}
+
+
+template <typename GeometryT>
+std::vector<uint> OsmLuaProcessing::intersectsQuery(const string &layerName, bool once, GeometryT &geom) const {
+	Box box; geom::envelope(geom, box);
+	std::vector<uint> ids = shpMemTiles.QueryMatchingGeometries(layerName, once, box,
+		[&](const RTree &rtree) { // indexQuery
+			vector<IndexValue> results;
+			rtree.query(geom::index::intersects(box), back_inserter(results));
+			return results;
+		},
+		[&](OutputObject &oo) { // checkQuery
+			return geom::intersects(geom, osmStore.retrieve<mmap::multi_polygon_t>(oo.handle));
+		}
+	);
+	return ids;
+}
+
+template <typename GeometryT>
+double OsmLuaProcessing::intersectsArea(const string &layerName, GeometryT &geom) const {
+	Box box; geom::envelope(geom, box);
+	double area = 0.0;
+	std::vector<uint> ids = shpMemTiles.QueryMatchingGeometries(layerName, false, box,
+		[&](const RTree &rtree) { // indexQuery
+			vector<IndexValue> results;
+			rtree.query(geom::index::intersects(box), back_inserter(results));
+			return results;
+		},
+		[&](OutputObject &oo) { // checkQuery
+			MultiPolygon tmp;
+			geom::intersection(geom, osmStore.retrieve<mmap::multi_polygon_t>(oo.handle), tmp);
+			area += multiPolygonArea(tmp);
+			return false;
+		}
+	);
+	return area;
+}
+
+template <typename GeometryT>
+std::vector<uint> OsmLuaProcessing::coveredQuery(const string &layerName, bool once, GeometryT &geom) const {
+	Box box; geom::envelope(geom, box);
+	std::vector<uint> ids = shpMemTiles.QueryMatchingGeometries(layerName, once, box,
+		[&](const RTree &rtree) { // indexQuery
+			vector<IndexValue> results;
+			rtree.query(geom::index::intersects(box), back_inserter(results));
+			return results;
+		},
+		[&](OutputObject &oo) { // checkQuery
+			if (oo.geomType!=OutputGeometryType::POLYGON) return false; // can only be covered by a polygon!
+			return geom::covered_by(geom, osmStore.retrieve<mmap::multi_polygon_t>(oo.handle));
+		}
+	);
+	return ids;
 }
 
 // Returns whether it is closed polygon
@@ -148,16 +220,7 @@ double OsmLuaProcessing::Area() {
 	geom::strategy::area::spherical<> sph_strategy(RadiusMeter);
 	if (isRelation) {
 		// Boost won't calculate area of a multipolygon, so we just total up the member polygons
-		double totalArea = 0;
-		MultiPolygon mp = multiPolygonCached();
-		for (MultiPolygon::const_iterator it = mp.begin(); it != mp.end(); ++it) {
-			geom::model::polygon<DegPoint> p;
-			geom::assign(p,*it);
-			geom::for_each_point(p, reverse_project);
-			double area = geom::area(p, sph_strategy);
-			totalArea += area;
-		}
-		return totalArea;
+		return multiPolygonArea(multiPolygonCached());
 	} else if (isWay) {
 		// Reproject back into lat/lon and then run Boo
 		geom::model::polygon<DegPoint> p;
@@ -176,6 +239,18 @@ double OsmLuaProcessing::Area() {
 	return 0;
 }
 
+double OsmLuaProcessing::multiPolygonArea(const MultiPolygon &mp) const {
+	geom::strategy::area::spherical<> sph_strategy(RadiusMeter);
+	double totalArea = 0;
+	for (MultiPolygon::const_iterator it = mp.begin(); it != mp.end(); ++it) {
+		geom::model::polygon<DegPoint> p;
+		geom::assign(p,*it);
+		geom::for_each_point(p, reverse_project);
+		totalArea += geom::area(p, sph_strategy);
+	}
+	return totalArea;
+}
+
 // Returns length
 double OsmLuaProcessing::Length() {
 	if (isWay) {
@@ -384,9 +459,8 @@ void OsmLuaProcessing::setNode(NodeID id, LatpLon node, const tag_map_t &tags) {
 	originalOsmID = id;
 	isWay = false;
 	isRelation = false;
-
-	setLocation(node.lon, node.latp, node.lon, node.latp);
-
+	lon = node.lon;
+	latp= node.latp;
 	currentTags = tags;
 
 	//Start Lua processing for node
@@ -418,8 +492,6 @@ void OsmLuaProcessing::setWay(WayID wayId, OSMStore::handle_t handle, const tag_
 	try {
 		auto const &nodeVecPtr = &indexStore->retrieve<WayStore::nodeid_vector_t>(nodeVecHandle);
 		isClosed = nodeVecPtr->front()==nodeVecPtr->back();
-		setLocation(indexStore->nodes_at(nodeVecPtr->front()).lon, indexStore->nodes_at(nodeVecPtr->front()).latp,
-				indexStore->nodes_at(nodeVecPtr->back()).lon, indexStore->nodes_at(nodeVecPtr->back()).latp);
 
 	} catch (std::out_of_range &err) {
 		std::stringstream ss;
@@ -500,9 +572,9 @@ void OsmLuaProcessing::setRelation(int64_t relationId, OSMStore::handle_t relati
 	originalOsmID = relationId;
 	isWay = true;
 	isRelation = true;
+	isClosed = true; // TODO: change this when we support route relations
 
 	this->relationHandle = relationHandle;
-	//setLocation(...); TODO
 
 	currentTags = tags;
 	/* currentTags.clear();

src/output_object.cpp

@@ -133,24 +133,6 @@ LatpLon buildNodeGeometry(OSMStore &osmStore, OutputObject const &oo, const Tile
 	throw std::runtime_error("Geometry type is not point");			
 }
 
-bool intersects(OSMStore &osmStore, OutputObject const &oo, Point const &p)
-{
-	switch(oo.geomType) {
-		case OutputGeometryType::POINT:
-			return boost::geometry::intersects(osmStore.retrieve<mmap::point_t>(oo.handle), p);
-
-		case OutputGeometryType::LINESTRING:
-			return boost::geometry::intersects(osmStore.retrieve<mmap::linestring_t>(oo.handle), p);
-
-
-		case OutputGeometryType::POLYGON:
-			return boost::geometry::intersects(osmStore.retrieve<mmap::multi_polygon_t>(oo.handle), p);
-
-		default:
-			throw std::runtime_error("Invalid output geometry");
-	}
-}
-
 // Find a value in the value dictionary
 // (we can't easily use find() because of the different value-type encoding - 
 //	should be possible to improve this though)