抽稀点要素,计算角度

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

//抽稀点要素,计算角度
bool VacuateDirection(const QString &pShp0, const QString &pShpRes, const int &interval)
{
	//注册
	GDALAllRegister();
	OGRRegisterAll();

	//读取原图层
	GDALDataset *oDS = (GDALDataset*)GDALOpenEx(pShp0.toStdString().c_str(), GDAL_OF_VECTOR, NULL, NULL, NULL);
	OGRLayer *oLayer = oDS->GetLayer(0);
	OGRSpatialReference *spatialReference = oLayer->GetSpatialRef();
	int featureCount = oLayer->GetFeatureCount();

#pragma region create new shp
	GDALDriver* driver = OGRSFDriverRegistrar::GetRegistrar()->GetDriverByName("ESRI Shapefile");

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

	OGRLayer *newLayer = newDataSet->CreateLayer("vacuatePointDirection", spatialReference, wkbPointZM, NULL);
	if (newLayer == NULL)
		return false;

#pragma region fields
	//FIDD
	OGRFieldDefn fieldFIDD("FIDD", OFTInteger);
	fieldFIDD.SetWidth(32);
	newLayer->CreateField(&fieldFIDD);
	//X
	OGRFieldDefn fieldX("X", OFTReal);
	fieldX.SetPrecision(16);
	newLayer->CreateField(&fieldX);
	//Y
	OGRFieldDefn fieldY("Y", OFTReal);
	fieldY.SetPrecision(16);
	newLayer->CreateField(&fieldY);
	//Z
	OGRFieldDefn fieldZ("Z", OFTReal);
	fieldZ.SetPrecision(16);
	newLayer->CreateField(&fieldZ);
	//direction
	OGRFieldDefn fieldDIRECTION("direction", OFTReal);
	fieldDIRECTION.SetPrecision(16);
	newLayer->CreateField(&fieldDIRECTION);
#pragma endregion

#pragma endregion

	//抽稀 计算角度
	for (int i = 0; i < (featureCount - interval); i += interval)
	{
#pragma region dir
		//P1
		OGRFeature *oFeature1 = oLayer->GetFeature(i);
		int fid1 = oFeature1->GetFID();
		OGRGeometry* oGeometry1 = oFeature1->GetGeometryRef();
		OGRPoint* oPt1 = (OGRPoint*)oGeometry1;
		double x1 = oPt1->getX();
		double y1 = oPt1->getY();
		double z1 = oPt1->getZ();

		//P2
		OGRFeature *oFeature2 = oLayer->GetFeature(i + interval);
		int fid2 = oFeature2->GetFID();
		OGRGeometry* oGeometry2 = oFeature2->GetGeometryRef();
		OGRPoint* oPt2 = (OGRPoint*)oGeometry2;
		double x2 = oPt2->getX();
		double y2 = oPt2->getY();
		double z2 = oPt2->getZ();

		//direction
		//double A = GetAnglePP(x1, y1, x2, y2);
		double ra = std::atan2(y2 - y1, x2 - x1);
		double A = ra * 180 / 3.1415926;
		if (A < 0)
		{
			A += 360;
		}
#pragma endregion

#pragma region create new feature
		//创建要素
		OGRFeature *newFeature = OGRFeature::CreateFeature(newLayer->GetLayerDefn());
		OGRPoint pt;
		pt.setX(x1);
		pt.setY(y1);
		pt.setZ(z1);
		newFeature->SetGeometry(&pt);

		//属性
		newFeature->SetField("FIDD", fid1);
		newFeature->SetField("X", x1);
		newFeature->SetField("Y", y1);
		newFeature->SetField("Z", z1);
		newFeature->SetField("direction", A);

		//判断
		if (newLayer->CreateFeature(newFeature) != OGRERR_NONE)
		{
			//QMessageBox::warning(0, ("prompt"), "Failed to create feature in shapefile!!!");
			return false;
		}
		OGRFeature::DestroyFeature(newFeature);
#pragma endregion
	}
	GDALClose(oDS);

	GDALClose(newDataSet);
}

ArcGIS设置旋转字段