I ran into sort of an odd problem with some rasterization code that changed
between 3.10.0 and 3.10.1. The problem also exists with 3.10.2. The core of
the problem is flushing data to disk but it seems that threading and
compression also have something to do with it.

Here's a way to reproduce it... simple made up geometry first, call it
box.geojson:

---
{
  "type": "FeatureCollection",
  "features": [
    {
      "type": "Feature",
      "geometry": {
        "type": "Polygon",
        "coordinates": [
          [
            [ 0.25, 0.25 ],
            [ 0.75, 0.25 ],
            [ 0.75, 0.75 ],
            [ 0.25, 0.75 ],
            [ 0.25, 0.25 ]
          ]
        ]
      }
    }
  ]
}
---

Here's the script:
---
from osgeo import gdal, osr

gdal.UseExceptions()
gdal.SetConfigOption("GDAL_NUM_THREADS", "8")

ds_aoi = gdal.OpenEx("box.geojson")

drv = gdal.GetDriverByName("GTiff")
ds_out = drv.Create("test.tif",
                    xsize=100,
                    ysize=100,
                    bands=1,
                    eType=gdal.GDT_Byte,
                    options=[
                        "COMPRESS=DEFLATE"
                    ])
ds_out.SetGeoTransform((0.0, 0.01, 0.0, 1.0, 0.0, -0.01))
srs = osr.SpatialReference()
srs.ImportFromEPSG(4326)
ds_out.SetSpatialRef(srs)
srs = None

# Uncomment this and it works...
# ds_out.FlushCache()

result = gdal.Rasterize(ds_out,
                        ds_aoi,
                        burnValues=[1],
                        bands=[1],
                        allTouched=True,
                        callback=gdal.TermProgress)
print(f"result = {result}")

ds_aoi = None
ds_out = None
---

Run that script, and it will generate a blank TIF and print "result = 0" if
using GDAL 3.10.1 or 3.10.2. But if using 3.10.0, it prints "result = 1"
and the output looks as expected, a TIF with a box of 1s in the middle of
it and 0s around the outside.

Side note... it's also weird that the gdal.Rasterize doesn't seem to throw
any exception, but I'm not even sure what the problem is. I tried
CPL_DEBUG=ON and don't see anything obvious:
---
GeoJSON: First pass: 100.00 %
GDAL: GDALOpen(box.geojson, this=0x356de0d0) succeeds as GeoJSON.
GDAL: GDALDriver::Create(GTiff,test.tif,100,100,1,Byte,0x357529a0)
GTiff: Using up to 8 threads for compression/decompression
GDAL: GDAL_CACHEMAX = 799 MB
GDAL: Rasterizer operating on 1 swaths of 100 scanlines.
0result = 0
GDAL: GDALClose(box.geojson, this=0x356de0d0)
GTIFF: Waiting for worker job to finish handling block 0
GDAL: GDALClose(test.tif, this=0x35ab1d80)
---

It seems that this combination is required to produce the error:
- Set GDAL_NUM_THREADS
- Specify COMPRESS in the GTiff creation options
- Do NOT flush the cache or close the file between creation and
rasterization

If any of those conditions are changed, it appears to work ok.

Is this expected? Is that FlushCache call supposed to be necessary? It's
odd that it wasn't required before, but it is now.
_______________________________________________
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Reply via email to