基于gdal的图形切分算法(线切面)

多条线切分面(包括折线)参数格式:

{ “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,&currentPoint);
				newLine.addPoint(&currentPoint);
				if (currentPoint == *splitPoint)
				{
					lines.addGeometry(newLine.clone());
					newLine = OGRLineString();
					newLine.addPoint(&currentPoint);
				}
			}
			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;
	}
}

猜你喜欢

转载自blog.csdn.net/u012453032/article/details/82886431