基于方格网法的填挖方量计算(C++)

目录

前言

一、功能设计

二、关键步骤

1.划分格网以及格网是否处于区域内

三、具体效果

总结

---

前言

填挖方介绍：填方指的是路基表面高于原地面时，从原地面填筑至路基表面部分的土石体积。土木工程施工时向地基或其他地方填充的土石方；挖方是指路基表面低于原地面时，从原地面至路基表面挖去部分的土石体积。

---

一、功能设计

引入微积分的概念，把填挖方区域按经纬度方向分成若干个微小的正方形格子，只要格子足够小，就可以不用考虑格子内自身的高程差，认为格子内所有点处于同一高程，此时只需要用格子的面积，乘以格子内任一点的高程到参考高程（或者是设计标高）的距离，即为该格子的填挖方量，正值为挖方，负值为填方，为0则不挖不填；对于边界情况，只需要判断你取的点是否在区域内即可，例如你去所有格子左下角点和格子中心点计算结果会存在差异。

二、关键步骤 1.划分格网以及格网是否处于区域内

1.首先获取多边形区域的外接矩形（长宽平行于经纬度）；

2.按格网边长进行格网划分；

注：有两种方案：第一种为从左下角坐标开始依次累加单个格网边长的距离，直到坐标值超过右上角坐标；第二种为先算出区域包围盒的长和宽，用二分法（即长宽各自循环除以2）直到长度小于格网边长；无论哪种方法都需要考虑投影问题，即地球表面是曲面，计算长度时需要进行改算。

3.判断格网（格网用于计算高程的点）是否处于区域内

void CEarthworkCompute::Form3DPoints()
{
    //vector points:表示闭合区域由这些点围成
    //计算外接矩形
	double minX = m_vec3D.at(0).x;
	double maxX = m_vec3D.at(0).x;
	double minY = m_vec3D.at(0).y;
	double maxY = m_vec3D.at(0).y;
	for ( unsigned int i = 1 ; i < m_vec3D.size() ; i++ )
	{
		CVector3D q = m_vec3D.at(i);
		minX = qMin( q.x, minX );
		maxX = qMax( q.x, maxX );
		minY = qMin( q.y, minY );
		maxY = qMax( q.y, maxY );
	}

	//计算横向距离
	CVector3D vLeftLowPoint(minX, minY, 0.0);
	CVector3D vRightLowPoint(maxX, minY, 0.0);
	if (m_pCoordConvertDistance == NULL)
	{
		return;
	}
	//CCoordConvertDistance* CoordConvertDistance = new CCoordConvertDistance(NULL);构造函数中new具体实现见另一篇博文根据地球表面两点坐标（经纬度）计算两点间距离
	m_pCoordConvertDistance->ConvertFunc(vLeftLowPoint, vRightLowPoint);
	double dX = m_pCoordConvertDistance->GetDistance();

	//计算纵向距离
	CVector3D vLeftTopPoint(minX, maxY, 0.0);
	m_pCoordConvertDistance->ConvertFunc(vLeftLowPoint, vLeftTopPoint);
	double dY = m_pCoordConvertDistance->GetDistance();

	//计算实际步长以及格网数量
	int nXcount = 0;
	int nYcount = 0;
	double dXTempDistance = dX;
	double dYTempDistance = dY;
    //按二分法划分格网
	//经度差
	while(dXTempDistance > m_designStep)
	{
		if (m_bStop)
		{
			return;
		}
		dXTempDistance = dXTempDistance/2;
		nXcount++;
	}
	int nXStepCount = pow(2.0, nXcount);
	double dActXStepDistance = (maxX-minX)/nXStepCount;

	//纬度差
	while(dYTempDistance > m_designStep)
	{
		if (m_bStop)
		{
			return;
		}
		dYTempDistance = dYTempDistance/2;
		nYcount++;
	}
	int nYStepCount = pow(2.0, nYcount);
	double dActYStepDistance = (maxY-minY)/nYStepCount;

	//格网面积
	m_GridArea = dXTempDistance * dYTempDistance;
	IAcMapModule* pIAcMapModule = GetGlobalAcMapModule();
	m_v2DList.clear();
	double dLeftDownX = minX;
	double dLeftDownY = minY;
	int nCount = 0;
	for(int i = 0; iCalculateAltitude(dLeftDownX, dLeftDownY, 15);//计算该点的高程值
				CVector3D *vTemp = new CVector3D(dLeftDownX,dLeftDownY,delevation);
				m_v3DList.push_back(vTemp);//获取处于区域内的带高程的点集
			}
			dLeftDownY += dActYStepDistance;
			nCount++;
			emit SigProgress(0, nCount, nXStepCount * nYStepCount);//进度条
		}
		dLeftDownY = minY;
		dLeftDownX += dActXStepDistance;
	}
}

判断一个点是否在几个点围成的区域内具体实现：

bool CIsPointInPolygon::IsPointInPolygonNew(const CVector3D& p, const std::vector& vec3D)
{

	int count = vec3D.size();

	if(count < 3)
	{
		return false;
	}

	bool result = false;

	for(int i = 0, j = count - 1; i < count; i++)
	{
		CVector3D p1 = vec3D.at(i);
		CVector3D p2 = vec3D.at(j);

		if(p1.x < p.x && p2.x >= p.x || p2.x < p.x && p1.x >= p.x)
		{
			if(p1.y + (p.x - p1.x) / (p2.x - p1.x) * (p2.y - p1.y) < p.y)
			{
				result = !result;
			}
		}
		j = i;
	}
	return result;
}

---

三、具体效果

 

---

总结

今天分享到这里！