OGR spatial analysis result error

[xiaose]

two layers A and B, firstly I intersect A and B , and then SymDifference operation, after the results seem to be wrong,  the same operation.ArcGis‘s Reslut will be What I want。

[xiaose]

AW and BW are input layers 22 is output layer

[Even Rouault]

Care to share also the code that was used ?

[xiaose]

Replying to rouault:

Care to share also the code that was used ?

int AEraseB(const char* A_Path/*输入*/, const char* B_Path/*输入*/, const char* C_Path/*输出*/)
{
	if (A_Path == "" || B_Path == "" || C_Path == "")
	{
		return -1;
	}
	//注册驱动
	OGRRegisterAll();

	CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO");//解决中文路径
	CPLSetConfigOption("SHAPE_ENCODING", "");//解决字段中文乱码

	OGRDataSource *poA_DS = NULL;
	OGRDataSource *poB_DS = NULL;
	OGRDataSource *poC_DS = NULL;
	OGRDataSource *poMem_DS = NULL;

	poA_DS = OGRSFDriverRegistrar::Open(A_Path, FALSE); //A的数据集
	poB_DS = OGRSFDriverRegistrar::Open(B_Path, FALSE); //B的数据集

	if (poA_DS == NULL || poB_DS == NULL)
	{
		return -1;
	}

	OGRLayer  *poLayer_A = NULL;
	OGRLayer  *poLayer_B = NULL;

	poLayer_A = poA_DS->GetLayer(0);
	poLayer_B = poB_DS->GetLayer(0);

	//OGRwkbGeometryType type = poLayer_B->GetGeomType();

	if (poLayer_A == NULL || poLayer_B == NULL)
	{
		OGRDataSource::DestroyDataSource(poA_DS);
		OGRDataSource::DestroyDataSource(poB_DS);
		return -1;
	}
    
	//创建C
	const char *pszDriverName = "ESRI Shapefile";
	const char *pszDriverMem = "Memory";

	OGRSFDriver *poDriver = NULL;
	OGRSFDriver *poDriverMem = NULL;
	poDriverMem = OGRSFDriverRegistrar::GetRegistrar()->GetDriverByName(pszDriverMem);
	poDriver = OGRSFDriverRegistrar::GetRegistrar()->GetDriverByName(pszDriverName);
	if (poDriver == NULL || poDriverMem == NULL)
	{
		OGRDataSource::DestroyDataSource(poA_DS);
		OGRDataSource::DestroyDataSource(poB_DS);
		return -1;
	}
	poDriver->DeleteDataSource(C_Path);
	poC_DS = poDriver->CreateDataSource(C_Path, NULL);
	poMem_DS = poDriverMem->CreateDataSource("", NULL);
	if (poC_DS == NULL || poMem_DS == NULL)
	{
		OGRDataSource::DestroyDataSource(poA_DS);
		OGRDataSource::DestroyDataSource(poB_DS);
		return -1;
	}
	OGRLayer *poLayer_C = NULL;
	OGRLayer *poLaer_Mem = NULL;
	OGRSpatialReference *  poSpatialRef = poLayer_A->GetSpatialRef();

	poLayer_C = poC_DS->CreateLayer("C", poSpatialRef, wkbPolygon, NULL);
	poLaer_Mem = poMem_DS->CreateLayer("C", poSpatialRef, wkbMultiPolygon, NULL);

	if (poLayer_C == NULL || poLaer_Mem == NULL)
	{
		OGRDataSource::DestroyDataSource(poMem_DS);
		OGRDataSource::DestroyDataSource(poA_DS);
		OGRDataSource::DestroyDataSource(poB_DS);
		OGRDataSource::DestroyDataSource(poC_DS);
		return -1;
	}
	OGRErr err =  poLayer_A->Intersection(poLayer_B, poLaer_Mem, NULL, GDALTermProgress, NULL);
	if (err != OGRERR_NONE)
	{
		OGRDataSource::DestroyDataSource(poMem_DS);
		OGRDataSource::DestroyDataSource(poA_DS);
		OGRDataSource::DestroyDataSource(poB_DS);
		OGRDataSource::DestroyDataSource(poC_DS);
		return -1;
	}

	err = poLayer_B->SymDifference(poLaer_Mem, poLayer_C, NULL, GDALTermProgress, NULL);

	if (err != OGRERR_NONE)
	{
		OGRDataSource::DestroyDataSource(poMem_DS);
		OGRDataSource::DestroyDataSource(poA_DS);
		OGRDataSource::DestroyDataSource(poB_DS);
		OGRDataSource::DestroyDataSource(poC_DS);
		return -1;
	}
	OGRDataSource::DestroyDataSource(poMem_DS);
	OGRDataSource::DestroyDataSource(poA_DS);
	OGRDataSource::DestroyDataSource(poB_DS);
	OGRDataSource::DestroyDataSource(poC_DS);

	return 0;
}

this is my code,is anything wrong?

[Even Rouault]

Ari, this seems to involve the OGR layer algebra API. (I didn't have a look, so it might be just a misuse rather than a bug)

[Ari Jolma]

What exactly seems to be wrong? At least the 22.shp is shifted down due to some projection related issue. I ran a symdifference like this

my $a = Geo::OGR::Open('AW.shp')->Layer;
my $b = Geo::OGR::Open('BW.shp')->Layer;
my $c = Geo::OGR::Driver('ESRI Shapefile')->Create('c.shp');
$a->SymDifference($b, $c->CreateLayer);

and the result is initially shifted similarly as 22.shp but when I set its projection to the same as AW and BW (using QGIS as a viewer) it is correct.

[Even Rouault]

I didn't notice any projection issue. All the .prj are identical. Can this ticket be closed then ?

[Ari Jolma]

The projection issue was something to do with QGIS.

I would like to see the ArcGIS result, which is mentioned in the original posting but otherwise this is works-for-me.