多条线切分面(包括折线)参数格式:
{ “type”: “Polygon”, “coordinates”: [ [ [ -0.581295811920313, 0.117927816793268 ], [ -0.16443779304273, 0.120529114882988 ], [ -0.330270546262362, 0.037287576011957 ], [ -0.120866050039925, -0.048555260948793 ], [ -0.365388070473577, -0.026444227186176 ], [ -0.362136447861428, -0.215038338690854 ], [ -0.535122770827788, -0.217639636780574 ], [ -0.487649080690403, -0.090176030384308 ], [ -0.576093215740873, -0.088875381339448 ], [ -0.581295811920313, 0.117927816793268 ] ] ] }
{ “type”: “LineString”, “coordinates”: [ [ -0.421315979402551, 0.1881628652157 ], [ -0.538374393439938, 0.059398609774575 ], [ -0.446028311254888, 0.056797311684855 ], [ -0.596903600458631, -0.087574732294589 ] ] }
{ “type”: “LineString”, “coordinates”: [ [ -0.17159136278946, 0.197267408529719 ], [ -0.486998756167974, -0.255358459081509 ] ] }
验证结果
{ “type”: “Polygon”, “coordinates”: [ [ [ -0.484617572197741, 0.118531113140991 ], [ -0.538374393439938, 0.059398609774575 ], [ -0.446028311254888, 0.056797311684855 ], [ -0.576614356781295, -0.068160024982656 ], [ -0.581295811920313, 0.117927816793268 ], [ -0.484617572197741, 0.118531113140991 ] ] ] }
{ “type”: “Polygon”, “coordinates”: [ [ [ -0.576614356781295, -0.068160024982656 ], [ -0.446028311254888, 0.056797311684855 ], [ -0.538374393439938, 0.059398609774575 ], [ -0.484617572197741, 0.118531113140991 ], [ -0.225330393185215, 0.120149129858698 ], [ -0.257680635550608, 0.073724864526092 ], [ -0.330270546262362, 0.037287576011957 ], [ -0.293558696408198, 0.022237997810871 ], [ -0.329729305025452, -0.029668690019251 ], [ -0.365388070473577, -0.026444227186176 ], [ -0.364472842373045, -0.079527457017034 ], [ -0.459926844917272, -0.216508870977408 ], [ -0.535122770827788, -0.217639636780574 ], [ -0.487649080690403, -0.090176030384308 ], [ -0.576093215740873, -0.088875381339448 ], [ -0.576614356781295, -0.068160024982656 ] ] ] }
{ “type”: “Polygon”, “coordinates”: [ [ [ -0.257680635550608, 0.073724864526092 ], [ -0.225330393185215, 0.120149129858698 ], [ -0.16443779304273, 0.120529114882988 ], [ -0.257680635550608, 0.073724864526092 ] ] ] }
{ “type”: “Polygon”, “coordinates”: [ [ [ -0.329729305025452, -0.029668690019251 ], [ -0.293558696408198, 0.022237997810871 ], [ -0.120866050039925, -0.048555260948793 ], [ -0.329729305025452, -0.029668690019251 ] ] ] }
{ “type”: “Polygon”, “coordinates”: [ [ [ -0.459926844917272, -0.216508870977408 ], [ -0.364472842373045, -0.079527457017034 ], [ -0.362136447861428, -0.215038338690854 ], [ -0.459926844917272, -0.216508870977408 ] ] ] }
std::vector<std::string> TDVectorCalculation::splitePolygonByMultipleLine(std::string geoPolygon, std::vector<std::string> singlelinevec)
{
auto pgeom = OGRGeometryFactory::createFromGeoJson(geoPolygon.c_str())->toPolygon();
std::vector<OGRGeometry *> ogrfirst, armgvec, tempvec;
std::vector<OGRPoint> ogrfirstPoint, armfirstPoint;
ogrfirst.push_back(pgeom);
std::for_each(singlelinevec.begin(), singlelinevec.end(), [&](std::string & linestr) {
armgvec.clear();
std::for_each(ogrfirst.begin(), ogrfirst.end(), [&](OGRGeometry * ppolygon) {
auto pgeomline = OGRGeometryFactory::createFromGeoJson(linestr.c_str())->toLineString();
tempvec.clear();
splitGeometry(*pgeomline, tempvec,false, armfirstPoint,ppolygon);
if (tempvec.size() == 0)
{
armgvec.push_back(ppolygon);
}
armgvec.insert(armgvec.end(), tempvec.begin(), tempvec.end());
auto pmutilePoint = ppolygon->getExteriorRing()->Difference(pgeomline);
//auto pmutilePoint = (OGRLineString *)OGR_G_Intersection(&ppolygon, pgeomline);
//auto hh = OGR_G_ExportToJson(pmutilePoint);
//if (NULL == pmutilePoint || pmutilePoint->getNumPoints() != 2)
//{
// //求交点如果没有焦点舍弃
// return;
//}
//int num = 0;
//OGRLinearRing plinearRing1, plinearRing2;
//OGRPolygon polygon1, polygon2;
//for (int i = 0; i < ppolygon.getExteriorRing()->getNumPoints() - 1; i++)
//{
// OGRPoint beginpoint, endPoint;
// ppolygon.getExteriorRing()->getPoint(i, &beginpoint);
// ppolygon.getExteriorRing()->getPoint(i + 1, &endPoint);
// OGRLineString templine;
// templine.addPoint(&beginpoint);
// templine.addPoint(&endPoint);
// auto point = (OGRGeometry *)OGR_G_Intersection(&templine, pmutilePoint);
// if (num == 0 && point->getGeometryType() != OGRwkbGeometryType::wkbPoint)
// {
// plinearRing1.setPoint(plinearRing1.getNumPoints(), &beginpoint);
// continue;
// }
// if (num == 0 && point->getGeometryType() == OGRwkbGeometryType::wkbPoint)
// {
// plinearRing1.setPoint(plinearRing1.getNumPoints(), &beginpoint);
// plinearRing1.setPoint(plinearRing1.getNumPoints(), (OGRPoint *)point);
// plinearRing2.setPoint(plinearRing2.getNumPoints(), (OGRPoint *)point);
// num++;
// continue;
// }
// if (num == 1 && point->getGeometryType() != OGRwkbGeometryType::wkbPoint)
// {
// plinearRing2.setPoint(plinearRing2.getNumPoints(), &beginpoint);
// continue;
// }
// if (num == 1 && point->getGeometryType() == OGRwkbGeometryType::wkbPoint)
// {
// plinearRing2.setPoint(plinearRing2.getNumPoints(), &beginpoint);
// plinearRing2.setPoint(plinearRing2.getNumPoints(), (OGRPoint *)point);
// plinearRing1.setPoint(plinearRing1.getNumPoints(), (OGRPoint *)point);
// num++;
// continue;
// }
// if (num == 2 && point->getGeometryType() != OGRwkbGeometryType::wkbPoint)
// {
// plinearRing1.setPoint(plinearRing1.getNumPoints(), &beginpoint);
// num++;
// continue;
// }
//}
//plinearRing1.closeRings();
//plinearRing2.closeRings();
//polygon1.addRing(&plinearRing1);
//polygon2.addRing(&plinearRing2);
//armgvec.push_back(polygon1);
//armgvec.push_back(polygon2);
});
//ogrfirst.add(armgvec);
ogrfirst = armgvec;
ogrfirstPoint.insert(ogrfirstPoint.end(), armfirstPoint.begin(), armfirstPoint.end());
});
std::vector<std::string> tempstring;
std::for_each(ogrfirst.begin(), ogrfirst.end(), [&](OGRGeometry * ppolygon) {
tempstring.push_back(OGR_G_ExportToJson(ppolygon));
});
return tempstring;
}
bool TDVectorCalculation::topologicalTestPointsSplit(const OGRLineString *splitLine, std::vector<OGRPoint> &testPoints, OGRGeometry * geoms)
{
testPoints.clear();
auto intersectionGeom = geoms->Intersection(splitLine);
if (!intersectionGeom)
return false;
bool simple = false;
int nIntersectGeoms = 1;
if (intersectionGeom->getGeometryType() == OGRwkbGeometryType::wkbLineString
|| intersectionGeom->getGeometryType() == OGRwkbGeometryType::wkbPoint)
simple = true;
if (!simple)
nIntersectGeoms = intersectionGeom->toGeometryCollection()->getNumGeometries();//GEOSGetNumGeometries_r(geosinit.ctxt, intersectionGeom.get());
for (int i = 0; i < nIntersectGeoms; ++i)
{
const OGRGeometry *currentIntersectGeom = nullptr;
if (simple)
currentIntersectGeom = intersectionGeom;
else
currentIntersectGeom = intersectionGeom->toGeometryCollection()->getGeometryRef(i);// GEOSGetGeometryN_r(geosinit.ctxt, intersectionGeom.get(), i);
unsigned int sequenceSize = 0;
if (currentIntersectGeom->toPoint()->getGeometryType() == wkbPoint)
{
testPoints.push_back(*currentIntersectGeom->toPoint());
}
if (currentIntersectGeom->toLineString()->getGeometryType() == wkbLineString)
{
auto plinestring = currentIntersectGeom->toLineString();
for (int i = 0; i < plinestring->getNumPoints(); i++)
{
OGRPoint tempPoint;
plinestring->getPoint(i, &tempPoint);
testPoints.push_back(tempPoint);
}
}
}
return true;
}
OGRGeometry * TDVectorCalculation::linePointDifference(OGRGeometry *GEOSsplitPoint, OGRGeometry *mGeos)
{
int type = mGeos->getGeometryType();
OGRMultiCurve * multiCurve;
if (type == wkbMultiLineString)
{
multiCurve = mGeos->clone()->toMultiCurve();
}
else if (type == GEOS_LINESTRING)
{
multiCurve = new OGRMultiCurve();
multiCurve->addGeometry(mGeos->clone());
}
else
{
return nullptr;
}
if (!multiCurve)
{
return nullptr;
}
OGRGeometry * splitGeom = GEOSsplitPoint->clone();
OGRPoint *splitPoint = splitGeom->toPoint();
if (!splitPoint)
{
return nullptr;
}
OGRMultiCurve lines;
//For each part
for (int i = 0; i < multiCurve->getNumGeometries(); ++i)
{
const OGRLineString *line = dynamic_cast<const OGRLineString *>(multiCurve->getGeometryRef(i));
if (line)
{
//For each segment
OGRLineString newLine;
OGRPoint ogrpoint;
line->getPoint(0,&ogrpoint);
newLine.addPoint(&ogrpoint);
int nVertices = line->getNumPoints();
for (int j = 1; j < (nVertices - 1); ++j)
{
OGRPoint currentPoint;
line->getPoint(j,¤tPoint);
newLine.addPoint(¤tPoint);
if (currentPoint == *splitPoint)
{
lines.addGeometry(newLine.clone());
newLine = OGRLineString();
newLine.addPoint(¤tPoint);
}
}
line->getPoint(nVertices - 1, &ogrpoint);
newLine.addPoint(&ogrpoint);
lines.addGeometry(newLine.clone());
}
}
return lines.clone();
}
void TDVectorCalculation::splitLinearGeometry(const OGRLineString &splitLine, std::vector<OGRGeometry *> &newGeometries, OGRGeometry * geoms)
{
if (!geoms->Intersects(&splitLine))
return;
int splitGeomType = splitLine.getGeometryType();// GEOSGeomTypeId_r(geosinit.ctxt, splitLine);
OGRGeometry * splitGeom = NULL;
if (splitGeomType == wkbPoint)
{
splitGeom = linePointDifference(splitLine.clone(),geoms);
}
else
{
splitGeom = geoms->Difference(&splitLine);
}
std::vector<OGRGeometry *> lineGeoms;
int splitType = splitGeom->getGeometryType();// GEOSGeomTypeId_r(geosinit.ctxt, splitGeom.get());
if (splitType == wkbMultiLineString)
{
int nGeoms = splitGeom->toGeometryCollection()->getNumGeometries();// GEOSGetNumGeometries_r(geosinit.ctxt, splitGeom.get());
lineGeoms.reserve(nGeoms);
for (int i = 0; i < nGeoms; ++i)
lineGeoms.push_back(splitGeom->toGeometryCollection()->getGeometryRef(i)->clone());//GEOSGeom_clone_r(geosinit.ctxt, GEOSGetGeometryN_r(geosinit.ctxt, splitGeom.get(), i));
}
else
{
lineGeoms.push_back(splitGeom->clone());
}
mergeGeometriesMultiTypeSplit(lineGeoms,geoms);
for (int i = 0; i < lineGeoms.size(); ++i)
{
newGeometries.push_back(lineGeoms[i]->clone());
delete lineGeoms[i];
lineGeoms[i] = NULL;
}
}
OGRGeometry * TDVectorCalculation::createGeosCollection(int typeId, const std::vector<OGRGeometry *> &geoms)
{
int nNullGeoms = count(geoms.begin(), geoms.end(), nullptr);
int nNotNullGeoms = geoms.size() - nNullGeoms;
int i = 0;
OGRGeometryCollection * geom = new OGRGeometryCollection();
std::vector<OGRGeometry *>::const_iterator geomIt = geoms.cbegin();
for (; geomIt != geoms.cend(); ++geomIt)
{
if (*geomIt)
{
geom->addGeometry(*geomIt);
++i;
}
}
return geom;
}
int TDVectorCalculation::mergeGeometriesMultiTypeSplit(std::vector<OGRGeometry *> &splitResult,OGRGeometry * mGeos)
{
if (!mGeos)
return 1;
//convert mGeos to geometry collection
int type = mGeos->getGeometryType();// GEOSGeomTypeId_r(geosinit.ctxt, mGeos.get());
if (type != wkbGeometryCollection &&
type != wkbMultiLineString &&
type != wkbMultiPolygon &&
type != wkbMultiPoint)
return 0;
std::vector<OGRGeometry *> copyList = splitResult;
splitResult.clear();
//collect all the geometries that belong to the initial multifeature
std::vector<OGRGeometry *> unionGeom;
for (int i = 0; i < copyList.size(); ++i)
{
//is this geometry a part of the original multitype?
bool isPart = false;
for (int j = 0; j < mGeos->toGeometryCollection()->getNumGeometries();j++)
{
if (copyList[i]->Equals(mGeos->toGeometryCollection()->getGeometryRef(i)))//GEOSEquals_r(geosinit.ctxt, copyList[i], GEOSGetGeometryN_r(geosinit.ctxt, mGeos.get(), j)))
{
isPart = true;
break;
}
}
if (isPart)
{
unionGeom.push_back(copyList[i]);
}
else
{
std::vector<OGRGeometry *> geomVector;
geomVector.push_back(copyList[i]);
if (type == GEOS_MULTILINESTRING)
splitResult.push_back(createGeosCollection(GEOS_MULTILINESTRING, geomVector));
else if (type == GEOS_MULTIPOLYGON)
splitResult.push_back(createGeosCollection(GEOS_MULTIPOLYGON, geomVector));
else
//GEOSGeom_destroy_r(geosinit.ctxt, copyList[i]);
delete copyList[i];
}
}
if (!unionGeom.empty())
{
if (type == GEOS_MULTILINESTRING)
splitResult.push_back(createGeosCollection(GEOS_MULTILINESTRING, unionGeom));
else if (type == GEOS_MULTIPOLYGON)
splitResult.push_back(createGeosCollection(GEOS_MULTIPOLYGON, unionGeom));
}
else
{
unionGeom.clear();
}
return 0;
}
void TDVectorCalculation::splitGeometry(const OGRLineString &splitLine, std::vector<OGRGeometry *>&newGeometries, bool topological, std::vector<OGRPoint> &topologyTestPoints, OGRGeometry * geoms)
{
if (splitLine.getNumPoints() < 1 || splitLine.getNumPoints() < 2)
{
return;
}
newGeometries.clear();
OGRGeometry * splitLineGeos = NULL;
if (splitLine.getNumPoints() > 1)
{
splitLineGeos = splitLine.clone();
}
else if (splitLine.getNumPoints() == 1)
{
splitLineGeos = new OGRPoint(splitLine.getX(0), splitLine.getY(0));
}
else
{
return;
}
if (!splitLineGeos->IsValid() || !splitLineGeos->IsSimple())
{
return;
}
if (topological)
{
//find out candidate points for topological corrections
if (!topologicalTestPointsSplit(&splitLine, topologyTestPoints, geoms))
{
return; // TODO: is it really an invalid input?
}
}
//call split function depending on geometry type
if (geoms->getDimension() == 1)
{
splitLinearGeometry(*splitLineGeos->clone()->toLineString(),newGeometries,geoms);
}
else if (geoms->getDimension() == 2)
{
splitPolygonGeometry(splitLineGeos->clone(), newGeometries, geoms);
}
else
{
return;
}
}
OGRGeometry* TDVectorCalculation::nodeGeometries(const OGRGeometry *splitLine, const OGRGeometry *geom)
{
if (!splitLine || !geom)
return nullptr;
OGRGeometry * geometryBoundary;
if (geom->getGeometryType() == wkbPolygon || geom->getGeometryType() == wkbMultiPolygon)
geometryBoundary = geom->getBoundary();
else
geometryBoundary = geom->clone();
OGRGeometry * splitLineClone = splitLine->clone();
OGRGeometry * unionGeometry = splitLineClone->Union(geometryBoundary);
return unionGeometry;
}
int TDVectorCalculation::numberOfGeometries(OGRGeometry *g)
{
if (!g)
return 0;
int geometryType = g->getGeometryType();
if (geometryType == wkbPoint || geometryType == wkbLineString || geometryType == wkbLinearRing
|| geometryType == wkbPolygon)
return 1;
//calling GEOSGetNumGeometries is save for multi types and collections also in geos2
return g->toGeometryCollection()->getNumGeometries();
}
void TDVectorCalculation::splitPolygonGeometry(OGRGeometry *splitLine, std::vector <OGRGeometry *> &newGeometries, OGRGeometry *mGeos)
{
if (!splitLine)
return;
if (!mGeos->Intersects(splitLine))
return;
//first union all the polygon rings together (to get them noded, see JTS developer guide)
OGRGeometry * nodedGeometry = nodeGeometries(splitLine, mGeos);
if (!nodedGeometry)
return;
const OGRGeometry *noded = nodedGeometry;
auto polygons = noded->Polygonize();
if (!polygons || numberOfGeometries(polygons) == 0)
{
return;
}
//test every polygon if contained in original geometry
//include in result if yes
std::vector<OGRGeometry *> testedGeometries;
OGRGeometry * intersectGeometry;
//ratio intersect geometry / geometry. This should be close to 1
//if the polygon belongs to the input geometry
for (int i = 0; i < numberOfGeometries(polygons); i++)
{
auto polygon = polygons->toGeometryCollection()->getGeometryRef(i);
intersectGeometry = mGeos->Intersection(polygon);
if (!intersectGeometry)
{
continue;
}
double intersectionArea = OGR_G_GetArea(intersectGeometry);
double polygonArea = OGR_G_GetArea(polygon);;
const double areaRatio = intersectionArea / polygonArea;
if (areaRatio > 0.99 && areaRatio < 1.01)
testedGeometries.push_back(polygon->clone());
}
int nGeometriesThis = numberOfGeometries(mGeos->clone()); //original number of geometries
if (testedGeometries.empty() || testedGeometries.size() == nGeometriesThis)
{
//no split done, preserve original geometry
for (int i = 0; i < testedGeometries.size(); ++i)
{
delete testedGeometries[i];
testedGeometries[i] = NULL;
}
}
mergeGeometriesMultiTypeSplit(testedGeometries,mGeos->clone());
int i;
for (i = 0; i < testedGeometries.size() && testedGeometries[i]->IsValid(); ++i)
;
if (i < testedGeometries.size())
{
for (i = 0; i < testedGeometries.size(); ++i)
{
delete testedGeometries[i];
testedGeometries[i] = NULL;
}
return;
}
for (i = 0; i < testedGeometries.size(); ++i)
{
newGeometries.push_back(testedGeometries[i]->clone());
delete testedGeometries[i];
testedGeometries[i] = NULL;
}
}