[gdal-dev] Hidden regression in GDAL 1.11.0 - gdal_rasterize with postgis silently fails to render many polygons.

[Graeme B. Bell]

TLDR;

gdal_rasterize has begun silently producing wrong output at least since GDAL 1.11.0. 
Whenever  >500 polygons are present in the database query result (and perhaps under other circumstances), only 500 will be drawn. 


Hi,

I am rasterizing many maps and areas of maps in SRID 25833 using gdal_rasterize with a postgis database. 

When I use GDAL 1.11, the output is missing much of the data, despite gdal_rasterize reporting completing successfully. 
With GDAL 1.09 and GDAL 1.10, the output is correct and complete and identical.

I'm posting here to gather more information about the bug and also to raise an alarm for anyone using postgis/gdal_rasterize (1.11) in their data workflow over the last year. 

Background of bug
--------

Versions:  

POSTGIS="2.1.5 r13152" GEOS="3.4.2-CAPI-1.8.2 r3921" PROJ="Rel. 4.8.0, 6 March 
2012" GDAL="GDAL 1.11.0, released 2014/04/16"    (from EPEL RPM)

A. The source table in postgis is registered in geometry_columns correctly. 

B. I run this command in Centos linux. 

gdal_rasterize -te 1060000 7865000 1120000 7940000 -a jord -init 255 -tr 100 100 -a_nodata 255 -a_srs 'EPSG:25833' -co 'TILED=YES' -co 'BIGTIFF=IF_SAFER' -co 'COMPRESS=PACKBITS' -co 'PREDICTOR=1' -co 'ZLEVEL=1' -ot Byte PG:"host=XXXX'' dbname='XXXX' user='XXXX' password='XXXX' tables='XXXX'" -sql "SELECT xxxx FROM xxxx WHERE st_intersects(geo, ST_GeometryFromText('POLYGON((1060000 7865000,1060000 7940000,1120000 7940000, 1120000 7865000, 1060000 7865000))',25833)) " xxxx.tif

C. The SQL query (-sql), when run directly on postgres with COUNT(), produces this output. 

 count 
-------
  3959
(1 row)

D. GDAL 1.9.2 and GDAL 1.10.0 produce a TIF file of 12238 bytes which is visually complete. 

E. GDAL 1.11.0 produces a TIF file of 10920 bytes which is missing many polygons.

F. Console output from gdal_rasterize in all 3 cases is as follows;
0...10...20...30...40...50...60...70...80...90...100 - done.

G. PG server log output for GDAL 1.9.2 and GDAL 1.10.0

(nothing...)

PG server log output for GDAL 1.11.0

< 2015-02-09 13:52:11.720 CET >STATEMENT:  FETCH 500 in OGRPGLayerReader0xd5ce50
< 2015-02-09 13:52:11.720 CET >ERROR:  cursor "ogrpglayerreader0xd5ce50" does not exist
< 2015-02-09 13:52:11.720 CET >STATEMENT:  CLOSE OGRPGLayerReader0xd5ce50
< 2015-02-09 13:52:18.657 CET >ERROR:  cursor "ogrpglayerreader0x1fd8e50" does not exist

H. Debug mode output for GDAL 1.11.0

OGR: OGROpen(...) succeeded as PostgreSQL.
GDAL: QuietDelete(xxxx.tif) invoking Delete()
GDAL: GDALOpen(xxxx.tif, this=0x259a8f0) succeeds as GTiff.
GDAL: GDALDefaultOverviews::OverviewScan()
GDAL: GDALClose(xxxx.tif, this=0x259a8f0)
GDAL: GDALDriver::Create(GTiff,xxxx.tif,600,750,1,Byte,0x25744e0)
PG: Primary key name (FID): gid
PG: Using column 'gid' as FID for table 'xxxx_25833'
GDAL: Rasterizer operating on 1 swaths of 750 scanlines.
0...10...20...30...40...50...60...70...80...90...100 - done.
PG: 500 features read on layer 'sql_statement'.
GDAL: GDALClose(xxxx.tif, this=0x2598930)

I. Debug mode output for GDAL 1.10.0

OGR: OGROpen(xxxx) succeeded as PostgreSQL.
GDAL: QuietDelete(xxxx) invoking Delete()
GDAL: GDALOpen(xxxx) succeeds as GTiff.
GDAL: GDALDefaultOverviews::OverviewScan()
GDAL: GDALClose(xxxx.tif, this=0xeb5030)
GDAL: GDALDriver::Create(GTiff,xxxx.tif,600,750,1,Byte,0xe8ede0)
PG: Primary key name (FID): gid
PG: Using column 'gid' as FID for table 'dmk_25833'
PG: Layer 'xxxx_25833' geometry type: GEOMETRY:Unknown (any), Dim=2
0...10...20...30...40...50...60...70...80...90...100 - done.
PG: 3959 features read on layer 'sql_statement'.
GDAL: GDALClose(xxxx.tif, this=0xeb2fd0)


As well as local unix sockets I've tested TCP connections to a remote PG database from GDAL1.10 and GDAL1.11 clients, and I've tested with other data layers. Every time, I'm getting the same problem e.g. 500 features out 8000.

Has anyone else experienced anything similar? 

I've spent some time looking through the source but I wasn't able to determine exactly where the problem is occuring. 
It seems to somehow be related to the select page size (500 rows) being used in place of the total feature count; or perhaps some logic mistake causing the reads to finish after 1 page is complete.

Graeme Bell