基于GDAL的OGRPolygon网格化
在使用GDAL做开发的过程中，我需要对有的面进行网格化，如：建筑物面等；在shp文件中建筑物面都是以多边形的形式进行描述的，使用GDAL读取SHP文件中的建筑物面会得到一个个OGRPolygon对象，依据这个对象进行网格化。
需要注意的是，我的这种网格化的方法只针对平面坐标系统有效果，我是基于QT5进行开发的。
首先第一步需要根据多边形建立最大矩形，所以我们需要获得该多边形的最大坐标和最小坐标，从而建立一个最大矩形。
1.获取多边形的外环
OGRLinearRing *poLR = mPolygon->getExteriorRing();2.分别获取横纵坐标列表
QList<double> x;
QList<double> y;
for (int i = 0; i < poLR->getNumPoints(); i++)
{
OGRPoint tempPoint;
poLR->getPoint(i, &tempPoint);
x.append(tempPoint.getX());
y.append(tempPoint.getY());
}3.给横纵坐标列表进行排序，然后获取最大和最小值
qSort(x.begin(), x.end());
qSort(y.begin(), y.end());
double Xmin = x.first();
double Ymin = y.first();
double Xmax = x.last();
double Ymax = y.last();然后，根据最大矩形，以给定的长宽进行网格化
//网格化最大矩形
QList<double> girdX;
QList<double> girdY;
double gx, gy;
for (gx = Xmin; gx <= Xmax; gx += girdWidth)
{
girdX.append(gx);
}
gx -= girdWidth;
if (gx < Xmax) girdX.append(Xmax);
for (gy = Ymin; gy <= Ymax; gy += girdHeight)
{
girdY.append(gy);
}
gy -= girdHeight;
if (gy < Ymax) girdY.append(Ymax);最后，实例化网格列表，并删除无用的网格
for (int xIndex = 1; xIndex < girdX.size(); xIndex++)
{
for (int yIndex = 1; yIndex < girdY.size(); yIndex++)
{
//OGRPoint minPoint(girdX[xIndex-1],girdY[yIndex-1]);
//OGRPoint maxPoint(girdX[xIndex],girdY[yIndex]);
//qDebug() << xIndex - 1 << yIndex - 1 << xIndex << yIndex;
OGRPoint lbottomPoint(girdX[xIndex - 1], girdY[yIndex - 1]);
OGRPoint ltopPoint(girdX[xIndex - 1], girdY[yIndex]);
OGRPoint rtopPoint(girdX[xIndex], girdY[yIndex]);
OGRPoint rbottomPoint(girdX[xIndex], girdY[yIndex-1]);
OGRPolygon *tempPolygon = new OGRPolygon; //多边形
OGRLinearRing *mLineRing = new OGRLinearRing; //多边形中的ring
mLineRing->addPoint(&lbottomPoint);
mLineRing->addPoint(<opPoint);
mLineRing->addPoint(&rtopPoint);
mLineRing->addPoint(&rbottomPoint);
tempPolygon->addRingDirectly(mLineRing);
tempPolygon->closeRings();
//筛选网格，求交叉
OGRPolygon *iPolygon =(OGRPolygon*)mPolygon->Intersection(tempPolygon);
OGRLinearRing *plRing = iPolygon->getExteriorRing();
if (plRing != nullptr)
{
double mArea = iPolygon->get_Area();
qDebug() << mArea;
int num = plRing->getNumPoints();
QPolygonF mPolygonF;
for (int i = 0; i < num; i++)
{
OGRPoint mPoint;
plRing->getPoint(i, &mPoint);
QPointF mPointf(mPoint.getX(), mPoint.getY());
mPolygonF.append(mPointf);
}
tempGirdList.append(mPolygonF);
//qDebug() << mPolygonF;
}
}
}网格化示例图
以下是全部代码示例
QList<QPolygonF> LayerOper::girdPolygon(OGRPolygon *mPolygon, double girdWidth, double girdHeight)
{
//获取最小和最大坐标，建立最大矩形
QList<QPolygonF> tempGirdList;
OGRLinearRing *poLR = mPolygon->getExteriorRing();
QList<double> x;
QList<double> y;
for (int i = 0; i < poLR->getNumPoints(); i++)
{
OGRPoint tempPoint;
poLR->getPoint(i, &tempPoint);
x.append(tempPoint.getX());
y.append(tempPoint.getY());
}
qSort(x.begin(), x.end());
qSort(y.begin(), y.end());
double Xmin = x.first();
double Ymin = y.first();
double Xmax = x.last();
double Ymax = y.last();
//网格化最大矩形
QList<double> girdX;
QList<double> girdY;
double gx, gy;
for (gx = Xmin; gx <= Xmax; gx += girdWidth)
{
girdX.append(gx);
}
gx -= girdWidth;
if (gx < Xmax) girdX.append(Xmax);
for (gy = Ymin; gy <= Ymax; gy += girdHeight)
{
girdY.append(gy);
}
gy -= girdHeight;
if (gy < Ymax) girdY.append(Ymax);
for (int xIndex = 1; xIndex < girdX.size(); xIndex++)
{
for (int yIndex = 1; yIndex < girdY.size(); yIndex++)
{
//OGRPoint minPoint(girdX[xIndex-1],girdY[yIndex-1]);
//OGRPoint maxPoint(girdX[xIndex],girdY[yIndex]);
//qDebug() << xIndex - 1 << yIndex - 1 << xIndex << yIndex;
OGRPoint lbottomPoint(girdX[xIndex - 1], girdY[yIndex - 1]);
OGRPoint ltopPoint(girdX[xIndex - 1], girdY[yIndex]);
OGRPoint rtopPoint(girdX[xIndex], girdY[yIndex]);
OGRPoint rbottomPoint(girdX[xIndex], girdY[yIndex-1]);
OGRPolygon *tempPolygon = new OGRPolygon; //多边形
OGRLinearRing *mLineRing = new OGRLinearRing; //多边形中的ring
mLineRing->addPoint(&lbottomPoint);
mLineRing->addPoint(<opPoint);
mLineRing->addPoint(&rtopPoint);
mLineRing->addPoint(&rbottomPoint);
tempPolygon->addRingDirectly(mLineRing);
tempPolygon->closeRings();
//筛选网格，求交叉
OGRPolygon *iPolygon =(OGRPolygon*)mPolygon->Intersection(tempPolygon);
OGRLinearRing *plRing = iPolygon->getExteriorRing();
if (plRing != nullptr)
{
double mArea = iPolygon->get_Area();
qDebug() << mArea;
int num = plRing->getNumPoints();
QPolygonF mPolygonF;
for (int i = 0; i < num; i++)
{
OGRPoint mPoint;
plRing->getPoint(i, &mPoint);
QPointF mPointf(mPoint.getX(), mPoint.getY());
mPolygonF.append(mPointf);
}
tempGirdList.append(mPolygonF);
//qDebug() << mPolygonF;
}
}
}
return tempGirdList;
}经过测试，我发现了一个BUG,多怪我考虑得不够细致，按照上面的多边形划分的方法，可能会出现下面这种情况：
如红色标记的网格，它与多边形的交叉部分有两个区域，所以上面的写法就不正确了，那只获取了一个多边形的区域，要获取多个交叉区域怎么办呢？看下面代码：
OGRGeometry *pRet = mPolygon->Intersection(tempPolygon);
if (pRet != NULL)
{
OGRwkbGeometryType pGeoType = pRet->getGeometryType();
if (pGeoType == wkbPolygon)//这里就是多边形判断
{
OGRPolygon *rdPolygon = (OGRPolygon*)pRet->clone();
QPolygonF mPolygonF;
bool ret = ogrPolygon2qpolygonf(rdPolygon, mPolygonF);
if (ret)
{
tempGirdList.append(mPolygonF);
}
}
else if (pGeoType == wkbMultiPolygon) //这里是多个相交区域获取方法
{
OGRMultiPolygon * multiPolygon = (OGRMultiPolygon *)pRet;
int num = multiPolygon->getNumGeometries();
for (int i = 0; i < num; i++)
{
OGRPolygon *mOGRPolygon = (OGRPolygon *)multiPolygon->getGeometryRef(i);
QPolygonF mPolygonF;
bool ret = ogrPolygon2qpolygonf(mOGRPolygon, mPolygonF);
if (ret)
{
tempGirdList.append(mPolygonF);
}
}
}
}我把OGRPolygon转QPolygon类型封装成了一个函数：
bool LayerOper::ogrPolygon2qpolygonf(OGRPolygon *mPolygon, QPolygonF &targetPolygon)
{
OGRLinearRing *plRing = mPolygon->getExteriorRing();
if (plRing != nullptr)
{
int num = plRing->getNumPoints();
for (int i = 0; i < num; i++)
{
OGRPoint mPoint;
plRing->getPoint(i, &mPoint);
QPointF mPointf(mPoint.getX(), mPoint.getY());
targetPolygon.append(mPointf);
}
return true;
}
else
{
return false;
}
}欢迎大家指正！
转载自：https://blog.csdn.net/octdream/article/details/77097909
