使用GDAL/GEOS求面特征的并集

最新推荐文章于 2023-02-21 22:50:55 发布

charlee44

最新推荐文章于 2023-02-21 22:50:55 发布

阅读量699

点赞数

分类专栏： GEOS # GDAL 文章标签： GDAL GEOS Union 求交

本文链接：https://blog.csdn.net/charlee44/article/details/100864218

版权

GDAL 同时被 2 个专栏收录

19 篇文章 6 订阅

订阅专栏

GEOS

1 篇文章 0 订阅

订阅专栏

存在这样一个示例的矢量文件，包含了两个重叠的面特征：
在这里插入图片描述
一个很常见的需求是求取这个矢量中所有面元素的并集，通过GDAL/GEOS很容易实现这个功能，具体代码如下：

#include <iostream>

#include <gdal/ogrsf_frmts.h>

using namespace std;

bool WritePolygon(string filePath, OGRPolygon *pOgrMerged)
{
	//创建
	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);
	
	//创建特征
	{
		OGRFeature *poFeature = new OGRFeature(poLayer->GetLayerDefn());
		poFeature->SetGeometry(pOgrMerged);

		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:/Work/OSGWork/shpTest/test/src.shp";
	GDALDataset *poDS = (GDALDataset*)GDALOpenEx(filePath.c_str(), GDAL_OF_VECTOR, NULL, NULL, NULL);
	if (!poDS)
	{
		printf("无法读取该文件，试检查格式是否正确！");
		return 1;
	}
	if (poDS->GetLayerCount()<1)
	{
		printf("该文件的层数小于1，试检查格式是否正确！");
		return 1;
	}

	OGRLayer  *poLayer = poDS->GetLayer(0); //读取层
	poLayer->ResetReading();

	std::unique_ptr<OGRPolygon> pOgrMerged(new OGRPolygon());

	OGRFeature *poFeature;
	while ((poFeature = poLayer->GetNextFeature()) != NULL)
	{
		//
		OGRGeometry *pGeo = poFeature->GetGeometryRef();
		OGRwkbGeometryType pGeoType = pGeo->getGeometryType();

		if (pGeoType == wkbPolygon)
		{
			OGRPolygon  *pPolygon = (OGRPolygon*)pGeo;
			if (!pPolygon)
			{
				continue;
			}
						
			OGRPolygon* pTemp = static_cast<OGRPolygon*>(pOgrMerged->Union(pPolygon));
			if (pTemp)
			{
				pOgrMerged.reset(pTemp);
			}		
		}

		OGRFeature::DestroyFeature(poFeature);
	}

	GDALClose(poDS);
	poDS = nullptr;
			
	if (pOgrMerged && pOgrMerged->IsValid() && pOgrMerged->getExteriorRing())
	{
		string path = "D:/Work/OSGWork/shpTest/test/dst.shp";
		WritePolygon(path, pOgrMerged.get());
	}

	return 0;
}

在这段代码中，遍历了示例矢量文件中的每个面元素，求取了所有面元素的并集，得到最终一个面元素，并将这个面元素保存成新的矢量文件。这个矢量文件用ArcMap打开显示如下：

在这里插入图片描述