Opened 9 years ago

Closed 5 years ago

#5836 closed defect (wontfix)

gdal_rasterize cannot draw straight polygons

Reported by: Even Rouault Owned by: warmerdam
Priority: normal Milestone: closed_because_of_github_migration
Component: Algorithms Version: unspecified
Severity: normal Keywords:
Cc:

Description

Let's create a layer with

from osgeo import ogr
from osgeo import osr

ds = ogr.Open('pg:dbname=autotest')
srs = osr.SpatialReference()
srs.ImportFromEPSG(4326)
lyr = ds.CreateLayer('world_rectangles', geom_type = ogr.wkbPolygon, srs = srs, options = ['OVERWRITE=YES'])
lyr.StartTransaction()
for y in range(180):
    for x in range(360):
        f = ogr.Feature(lyr.GetLayerDefn())
        x0 = x - 180
        y0 = 90 - y - 1
        x1 = x0 + 1
        y1 = y0 + 1
        f.SetGeometry(ogr.CreateGeometryFromWkt('POLYGON((%d %d,%d %d,%d %d,%d %d,%d %d))' % (x0,y0,x1,y0,x1,y1,x1,y0,x0,y0)))
        lyr.CreateFeature(f)
lyr.CommitTransaction()

and then rasterize it with :

gdal_rasterize -sql "select * from world_rectangles where st_intersects(wkb_geometry, ST_GeometryFromText('POLYGON((-180 -90,-180 90,180 90,180 -90,-180 -90))',4326))"  "pg:dbname=autotest"  out.tif -burn 255 -ts 4000 2000  -a_srs EPSG:4326 -te -180 -90 180 90 --debug on

Nothing is burnt ! But if using diamond shaped polygons instead, it works

        x0 = x - 180 + 0.5
        y0 = 90 - y - 1
        f.SetGeometry(ogr.CreateGeometryFromWkt('POLYGON((%f %f,%f %f,%f %f,%f %f,%f %f))' % (x0,y0,x0-0.5,y0+0.5,x0,y0+1,x0+0.5,y0+0.5,x0,y0)))

Change History (2)

comment:1 by meier, 8 years ago

The issue is with the gdal_rasterize utility, as the correct(?) shapefiles are written in the database using the Python script both times. On trying the gdal_rasterize utility with the ALL_TOUCHED (-at) option, the correct result was generated as opposed to the pixels being filled with NaN values. Add the -at option to make it work. Checked against

$ gdal-config --version
2.0.1

The following works after running both the python scripts -

$ gdal_rasterize -sql "select * from world_rectangles where st_intersects(wkb_geometry, ST_GeometryFromText('POLYGON((-180 -90,-180 90,180 90,180 -90,-180 -90))',4326))" "pg:dbname=TestGIS host=localhost user=***** password=****** port=5433"  out.tif -burn 100 -ts 4000 2000 -at -a_srs EPSG:4326 -te -180 -90 180 90 --debug on
PG: DBName="TestGIS"
PG: PostgreSQL version string : 'PostgreSQL 9.4.5 on i686-pc-linux-gnu, compiled by gcc (Ubuntu 4.8.2-19ubuntu1) 4.8.2, 32-bit'
PG: PostGIS version string : '2.1 USE_GEOS=1 USE_PROJ=1 USE_STATS=1'
GDAL: GDALOpen(pg:dbname=TestGIS host=localhost user=postgis password=XXXXXXX port=5433, this=0x8392250) succeeds as PostgreSQL.
GDAL: GDALDriver::Create(GTiff,out.tif,4000,2000,1,Float64,(nil))
PG: Primary key name (FID): ogc_fid, type : int4
PG: Using column 'ogc_fid' as FID for table 'world_rectangles'
GDAL: Rasterizer operating on 7 swaths of 312 scanlines.
0...10...20...30...40...50...60...70...80...90...100 - done.
PG: 64800 features read on layer 'sql_statement'.
GDAL: GDALClose(pg:dbname=TestGIS host=localhost user=postgis password=XXXXXXX port=5433, this=0x8392250)
GDAL: GDALClose(out.tif, this=0x83b9620)

comment:2 by Even Rouault, 5 years ago

Milestone: closed_because_of_github_migration
Resolution: wontfix
Status: newclosed

This ticket has been automatically closed because Trac is no longer used for GDAL bug tracking, since the project has migrated to GitHub. If you believe this ticket is still valid, you may file it to https://github.com/OSGeo/gdal/issues if it is not already reported there.

Note: See TracTickets for help on using tickets.