GDAL：创建矢量线、矢量面数据

分享给有需要的人，代码质量勿喷。

C/C++——常规——附加包含目录——添加include目录

链接器——输入——附加依赖项——添加...\gdal_i.lib

#include "ogrsf_frmts.h"

一、创建矢量线图层（wkbLineStringZM） 添加单个要素

void xjCreateVectorLineByGDAL(QList<xjPoint> xjListNode, const QString &xjSavePath)
 {
	 GDALAllRegister();
	 OGRRegisterAll();
	 const char *xjDriverName = "ESRI Shapefile";
	 GDALDriver *xjDriver = OGRSFDriverRegistrar::GetRegistrar()->GetDriverByName(xjDriverName);
	 if (xjDriver == NULL)
	 {
		 QMessageBox::warning(0, ("prompt"), QString(xjDriverName) + "driver not available!!!");
		 return;
	 }

	 GDALDataset *xjDataset = xjDriver->Create(xjSavePath.toStdString().c_str(), 0, 0, 0, GDT_Unknown, NULL);
	 if (xjDataset == NULL)
	 {
		 QMessageBox::warning(0, ("prompt"), "Creation of output file failed!!!");
		 return;
	 }
     
     //空间参考
     OGRSpatialReference sf = CreateProjectionReference(120);

	 OGRLayer *xjLayer = xjDataset->CreateLayer("point_out", sf/* NULL */, wkbLineStringZM, NULL);
	 if (xjLayer == NULL)
	 {
		 QMessageBox::warning(0, ("prompt"), "Layer creation failed!!!");
		 return;
	 }

/* GDALDataset* xjDataset = (GDALDataset*)GDALOpenEx(xjShpPath.toStdString().c_str(), GDAL_OF_UPDATE, NULL, NULL, NULL);
	OGRLayer *xjLayer = xjDataset->GetLayer(0); */

	 //创建属性字段
	 OGRFieldDefn fieldID("ID", OFTInteger);
	 fieldID.SetWidth(32);
	 xjLayer->CreateField(&fieldID);
	 OGRFieldDefn fieldNAME("NAME", OFTString);
	 fieldNAME.SetWidth(32);
	 xjLayer->CreateField(&fieldNAME);
	 OGRFieldDefn fieldLENGTH("LENGTH", OFTReal);
	 fieldLENGTH.SetPrecision(16);
	 xjLayer->CreateField(&fieldLENGTH);

	 //创建要素
	 OGRFeature *xjFeature = OGRFeature::CreateFeature(xjLayer->GetLayerDefn());
	 OGRLineString xjLine;
	 int pointCount = xjListNode.size();
	 for (int i = 0; i < pointCount; i++)
	 {
		 xjPoint Pd = xjListNode.at(i);
		 OGRPoint pt;
		 pt.setX(Pd.x);
		 pt.setY(Pd.y);
		 pt.setZ(Pd.z);
		 xjLine.addPoint(&pt);
	 }
	 xjFeature->SetGeometry(&xjLine);

	 //属性
	 xjFeature->SetField("ID", "id");
	 xjFeature->SetField("NAME", "name");
	 xjFeature->SetField("LENGTH", 2.3);

	 //判断
	 if (xjLayer->CreateFeature(xjFeature) != OGRERR_NONE)
	 {
		 QMessageBox::warning(0, ("prompt"), "Failed to create feature in shapefile!!!");
		 return;
	 }
	 OGRFeature::DestroyFeature(xjFeature);

	 GDALClose(xjDataset);
 }

二、创建矢量线图层（wkbMultiLineStringZM，即合并wkbLineStringZM中的多个要素） 添加单个要素

bool ImageCrackThread::xjCreateCrackShapeFile()
{
	QMutexLocker locker(&xjMutex);

#pragma region create shape file
	QString shpPath = "F:/testRes.shp";

	GDALAllRegister();
	OGRRegisterAll();
	const char *xjDriverName = "ESRI Shapefile";
	GDALDriver *xjDriver = OGRSFDriverRegistrar::GetRegistrar()->GetDriverByName(xjDriverName);
	if (xjDriver == NULL)
		return false;

	GDALDataset *xjDataset = xjDriver->Create(shpPath.toStdString().c_str(), 0, 0, 0, GDT_Unknown, NULL);
	if (xjDataset == NULL)
		return false;

	OGRSpatialReference sf = CreateProjectionReference(paraCamera, false);
	OGRLayer *pLayer = xjDataset->CreateLayer("crack", &sf, wkbMultiLineStringZM, NULL);
	if (pLayer == NULL)
		return false;
#pragma endregion

	OGRFeature *xjFeature = OGRFeature::CreateFeature(pLayer->GetLayerDefn());
	OGRMultiLineString xjMultiLine;
	for (int i = 0; i < 10; i += 2)
	{
		OGRLineString xjLine;
		for (int j = 0; j < 5; j++)
		{
			OGRPoint pt;
			pt.setX(i*j + j);
			pt.setY(i*j + j * j);
			pt.setZ(0);
			xjLine.addPoint(&pt);
		}
		xjMultiLine.addGeometry(&xjLine);
	}

	xjFeature->SetGeometry(&xjMultiLine);

	if (pLayer->CreateFeature(xjFeature) != OGRERR_NONE)
		return false;
	OGRFeature::DestroyFeature(xjFeature);

	GDALClose(xjDataset);

	return true;
}

三、创建矢量面数据（wkbPolygonZM） 多个要素

 void xjCreateVectorPolygonByGDAL(const QList<QList<xjPoint>> &xjListPoint, const QString &xjSavePath)
 {
	 GDALAllRegister();
	 OGRRegisterAll();
	 const char *xjDriverName = "ESRI Shapefile";
	 GDALDriver *xjDriver = OGRSFDriverRegistrar::GetRegistrar()->GetDriverByName(xjDriverName);
	 if (xjDriver == NULL)
	 {
		 QMessageBox::warning(0, "prompt", "error of GDAL: Driver not available!!!");
		 return;
	 }

	 GDALDataset *xjDataset = xjDriver->Create(xjSavePath.toStdString().c_str(), 0, 0, 0, GDT_Unknown, NULL);
	 if (xjDataset == NULL)
	 {
		 QMessageBox::warning(0, "prompt", "error of GDAL: Creation of output file failed!!!");
		 return;
	 }

     //空间参考
     OGRSpatialReference sf = CreateProjectionReference(120);

	 OGRLayer *xjLayer = xjDataset->CreateLayer("polygon1", sf/*NULL*/, wkbPolygonZM, NULL);
	 if (xjLayer == NULL)
	 {
		 QMessageBox::warning(0, "prompt", "error of GDAL: Layer creation failed!!!");
		 return;
	 }

	 //创建属性字段
	 OGRFieldDefn fieldID("ID", OFTInteger);
	 fieldID.SetWidth(32);
	 xjLayer->CreateField(&fieldID);
	 OGRFieldDefn fieldNAME("NAME", OFTString);
	 fieldNAME.SetWidth(32);
	 xjLayer->CreateField(&fieldNAME);
	 OGRFieldDefn fieldAREA("AREA", OFTReal);
	 fieldAREA.SetPrecision(16);
	 xjLayer->CreateField(&fieldAREA);

	 //创建要素
	 for (int i = 0; i < xjListPoint.size(); i++)
	 {
		 //创建要素
		 OGRFeature *xjFeature = OGRFeature::CreateFeature(xjLayer->GetLayerDefn());
		 //矢量面要素的边界是闭合环
		 OGRLinearRing xjRing;
		 QList<xjPoint> list = xjListPoint.at(i);
		 for (int j = 0; j < list.size(); j++)
		 {
			 xjRing.addPoint(list.at(j).x, list.at(j).y, list.at(j).z);
		 }
		 xjRing.closeRings();//首尾点重合形成闭合环

		 //图层添加要素
		 OGRPolygon xjPolygon;
		 xjPolygon.addRing(&xjRing);
		 xjFeature->SetGeometry(&xjPolygon);

		 //设置属性
		 xjFeature->SetFID(i);
		 xjFeature->SetField(0, i);
		 std::string pname = "name_" + to_string(i);
		 xjFeature->SetField(1, pname.c_str());
		 xjFeature->SetField("AREA", 2.2);

		 //判断
		 if (xjLayer->CreateFeature(xjFeature) != OGRERR_NONE)
		 {
			 QMessageBox::warning(0, "prompt", "error of GDAL: Failed to create feature in shapefile!!!");
			 continue;
		 }
		 OGRFeature::DestroyFeature(xjFeature);
	 }
	 GDALClose(xjDataset);
 }

四、自定义投影（也可忽略，自己记录）

/* create projection information */
OGRSpatialReference ImageCrackThread::CreateProjectionReference(const int &cL/* 3-40-120-4528 */)
{
	QString projName = "WGS_1984_GK_Zone_" + QString::number(cL / 3) + "_3";
	QString projInfo = "PROJCS[" + projName +
		", GEOGCS[\"GCS_WGS_1984\", DATUM[\"D_WGS_1984\", SPHEROID[\"WGS_1984\", 6378137.0, 298.257223563]]" +
		", PRIMEM[\"Greenwich\", 0.0], UNIT[\"Degree\", 0.0174532925199433]], PROJECTION[\"Gauss_Kruger\"]" +
		", PARAMETER[\"False_Easting\", 500000.0], PARAMETER[\"False_Northing\", 0.0]" +
		", PARAMETER[\"Central_Meridian\", " + QString::number(cL) + "]" +
		", PARAMETER[\"Scale_Factor\", 1.0], PARAMETER[\"Latitude_Of_Origin\", 0.0], UNIT[\"Meter\", 1.0]" +
		", AUTHORITY[\"EPSG\", " + QString::number(cL / 3 + 4488) + "]]";

	const char *pszWGS_GK_Zone = projInfo.toStdString().c_str();

	OGRSpatialReference sfWGS8484_GK_Zone;
	sfWGS8484_GK_Zone.SetFromUserInput(pszWGS_GK_Zone);
	return sfWGS8484_GK_Zone;
}

五、合并矢量文件（也可忽略，自己记录）

void ImageCrackDetectMulti::MergeShape(const QString &newShpPath, const QString &shpDir)
{
	GDALAllRegister();
	OGRRegisterAll();
	GDALDriver* driver = OGRSFDriverRegistrar::GetRegistrar()->GetDriverByName("ESRI Shapefile");

	//new file
	GDALDataset *newDataSet = driver->Create(newShpPath.toStdString().c_str(), 0, 0, 0, GDT_Unknown, NULL);
	if (newDataSet == NULL)
		return false;

	OGRSpatialReference sf = CreateProjectionReference(120);

	OGRLayer *pLayer = newDataSet->CreateLayer("point_out", &sf, wkbLineStringZM, NULL);
	if (pLayer == NULL)
		return false;

#pragma region fields
	//ID
	OGRFieldDefn fieldID("ID", OFTInteger);
	fieldID.SetWidth(32);
	pLayer->CreateField(&fieldID);
	//长度
	OGRFieldDefn fieldLENGTH("LENGTH(m)", OFTReal);
	fieldLENGTH.SetPrecision(16);
	pLayer->CreateField(&fieldLENGTH);
	//程度
	OGRFieldDefn fieldEXTENT("EXTENT", OFTString);
	fieldEXTENT.SetWidth(32);
	pLayer->CreateField(&fieldEXTENT);
#pragma endregion

	//filter files
	QFileInfoList fiList = FilterFileFromDir(shpDir, "*.shp");

	for (int j = 0; j < fiList.size(); j++)
	{
		double shpLength = 0;
		GDALDataset *poDS = (GDALDataset*)GDALOpenEx((fiList.at(j).absoluteFilePath()).toStdString().c_str(),
			GDAL_OF_VECTOR, NULL, NULL, NULL);
		OGRLayer *tmpLayer = poDS->GetLayer(0);
		OGRFeature *pFeature;
		while ((pFeature = tmpLayer->GetNextFeature()) != NULL)
		{
			pLayer->CreateFeature(pFeature);

			//OGRGeometry* pGeometry = pFeature->GetGeometryRef();
			//OGRLineString* pLineGeo = (OGRLineString*)pGeometry;
			//double staX = pLineGeo->getX(0);
			//double staY = pLineGeo->getY(0);
			shpLength += ((OGRLineString*)(pFeature->GetGeometryRef()))->get_Length();
		}
		GDALClose(poDS);
	}
	GDALClose(newDataSet);
}

六、从文件夹中，根据扩展名，筛选所需文件（也可忽略，自己记录）

//filter file from directory 
QFileInfoList ImageCrackDetectMulti::FilterFileFromDir(const QString &directoryPath, const QString &filterName)
{
	QStringList filters;
	filters << filterName;

	QDir dir;
	dir.setPath(directoryPath);
	dir.setNameFilters(filters);

	QList<QFileInfo> fileInfo = QList<QFileInfo>(dir.entryInfoList(filters, QDir::AllEntries, QDir::NoSort));

	QDirIterator iter(dir, QDirIterator::Subdirectories);
	QFileInfoList fiList;
	while (iter.hasNext())
	{
		iter.next();
		QFileInfo info = iter.fileInfo();
		if (info.isFile())
		{
			fiList.append(info);
		}
	}

	return fiList;
}