Retrieval through SQL from PG DB numerical bug

[thomw00@…]

Retrieval from DB numerical bug:

Step 1: Import test3.shp through:
C:\Program Files\FWTools0.9.9>ogr2ogr -f "PostgreSQL" "PG:user=postgres 
password=thomas! dbname=dcmms host=localhost port=5432" test3.shp

Points in the shapefile have one integer attribute (chloroph) and 2D-
coordinates:
//Point 1: 118; 2.2217631099999999, -1.9194123000000001
//Point 2: 118; 2.1262839100000002, -2.0387613
//Point 3: 118; 2.35702532, -2.3888517199999999

Step 2: retrieve entire dataset:
118.000,2.2217631099999999,-1.9194123000000001
118.000,2.1262839100000002,-2.0387613000000000
118.000,2.3570253200000000,-2.3888517199999999
(correct result)

and through an SQL statement (WHERE chlorph=118; returns all three points in 
this case)
118.000,2.2217631099999999,-1.9194117039444574
118.000,2.1259787342187502,-2.0387613000000000
118.000,2.3570253200000000,-2.3888516454941939
(numerical variations)

Numerical difference is evident in all three points, but most noticeable in x-
coordinate of Point 2.

using: PostgreSQL1.2.1/PostGIS and gdal 1.3.0 through FWTools 0.9.9

Initially these numerical variations were noticed on a buffer(geometry) 
statement using an -sql flag in ogr2ogr, where it produced quite noticeable 
spikes. Possibilities of geos, postgis, postgres causing these numerical 
variations were explored and eliminated through sql retrievals AsText in 
postgres' pgAdminIII (which showed none of these numerical variations). Thus, 
leaving the assumption the numerical variations originate in ogr through the 
translation routines (maybe hex?, as suggested by Frank).

Code used for retrieval (C++, following ogr example on web with addition of 
execute command and output):

int _tmain(int argc, _TCHAR* argv[])
{
	OGRRegisterAll();
	ofstream examplefile ("example.txt"); 

    OGRDataSource       *poDS;

    poDS = OGRSFDriverRegistrar::Open( "PG:user=postgres password=xxx 
dbname=xxx host=localhost port=5432", FALSE );
    if( poDS == NULL )
    {
        printf( "Open failed.\n%s" );
        exit( 1 );
    }


    OGRLayer  *poLayer;

    poLayer = poDS->GetLayerByName( "test3" );

    OGRFeature *poFeature;

    poLayer->ResetReading();
    while( (poFeature = poLayer->GetNextFeature()) != NULL )
    {
        OGRFeatureDefn *poFDefn = poLayer->GetLayerDefn();
        int iField;

		for( iField = 0; iField < poFDefn->GetFieldCount(); iField++ )
        {
            OGRFieldDefn *poFieldDefn = poFDefn->GetFieldDefn( iField );

			
			if( poFieldDefn->GetType() == OFTInteger )
			{
                printf( "%d,", poFeature->GetFieldAsInteger( iField ) );
				examplefile << poFeature->GetFieldAsInteger( 
iField ) << "\n";
			}
            else if( poFieldDefn->GetType() == OFTReal )
			{
                printf( "%.3f,", poFeature->GetFieldAsDouble(iField) );
				examplefile << poFeature->GetFieldAsDouble
(iField) << "\n";
			}
            else if( poFieldDefn->GetType() == OFTString )
			{
                printf( "%s,", poFeature->GetFieldAsString(iField) );
				examplefile << poFeature->GetFieldAsString
(iField) << "\n";
			}
            else
			{
                printf( "%s,", poFeature->GetFieldAsString(iField) );
				examplefile << poFeature->GetFieldAsString
(iField) << "\n";
			}
        }

        OGRGeometry *poGeometry;

        poGeometry = poFeature->GetGeometryRef();


        if( poGeometry != NULL 
            && wkbFlatten(poGeometry->getGeometryType()) == wkbPoint )
        {
            OGRPoint *poPoint = (OGRPoint *) poGeometry;

            printf( "%.16f,%.16f\n", poPoint->getX(), poPoint->getY() );


			examplefile << setprecision (18) << poPoint->getX() 
<< "\n" << poPoint->getY() << "\n"; 
			
        }
        else
        {
            printf( "no point geometry\n" );
        }       


        OGRFeature::DestroyFeature( poFeature );
	}


	examplefile.close();




	OGRLayer *poResultSet;	const char  *pszSQLStatement = NULL;
	pszSQLStatement = "SELECT * FROM test3 WHERE 
chloroph=118"; //pszSQLStatement = "SELECT * FROM test INNER JOIN ksample ON 
test.c_id = ksample.chloroph WHERE something = 467";
        poResultSet = poDS->ExecuteSQL( pszSQLStatement, NULL, NULL );
        if( poResultSet != NULL )
		{
			ofstream resultfile ("results.txt");
			OGRFeature *poFeature2;

			poResultSet->ResetReading();
			while( (poFeature2 = poResultSet->GetNextFeature()) != 
NULL )
			{
				OGRFeatureDefn *poFDefn2 = poResultSet-
>GetLayerDefn();
				int iField;

				for( iField = 0; iField < poFDefn2-
>GetFieldCount(); iField++ )
				{
					OGRFieldDefn *poFieldDefn2 = poFDefn2-
>GetFieldDefn( iField );
				
					if( poFieldDefn2->GetType() == 
OFTInteger )
					{
						printf( "%d,", poFeature2-
>GetFieldAsInteger( iField ) );
						resultfile << poFeature2-
>GetFieldAsInteger( iField ) << "\n"; 
					}
					else if( poFieldDefn2->GetType() == 
OFTReal )
					{
						printf( "%.3f,", poFeature2-
>GetFieldAsDouble(iField) );
						resultfile << poFeature2-
>GetFieldAsDouble(iField) << "\n";
					}
					else if( poFieldDefn2->GetType() == 
OFTString )
					{
						printf( "%s,", poFeature2-
>GetFieldAsString(iField) );
						resultfile << poFeature2-
>GetFieldAsString(iField) << "\n";
					}
					else
					{
						printf( "%s,", poFeature2-
>GetFieldAsString(iField) );
						resultfile << poFeature2-
>GetFieldAsString(iField) << "\n";
					}
				}

				OGRGeometry *poGeometry2;

				poGeometry2 = poFeature2->GetGeometryRef();

				
				if( poGeometry2 != NULL 
					&& wkbFlatten(poGeometry2-
>getGeometryType()) == wkbPoint )
				{
					OGRPoint *poPoint2 = (OGRPoint *) 
poGeometry2;

					printf( "%.16f,%.16f\n", poPoint2->getX
(), poPoint2->getY() );


					resultfile << setprecision (18) << 
poPoint2->getX() << "\n" << poPoint2->getY() << "\n"; 
					
				}
				else
				{
					printf( "no point geometry\n" );
				}       


				OGRFeature::DestroyFeature( poFeature2 );
			}
			resultfile.close();
		}

    poDS->ReleaseResultSet( poResultSet );
    OGRDataSource::DestroyDataSource( poDS );
}

Additional occurance (requires buffer()):
Numerical differences in, for example, buffer:

I noticed numerical differences between an sql statement through the ogr 
library versus directly through postgreSQL (I thought the sqlexecute goes 
through?). Has anyone noticed these differences as well? In the following 
dataset these numerical differences produce sort of ears that are equal to the 
radius of the buffer (so not really negligible and quite visible).

Here is an example (other number combinations are less visible, but still 
produce errors):

Buffer of Point(3.0 -3.7) results:

C:\Program Files\FWTools0.9.9>ogr2ogr -f "ESRI Shapefile" -sql "SELECT buffer(Ge

omFromText('POINT(3.0 -3.7)'), 1.0,4)" buffer.shp "PG:user=postgres 
password=xxx dbname=xxx host=localhost port=5432"


In numbers:

4.0000000000000000,-3.6999999999999957

3.9238604590292496,-4.0820349338195641

3.6289817765446712,-4.2508567811706746

3.3826834323961990,-4.6238413855257363

3.0000000000437552,-4.7000000000000002

2.6173153755938636,-4.6238413855600982

2.2928932188610727,-4.2508567812341678

2.0761204675187801,-4.0820349339023796

2.0000000000000000,-3.7000000000716033

2.0761204674499845,-2.0673165677282341 
 ear 1 
 should be around 2.07.., -
3.06..

2.2928932187340854,-2.9928932188928199

2.6173153754279466,-2.7761192754430666

2.9999999998653095,-2.6999999999999957

3.3826834322302819,-2.7761192753399788

3.6289817810754332,-2.9928932187023385

3.9238604586694903,-2.0673165674793625 
 ear 2 
 should be around 3.92.., -
3.06..

4.0000000000000000,-3.6999999998204118

4.0000000000000000,-3.6999999999999957

Initially, I thought this might be a PostGIS or Geos issue, so I tested the 
same sql statement directly in PostgreSQL using:

SELECT AsText(buffer(GeomFromText('POINT(3.0 -3.7)'), 1.0,4))

"POLYGON((

4 -3.7

3.92387953251558 -4.08268343235472

3.70710678120242 -4.40710678117067

3.3826834323962   -4.6238795324984

3.0000000000449   -4.7

2.61731656768676 -4.62387953253276

2.29289321886107 -4.40710678123417

2.07612046751878 -4.08268343243768

2                              -3.70000000008979

2.07612046745006 -3.31731656772824

2.29289321873409 -2.99289321889282

2.61731656752084 -2.77612046753596

2.99999999986531 -2.7

3.38268343223028 -2.77612046743287

3.70710678107543 -2.99289321870234

3.92387953244686 -3.31731656747936

4                              -3.69999999982041

4                              -3.7

))"


The PostGIS buffer function did not reproduce the ears on the buffer. Also, 
overall there are some differences between the datasets created through ogr2ogr 
versus directly in PostgreSQL. The shapefile is not really an issue since I 
also tested reading the results directly through ogr ...

Thanks,
Thomas

[thomw00@…]

zipped shapefile with 3 points

[warmerdam]

Sorry for dropping the ball on this. 

I have now isolated the error to be failure to translate 'a' or 'A' in
the hex to binary conversion.  I have fixed it in 

gdal/ogr/ogrsf_frmts/pg/ogrpglayer.cpp

and

gdal/port/cpl_string.cpp

Though the one in ogrpglayer.cpp is all that matters for your particular case.

The fix will make it into the GDAL 1.3.1 final release, though it is not
in the GDAL 1.3.1a1 I just published.

Thanks for the detailed bug report!

PS. If you want to patch your local version in the meantime, you can
replace the OGRPGLayer::HEXToGeometry() method in
 gdal/ogr/ogrsf_frmts/pg/ogrpglayer.cpp with the following:

/************************************************************************/
/*                           HEXToGeometry()                            */
/************************************************************************/

OGRGeometry *OGRPGLayer::HEXToGeometry( const char *pszBytea )

{
    GByte       *pabyWKB;
    int iSrc=0, iDst=0;
    OGRGeometry *poGeometry;

    if( pszBytea == NULL )
        return NULL;

    pabyWKB = (GByte *) CPLMalloc(strlen(pszBytea)+1);
    while( pszBytea[iSrc] != '\0' )
    {
        if( pszBytea[iSrc] >= '0' && pszBytea[iSrc] <= '9' )
            pabyWKB[iDst] = pszBytea[iSrc] - '0';
        else if( pszBytea[iSrc] >= 'A' && pszBytea[iSrc] <= 'F' )
            pabyWKB[iDst] = pszBytea[iSrc] - 'A' + 10;
        else if( pszBytea[iSrc] >= 'a' && pszBytea[iSrc] <= 'f' )
            pabyWKB[iDst] = pszBytea[iSrc] - 'a' + 10;
        else 
            pabyWKB[iDst] = 0;

        pabyWKB[iDst] *= 16;

        iSrc++;

        if( pszBytea[iSrc] >= '0' && pszBytea[iSrc] <= '9' )
            pabyWKB[iDst] += pszBytea[iSrc] - '0';
        else if( pszBytea[iSrc] >= 'A' && pszBytea[iSrc] <= 'F' )
            pabyWKB[iDst] += pszBytea[iSrc] - 'A' + 10;
        else if( pszBytea[iSrc] >= 'a' && pszBytea[iSrc] <= 'f' )
            pabyWKB[iDst] += pszBytea[iSrc] - 'a' + 10;
        else 
            pabyWKB[iDst] += 0;

        iSrc++;
        iDst++;
    }

    poGeometry = NULL;
    OGRGeometryFactory::createFromWkb( pabyWKB, NULL, &poGeometry, iDst );

    CPLFree( pabyWKB );
    return poGeometry;
}