gdal 根据坐标点求服务区(服务区最小单位是米)

(1) 根据shp文件删选出来一坐标点为中心传入长度为半径的polygon所有拐点和道路信息
(2) 通过删选出来的道路信息构建Graph;
(3) 通过删选出来的拐点和中心点来求Graph最短路径
(4) 根据最短路径来求和半径的关系删选出在半径内的所有点
(5) 根据点求出最小凸包（凸包就是服务区范围）
一下是代码声明:

class TDServiceArea :public TDconversionBase
{
public:
	bool process(LayerPathInfoV1 layerpathinfo, LayerPathInfoV1 armlayerpathinfo,std::vector<std::string> vecinpout,std::vector<double> servicescope);
protected:
	void calsingleFeature(OGRFeature * ogrfeature, int flags_cal) override;
	void allcalsingleFeature() override;
private:
	std::map<int, std::vector<OGRPoint*>> m_RangePoint;
	std::map<int, std::vector<OGRFeature*>> m_RangeLine;
	std::vector<OGRPoint*> m_inputInfo;
	std::vector<double> m_params;
	std::map<int, GNMGraph> m_graph;
	std::map<int,std::map<TDOGRPoint, int>> m_mapPointindex;
	std::map<int,std::map<int, TDOGRPoint>> m_mapindexPoint;
	std::map<int,std::vector<OGRLineString>> m_tempOGRLine;
	double calPathLength(GNMPATH indexinfo, std::vector<OGRLineString> tempOGRLine);
	void makegraph();
	int getIndexData(OGRPoint * ppoint, int num);
};

实现:

#include "TDServiceArea.h"
#include "TDAPICommon.h"
#include <algorithm>
#include "TDAPICommon.h"
bool TDServiceArea::process(LayerPathInfoV1 layerpathinfo, LayerPathInfoV1 armlayerpathinfo, std::vector<std::string> vecinpout, std::vector<double> servicescope)
{
	std::for_each(vecinpout.begin(), vecinpout.end(), [&](std::string & ogrstr) {
		auto pgeoms = OGRGeometryFactory::createFromGeoJson(ogrstr.c_str());
		m_inputInfo.push_back(dynamic_cast<OGRPoint*>(pgeoms));
	});
	m_params = servicescope;
	return process_allwrite(layerpathinfo,armlayerpathinfo);
}
void TDServiceArea::calsingleFeature(OGRFeature * ogrfeature, int flags_cal)
{
	IF_NULL_RETRUN(ogrfeature)
	auto pgeoms = ogrfeature->GetGeometryRef()->toLineString();
	IF_NULL_RETRUN(pgeoms)
	int num = 0;
	std::for_each(m_inputInfo.begin(), m_inputInfo.end(), [&](OGRPoint * ppoint) {
		IF_NULL_RETRUN(ppoint)
		bool flags = false;
		for (int i = 0; i < pgeoms->getNumPoints(); i++)
		{
			OGRPoint *temppoint = new OGRPoint();
			pgeoms->getPoint(i,temppoint);
			if(TDAPICommon::computeDistanceBearing(*temppoint, *ppoint, pmoSourceSRS) <= m_params[num])
			{
				m_RangePoint[num].push_back(temppoint);
				flags = true;
			}
			//道路拐点在距离范围内
		}
		if (flags)
		{
			m_RangeLine[num].push_back(ogrfeature);
		}
		num++;
	});
}
void TDServiceArea::allcalsingleFeature()
{
	makegraph();
	std::for_each(m_RangePoint.begin(), m_RangePoint.end(), [&](std::map<int, std::vector<OGRPoint*>>::value_type &pointinfo) {
		auto ogrindex = getIndexData(m_inputInfo[pointinfo.first], pointinfo.first);
		OGRMultiPoint ogrmultipoint;
		std::for_each(pointinfo.second.begin(), pointinfo.second.end(), [&](OGRPoint * ppoint) {
			auto tempindex = getIndexData(ppoint, pointinfo.first);
			auto paths = m_graph[pointinfo.first].DijkstraShortestPath(ogrindex, tempindex);
			if (calPathLength(paths, m_tempOGRLine[pointinfo.first]) < m_params[pointinfo.first])
			{
				ogrmultipoint.addGeometry(ppoint);
			}
		});
		OGRFeature * feature = new OGRFeature(m_ogrfeaturedefn);
		feature->SetGeometry(ogrmultipoint.ConvexHull());
		//feature->SetGeometry(&ogrmultipoint);
		m_cachevector.push_back(feature);
	});
}
double TDServiceArea::calPathLength(GNMPATH indexinfo, std::vector<OGRLineString> tempOGRLine)
{
	double length = 0;
	std::for_each(indexinfo.begin(), indexinfo.end(), [&](EDGEVERTEXPAIR & indexit) {
		if (indexit.first != -1 && indexit.second != -1)
		{
			OGRLineString  ogrlinestring;
			ogrlinestring = tempOGRLine[indexit.second];
			auto geomNum = ogrlinestring.getNumPoints();
			for (int i = 0; i < geomNum -1; i++)
			{
				OGRPoint p1, p2;
				ogrlinestring.getPoint(i, &p1);
				ogrlinestring.getPoint(i+1, &p2);
				length += TDAPICommon::computeDistanceBearing(p1, p2, pmoSourceSRS);
			}
		}
	});
	return length;
}
int TDServiceArea::getIndexData(OGRPoint * ppoint,int num)
{
	int beginFID = 0;
	double beginlangPath = 0.0;
	auto fun2 = [&](double &ogrlang, OGRPoint & it, OGRPoint & it1, int innum, int & autnum)
	{
		if (TDAPICommon::computeDistanceBearing(it, it1, pmoSourceSRS) < ogrlang)
		{
			ogrlang = TDAPICommon::computeDistanceBearing(it, it1, pmoSourceSRS);
			autnum = innum;
		}
	};
	auto fun1 = [&](OGRPoint & it, int num, OGRPoint & beginPoint) {
		fun2(beginlangPath, beginPoint, it, num, beginFID);
	};
	auto fun = [&](OGRPoint indexPoint) {
		int indexnum = 0;
		std::for_each(m_mapindexPoint[num].begin(), m_mapindexPoint[num].end(), [&](std::map<int, TDOGRPoint>::value_type &it) {
			if (indexnum == 0)
			{
				beginlangPath = TDAPICommon::computeDistanceBearing(indexPoint, it.second, pmoSourceSRS);
			}
			fun1(it.second, it.first, indexPoint);
			indexnum++;
		});
	};
	fun(*ppoint);
	return beginFID;
}
void TDServiceArea::makegraph()
{
	int num = 0;
	std::for_each(m_RangeLine.begin(), m_RangeLine.end(), [&](std::map<int, std::vector<OGRFeature*>>::value_type & it) {
		int FDID = 0;
		int lineFDID = 0;
		std::for_each(it.second.begin(), it.second.end(), [&](OGRFeature* feature) {
			OGRMultiLineString  mpl;
			if (feature->GetGeometryRef()->getGeometryType() == wkbMultiLineString)
				mpl = (*feature->GetGeometryRef()->toMultiLineString());
			else if (feature->GetGeometryRef()->getGeometryType() == wkbLineString)
				mpl.addGeometry(feature->GetGeometryRef()->toLineString());
			auto geomNum = mpl.getNumGeometries();
			for (int i = 0; i < geomNum; i++)
			{
				OGRPoint pt; TDOGRPoint pts;
				auto lineString = mpl.getGeometryRef(i)->toLineString();
				auto numpoint = lineString->getNumPoints();
				for (int j = 0; j < numpoint; j++)
				{
					lineString->getPoint(j, &pt);
					pts.setX(pt.getX());
					pts.setY(pt.getY());
					if (m_mapPointindex[num].find(pts) == m_mapPointindex[num].end())
					{
						m_mapPointindex[num][pts] = FDID;
						m_mapindexPoint[num][FDID] = pts;
						m_graph[num].AddVertex(FDID);
						FDID++;
					}
				}
				for (int j = 0; j < numpoint - 1; j++)
				{
					OGRLineString ogrlinestring;
					OGRPoint *pt1 = new OGRPoint(), *pt2 = new OGRPoint();
					lineString->getPoint(j, pt1);
					ogrlinestring.addPoint(pt1);
					lineString->getPoint(j + 1, pt2);
					ogrlinestring.addPoint(pt2);
					m_tempOGRLine[num].push_back(ogrlinestring);
					lineFDID++;
				}
			}
		});
		lineFDID = 0;
		std::for_each(m_tempOGRLine[num].begin(), m_tempOGRLine[num].end(), [&](OGRLineString & ogrlinestring) {
			OGRPoint pt, pt1;
			ogrlinestring.getPoint(0, &pt);
			TDOGRPoint pts;
			pts.setX(pt.getX());
			pts.setY(pt.getY());
			auto beginindex = m_mapPointindex[num][pts];
			ogrlinestring.getPoint(1, &pt1);
			pts.setX(pt1.getX());
			pts.setY(pt1.getY());
			auto endindex = m_mapPointindex[num][pts];
			m_graph[num].AddEdge(lineFDID, beginindex, endindex, true, TDAPICommon::computeDistanceBearing(pt, pt1, pmoSourceSRS), TDAPICommon::computeDistanceBearing(pt, pt1, pmoSourceSRS));
			lineFDID++;
		});
		num++;
	});
}

坐标点椭球距离计算公式:

double TDAPICommon::computeDistanceBearing(const OGRPoint &p1, const OGRPoint &p2, OGRSpatialReference * psrs,double *course1, double *course2)
{
	if (CommonDoubleNear(p1.getX(), p2.getX()) && CommonDoubleNear(p1.getY(), p2.getY()))
		return 0;
	// ellipsoid
	double mSemiMajor = psrs->GetSemiMajor();
	double mSemiMinor = psrs->GetSemiMinor();
	double mInvFlattening = 1 / psrs->GetInvFlattening();

	double p1_lat = DEG2RAD(p1.getY()), p1_lon = DEG2RAD(p1.getX());
	double p2_lat = DEG2RAD(p2.getY()), p2_lon = DEG2RAD(p2.getX());

	double L = p2_lon - p1_lon;
	double U1 = atan((1 - mInvFlattening) * tan(p1_lat));
	double U2 = atan((1 - mInvFlattening) * tan(p2_lat));
	double sinU1 = sin(U1), cosU1 = cos(U1);
	double sinU2 = sin(U2), cosU2 = cos(U2);
	double lambda = L;
	double lambdaP = 2 * M_PI;

	double sinLambda = 0;
	double cosLambda = 0;
	double sinSigma = 0;
	double cosSigma = 0;
	double sigma = 0;
	double alpha = 0;
	double cosSqAlpha = 0;
	double cos2SigmaM = 0;
	double C = 0;
	double tu1 = 0;
	double tu2 = 0;

	int iterLimit = 20;
	while (abs(lambda - lambdaP) > 1e-12 && --iterLimit > 0)
	{
		sinLambda = sin(lambda);
		cosLambda = cos(lambda);
		tu1 = (cosU2 * sinLambda);
		tu2 = (cosU1 * sinU2 - sinU1 * cosU2 * cosLambda);
		sinSigma = sqrt(tu1 * tu1 + tu2 * tu2);
		cosSigma = sinU1 * sinU2 + cosU1 * cosU2 * cosLambda;
		sigma = atan2(sinSigma, cosSigma);
		alpha = asin(cosU1 * cosU2 * sinLambda / sinSigma);
		cosSqAlpha = cos(alpha) * cos(alpha);
		cos2SigmaM = cosSigma - 2 * sinU1 * sinU2 / cosSqAlpha;
		C = mInvFlattening / 16 * cosSqAlpha * (4 + mInvFlattening * (4 - 3 * cosSqAlpha));
		lambdaP = lambda;
		lambda = L + (1 - C) * mInvFlattening * sin(alpha) *
			(sigma + C * sinSigma * (cos2SigmaM + C * cosSigma * (-1 + 2 * cos2SigmaM * cos2SigmaM)));
	}

	if (iterLimit == 0)
		return -1;  // formula failed to converge

	double uSq = cosSqAlpha * (mSemiMajor * mSemiMajor - mSemiMinor * mSemiMinor) / (mSemiMinor * mSemiMinor);
	double A = 1 + uSq / 16384 * (4096 + uSq * (-768 + uSq * (320 - 175 * uSq)));
	double B = uSq / 1024 * (256 + uSq * (-128 + uSq * (74 - 47 * uSq)));
	double deltaSigma = B * sinSigma * (cos2SigmaM + B / 4 * (cosSigma * (-1 + 2 * cos2SigmaM * cos2SigmaM) -
		B / 6 * cos2SigmaM * (-3 + 4 * sinSigma * sinSigma) * (-3 + 4 * cos2SigmaM * cos2SigmaM)));
	double s = mSemiMinor * A * (sigma - deltaSigma);

	if (course1)
	{
		*course1 = atan2(tu1, tu2);
	}
	if (course2)
	{
		// PI is added to return azimuth from P2 to P1
		*course2 = atan2(cosU1 * sinLambda, -sinU1 * cosU2 + cosU1 * sinU2 * cosLambda) + M_PI;
	}
	return s;
}