使用GDAL/OGR讀寫矢量文件

感覺GIS中矢量相關內容還是挺龐雜的，並且由於版本迭代的關系，使用GDAL/OGR讀寫矢量的資料也有點不太一樣。這里總結了一個讀寫矢量的示例，實現代碼如下：

#include <iostream>

#include <gdal/ogrsf_frmts.h>

using namespace std;

bool ReadDXF(string filePath, vector<vector<OGRPoint>>& vertexPoint)
{	
	GDALDataset *poDS = (GDALDataset*)GDALOpenEx(filePath.c_str(), GDAL_OF_VECTOR, NULL, NULL, NULL);
	if (!poDS)
	{
		printf("無法讀取該文件，試檢查格式是否正確！");
		return false;
	}
	if (poDS->GetLayerCount()<1)
	{
		printf("該文件的層數小於1，試檢查格式是否正確！");
		return false;
	}

	OGRLayer  *poLayer = poDS->GetLayer(0); //讀取層
	poLayer->ResetReading();
	
	OGRFeature *poFeature;
	while ((poFeature = poLayer->GetNextFeature()) != NULL)
	{
		OGRGeometry *pGeo = poFeature->GetGeometryRef();
		OGRwkbGeometryType pGeoType = pGeo->getGeometryType();

		if (pGeoType == wkbLineString || pGeoType == wkbLineString25D)
		{
			OGRLinearRing  *pCurve = (OGRLinearRing*)pGeo;
			if (pCurve->getNumPoints() < 1)
			{
				continue;
			}

			vector<OGRPoint> pl;
			for (int i = 0; i<pCurve->getNumPoints(); i++)
			{
				OGRPoint point;
				pCurve->getPoint(i, &point);						
				pl.push_back(point);
			}
			vertexPoint.push_back(pl);
		}

		////		
		//OGRFeatureDefn *poFDefn = poLayer->GetLayerDefn();		
		//int n = poFDefn->GetFieldCount(); //獲得字段的數目，不包括前兩個字段（FID,Shape);
		//for (int iField = 0; iField <n; iField++)
		//{			
		//    //輸出每個字段的值
		//    cout << poFeature->GetFieldAsString(iField) << "    ";			
		//}
		//cout << endl;   

		OGRFeature::DestroyFeature(poFeature);
	}

	GDALClose(poDS);
	poDS = nullptr;

	return true;
}

bool WriteShp(string filePath, vector<vector<OGRPoint>> vertexPoint)
{
	//創建
	GDALDriver* driver = GetGDALDriverManager()->GetDriverByName("ESRI Shapefile");
	if (!driver)
	{
		printf("Get Driver ESRI Shapefile Error！\n");
		return false;
	}
	
	GDALDataset* dataset = driver->Create(filePath.c_str(), 0, 0, 0, GDT_Unknown, NULL);
	OGRLayer* poLayer = dataset->CreateLayer("houseType", NULL, wkbPolygon, NULL);
	
	//創建屬性字段
	{
		// 字符串
		OGRFieldDefn oField1("名稱", OFTString);
		oField1.SetWidth(8);
		if (poLayer->CreateField(&oField1) != OGRERR_NONE) {
			printf("Creating Name field failed.\n"); return FALSE;
		}

		// 浮點數
		OGRFieldDefn oField2("面積", OFTReal);
		oField2.SetPrecision(3);
		if (poLayer->CreateField(&oField2) != OGRERR_NONE) {
			printf("Creating Name field failed.\n"); return FALSE;
		}

		// 整型
		OGRFieldDefn oField3("結點數", OFTInteger);
		if (poLayer->CreateField(&oField3) != OGRERR_NONE) {
			printf("Creating Name field failed.\n"); return FALSE;
		}
	}

	//創建特征
	for (auto& iter : vertexPoint)
	{
		OGRFeature *poFeature = new OGRFeature(poLayer->GetLayerDefn());
		
		OGRLinearRing ogrring;
		int pNum = (int)iter.size();
		ogrring.setNumPoints(pNum);
		for (int i = 0; i < iter.size(); i++)
		{
			ogrring.setPoint(i, iter[i].getX(), iter[i].getY(), iter[i].getZ());
			//cout << iter[i].x() << '\t' << iter[i].y() << '\t' << iter[i].z() << endl;
		}
		//cout << "-----------------------------\n";

		OGRPolygon polygon;
		polygon.addRing(&ogrring);
		poFeature->SetGeometry(&polygon);		

		poFeature->SetField("名稱", "多邊形");
		poFeature->SetField("面積", polygon.get_Area());
		poFeature->SetField("結點數", pNum);

		if (poLayer->CreateFeature(poFeature) != OGRERR_NONE)
		{
			printf("Failed to create feature in shapefile.\n");
			return false;
		}
	}
	
	//釋放
	GDALClose(dataset);
	dataset = nullptr;
	//GDALDestroyDriverManager();

	return true;
}

int main()
{
	GDALAllRegister();
	CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO");  //支持中文路徑
	CPLSetConfigOption("SHAPE_ENCODING", "");  //解決中文亂碼問題
		
	string filePath = "D:/2.dxf";
	vector<vector<OGRPoint>> vertexPoint;
	if (!ReadDXF(filePath, vertexPoint))
	{
		return 1;
	}

	string newPath = "C:/Users/charlee/Desktop/SHP/dst.shp";
	WriteShp(newPath, vertexPoint);

	return 0;
}

在這個示例中，讀取一個DXF文件中的線（環）特征，將其轉換成面，然后保存在一個SHP中。同時，還給該SHP文件寫入了相應的屬性字段。

讀取的DXF文件：

創建並保存的SHP文件：