Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Make shapefile processing multi-threaded #614

Merged
merged 3 commits into from
Dec 15, 2023
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
11 changes: 8 additions & 3 deletions include/read_shp.h
Original file line number Diff line number Diff line change
Expand Up @@ -21,14 +21,19 @@ AttributeIndex readShapefileAttributes(DBFHandle &dbf, int recordNum,
std::unordered_map<int,std::string> &columnMap,
std::unordered_map<int,int> &columnTypeMap,
LayerDef &layer,
OsmLuaProcessing &osmLuaProcessing, int &minzoom);
OsmLuaProcessing &osmLuaProcessing, uint &minzoom);

/// Read shapefile, and create OutputObjects for all objects within the specified bounding box
void readShapefile(const Box &clippingBox,
class LayerDefinition &layers,
uint baseZoom, uint layerNum,
class ShpMemTiles &shpMemTiles,
OsmLuaProcessing &osmLuaProcessing);
uint threadNum,
class ShpMemTiles &shpMemTiles,
OsmLuaProcessing &osmLuaProcessing);

// Process an individual shapefile record
void processShapeGeometry(SHPObject* shape, AttributeIndex attrIdx, ShpMemTiles &shpMemTiles,
const Box &clippingBox, const LayerDef &layer, uint layerNum, bool hasName, const std::string &name);

#endif //_READ_SHP_H

215 changes: 106 additions & 109 deletions src/read_shp.cpp
Original file line number Diff line number Diff line change
@@ -1,12 +1,19 @@
#include "read_shp.h"

#include <boost/asio/thread_pool.hpp>
#include <boost/asio/post.hpp>

extern bool verbose;

using namespace std;
namespace geom = boost::geometry;

/*
Read shapefiles into Boost geometries
*/

std::mutex attributeMutex;

void fillPointArrayFromShapefile(vector<Point> *points, SHPObject *shape, uint part) {
uint start = shape->panPartStart[part];
uint end = (int(part)==shape->nParts-1) ? shape->nVertices : shape->panPartStart[part+1];
Expand Down Expand Up @@ -34,8 +41,9 @@ AttributeIndex readShapefileAttributes(
DBFHandle &dbf,
int recordNum, unordered_map<int,string> &columnMap, unordered_map<int,int> &columnTypeMap,
LayerDef &layer,
OsmLuaProcessing &osmLuaProcessing, int &minzoom) {
OsmLuaProcessing &osmLuaProcessing, uint &minzoom) {

std::lock_guard<std::mutex> lock(attributeMutex);
AttributeStore& attributeStore = osmLuaProcessing.getAttributeStore();

AttributeSet attributes;
Expand Down Expand Up @@ -99,16 +107,15 @@ AttributeIndex readShapefileAttributes(

// Read shapefile, and create OutputObjects for all objects within the specified bounding box
void readShapefile(const Box &clippingBox,
class LayerDefinition &layers,
class LayerDefinition &layers,
uint baseZoom, uint layerNum,
class ShpMemTiles &shpMemTiles,
OsmLuaProcessing &osmLuaProcessing)
uint threadNum,
class ShpMemTiles &shpMemTiles,
OsmLuaProcessing &osmLuaProcessing)
{
LayerDef &layer = layers.layers[layerNum];
const string &filename = layer.source;
const vector<string> &columns = layer.sourceColumns;
const string &layerName = layer.name;
bool isIndexed = layer.indexed;
const string &indexName = layer.indexName;

// open shapefile
Expand All @@ -117,7 +124,6 @@ void readShapefile(const Box &clippingBox,
if(shp == nullptr || dbf == nullptr)
return;
int numEntities=0, shpType=0;
vector<Point> points;
double adfMinBound[4], adfMaxBound[4];
SHPGetInfo(shp, &numEntities, &shpType, adfMinBound, adfMaxBound);

Expand All @@ -142,13 +148,11 @@ void readShapefile(const Box &clippingBox,
int indexField=-1;
if (indexName!="") { indexField = DBFGetFieldIndex(dbf,indexName.c_str()); }

boost::asio::thread_pool pool(threadNum);
for (int i=0; i<numEntities; i++) {
SHPObject* shape = SHPReadObject(shp, i);
if(shape == nullptr) { cerr << "Error loading shape from shapefile" << endl; continue; }

int shapeType = shape->nSHPType; // 1=point, 3=polyline, 5=(multi)polygon [8=multipoint, 11+=3D]
int minzoom = layer.minzoom;

// Check shape is in clippingBox
Box shapeBox(Point(shape->dfXMin, lat2latp(shape->dfYMin)), Point(shape->dfXMax, lat2latp(shape->dfYMax)));
if (shapeBox.min_corner().get<0>() > clippingBox.max_corner().get<0>() ||
Expand All @@ -159,117 +163,110 @@ void readShapefile(const Box &clippingBox,
continue;
}

if (shapeType==1) {
// Points
Point p( shape->padfX[0], lat2latp(shape->padfY[0]) );
if (geom::within(p, clippingBox)) {

string name;
bool hasName = false;
if (indexField>-1) { name=DBFReadStringAttribute(dbf, i, indexField); hasName = true;}

AttributeIndex attrIdx = readShapefileAttributes(dbf, i, columnMap, columnTypeMap, layer, osmLuaProcessing, minzoom);
shpMemTiles.StoreShapefileGeometry(layerNum, layerName, POINT_, p, isIndexed, hasName, name, minzoom, attrIdx);
}
boost::asio::post(pool, [&, shape]() {
// process attributes
string name;
bool hasName = false;
if (indexField>-1) { name=DBFReadStringAttribute(dbf, i, indexField); hasName = true;}
AttributeIndex attrIdx = readShapefileAttributes(dbf, i, columnMap, columnTypeMap, layer, osmLuaProcessing, layer.minzoom);
// process geometry
processShapeGeometry(shape, attrIdx, shpMemTiles, clippingBox, layer, layerNum, hasName, name);
SHPDestroyObject(shape);
});
}
pool.join();
SHPClose(shp);
DBFClose(dbf);
}

} else if (shapeType==3) {
// (Multi)-polylines
// Due to https://svn.boost.org/trac/boost/ticket/11268, we can't clip a MultiLinestring with Boost 1.56-1.58,
// so we need to create everything as polylines and clip individually :(
for (int j=0; j<shape->nParts; j++) {
Linestring ls;
fillPointArrayFromShapefile(&points, shape, j);
geom::assign_points(ls, points);
MultiLinestring out;
geom::intersection(ls, clippingBox, out);
for (MultiLinestring::const_iterator it = out.begin(); it != out.end(); ++it) {
void processShapeGeometry(SHPObject* shape, AttributeIndex attrIdx, ShpMemTiles &shpMemTiles,
const Box &clippingBox, const LayerDef &layer, uint layerNum, bool hasName, const string &name) {
int shapeType = shape->nSHPType; // 1=point, 3=polyline, 5=(multi)polygon [8=multipoint, 11+=3D]
int minzoom = layer.minzoom;

string name;
bool hasName = false;
if (indexField>-1) { name=DBFReadStringAttribute(dbf, i, indexField); hasName = true;}
if (shapeType==1) {
// Points
Point p( shape->padfX[0], lat2latp(shape->padfY[0]) );
if (geom::within(p, clippingBox)) {
shpMemTiles.StoreShapefileGeometry(layerNum, layer.name, POINT_, p, layer.indexed, hasName, name, minzoom, attrIdx);
}

AttributeIndex attrIdx = readShapefileAttributes(dbf, i, columnMap, columnTypeMap, layer, osmLuaProcessing, minzoom);
shpMemTiles.StoreShapefileGeometry(layerNum, layerName, LINESTRING_, *it, isIndexed, hasName, name, minzoom, attrIdx);
}
} else if (shapeType==3) {
// (Multi)-polylines
// Due to https://svn.boost.org/trac/boost/ticket/11268, we can't clip a MultiLinestring with Boost 1.56-1.58,
// so we need to create everything as polylines and clip individually :(
vector<Point> points;
for (int j=0; j<shape->nParts; j++) {
Linestring ls;
fillPointArrayFromShapefile(&points, shape, j);
geom::assign_points(ls, points);
MultiLinestring out;
geom::intersection(ls, clippingBox, out);
for (MultiLinestring::const_iterator it = out.begin(); it != out.end(); ++it) {
shpMemTiles.StoreShapefileGeometry(layerNum, layer.name, LINESTRING_, *it, layer.indexed, hasName, name, minzoom, attrIdx);
}
}

} else if (shapeType==5) {
// (Multi)-polygons
MultiPolygon multi;
Polygon poly;
Ring ring;
int nInteriorRings = 0;
} else if (shapeType==5) {
// (Multi)-polygons
MultiPolygon multi;
Polygon poly;
Ring ring;
int nInteriorRings = 0;
vector<Point> points;

// To avoid expensive computations, we assume the shapefile has been pre-processed
// such that each polygon's exterior ring is immediately followed by its interior rings.
for (int j=0; j<shape->nParts; j++) {
fillPointArrayFromShapefile(&points, shape, j);
// Read points into a ring
ring.clear();
geom::append(ring, points);
// To avoid expensive computations, we assume the shapefile has been pre-processed
// such that each polygon's exterior ring is immediately followed by its interior rings.
for (int j=0; j<shape->nParts; j++) {
fillPointArrayFromShapefile(&points, shape, j);
// Read points into a ring
ring.clear();
geom::append(ring, points);

if (j == 0) {
// We assume the first part is an exterior ring of the first polygon.
geom::append(poly, ring);
}
else if (geom::area(ring) > 0.0) {
// This part has clockwise orientation - an exterior ring.
// Start a new polygon.
multi.push_back(poly);
poly.clear();
nInteriorRings = 0;
geom::append(poly, ring);
} else {
// This part has anti-clockwise orientation.
// Add another interior ring to the current polygon.
nInteriorRings++;
geom::interior_rings(poly).resize(nInteriorRings);
geom::append(poly, ring, nInteriorRings - 1);
}
if (j == 0) {
// We assume the first part is an exterior ring of the first polygon.
geom::append(poly, ring);
}
// All parts read. Add the last polygon.
multi.push_back(poly);
geom::remove_spikes(multi);

string reason;
#if BOOST_VERSION >= 105800
if (!geom::is_valid(multi, reason)) {
cerr << "Shapefile entity #" << i << " type " << shapeType << " is invalid. Parts:" << shape->nParts << ". Reason:" << reason;

// Perform make_valid operation
make_valid(multi);

if (geom::is_valid(multi, reason)) {
cerr << "... corrected";
} else {
cerr << "... failed to correct. Reason: " << reason;
}
cerr << endl;
else if (geom::area(ring) > 0.0) {
// This part has clockwise orientation - an exterior ring.
// Start a new polygon.
multi.push_back(poly);
poly.clear();
nInteriorRings = 0;
geom::append(poly, ring);
} else {
// This part has anti-clockwise orientation.
// Add another interior ring to the current polygon.
nInteriorRings++;
geom::interior_rings(poly).resize(nInteriorRings);
geom::append(poly, ring, nInteriorRings - 1);
}
#else
if (!geom::is_valid(multi)) { geom::correct(multi); geom::remove_spikes(multi); }
#endif
// clip to bounding box
MultiPolygon out;
geom::intersection(multi, clippingBox, out);
if (boost::size(out)>0) {
}

string name;
bool hasName = false;
if (indexField>-1) { name=DBFReadStringAttribute(dbf, i, indexField); hasName = true;}
// All parts read. Add the last polygon.
multi.push_back(poly);
geom::remove_spikes(multi);

// create OutputObject
AttributeIndex attrIdx = readShapefileAttributes(dbf, i, columnMap, columnTypeMap, layer, osmLuaProcessing, minzoom);
shpMemTiles.StoreShapefileGeometry(layerNum, layerName, POLYGON_, out, isIndexed, hasName, name, minzoom, attrIdx);
// Make valid if needs be
string reason;
if (!geom::is_valid(multi, reason)) {
if (verbose) cerr << "Shapefile entity " << shape->nShapeId << " type " << shapeType << " is invalid. Parts:" << shape->nParts << ". Reason:" << reason;
make_valid(multi);
if (verbose) {
if (geom::is_valid(multi, reason)) { cerr << "... corrected"; }
else { cerr << "... failed to correct. Reason: " << reason; }
cerr << endl;
}

} else {
// Not supported
cerr << "Shapefile entity #" << i << " type " << shapeType << " not supported" << endl;
}
SHPDestroyObject(shape);
}
// clip to bounding box
MultiPolygon out;
geom::intersection(multi, clippingBox, out);
if (boost::size(out)>0) {
shpMemTiles.StoreShapefileGeometry(layerNum, layer.name, POLYGON_, out, layer.indexed, hasName, name, minzoom, attrIdx);
}

SHPClose(shp);
DBFClose(dbf);
} else {
// Not supported
cerr << "Shapefile entity #" << shape->nShapeId << " type " << shapeType << " not supported" << endl;
}
}
3 changes: 2 additions & 1 deletion src/tilemaker.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -346,7 +346,8 @@ int main(int argc, char* argv[]) {
readShapefile(clippingBox,
layers,
config.baseZoom, layerNum,
shpMemTiles, osmLuaProcessing);
threadNum,
shpMemTiles, osmLuaProcessing);
}
}
shpMemTiles.reportSize();
Expand Down
Loading