使用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文件:

posted @ 2019-09-14 21:33  charlee44  阅读(1190)  评论(0编辑  收藏