Platform: Win7 Enterprise; PostgreSQL 9.6.1 with PostGIS 2.3.1; GDAL 2.1.3 and Python 2.7.11.

Like most people who try to do Python things Pythonically, I was pleased when it seemed that there was no more need to run gdalwarp and/or gdal_translate through os.system() or subprocess(). That was short lived once I attempted to fail my way to a working implementation of the Python API calls.

Every now and then I take a 10cm/px aerial of a major city, and snip out an image for every single property boundary so that I can run the images through an image classifier (in order to do some land-use change stuff).

Making a bunch of calls to the command line from within Python always felt kludgy and half-baked, and memory leaks start to add up to a big number.

So I spent some time thrashing through the documentation for key aspects of each GDAL Python API call (for example: creationOptions - is its argument an array? is it a dict of key: value pairs? or is it the worst possible choice - ["option=value","option=value"]? Nobody seemed to want to write down which it was.)

I have started the process with one chunk of the gigantic ECW (the full thing is 131GB, so the task gets done in stages). I have chopped the chunk up into ~3million TIFFs using the BB of the property boundaries. That bit went according to plan; the code (the bit that repeats ~3m times) looks like this

    gdal.Warp(ofn, # the identifier of the parcel, plus '.tif'
          df,      # the first chunk of the behemoth ECW 
          cutlineDSName='PG:' + cs, # cs is the PG connection string
          cutlineSQL=sql_01,        # cut geometry via PG is faster
          cropToCutline=True,  # <- obviously

You might think that setColorInterpretation maybe ought to be true: that's what I thought, but that guaranteed a 100% failure rate.

So far so good.

The usual next step is to run another loop which

  • masks out NODATA using the alpha layer (leaving a 'tilted' image);
  • finds the MBR for the tilted image (which has, as its second element, the dangle-angle);
  • tilts the image 'upright' by the reverse of the dangle-angle (i.e., so that the longest side of the image is one of the two vertical sides); and
  • chucks CUDA at the 3 million 'regularised' images so that a CNN can categorise them faster and more accurately.

As you might glean from the foregoing, I do not need to retain the spatial information about the sub-images: I have all that in the PostGIS data that drives the process that cuts the big ECW into bits.

If I remain within Python, the process stops after the chopping-up-the-raster stage, because something prevents GDAL from recognising the fourth layer.

QGIS can see the 4th layer, and the 4th layer is exactly what I want it to be: 255 inside the image, and 0 outside.

CV2 has the same problem with the image/s as GDAL. I am pretty sure that it means that internally, its code just checks the colorspace string rather than looking for the number of child layers.

This is a process that I've done two dozen times with kludgy os calls, and it consumes significant resources.

I could probably get the next stage done with gdal_translate through subprocess - that's what has worked in the past, and the image has metadata that is not structurally different from past versions of the same task.

The API is clearly more efficient - the gdal.Warp (slice-and-dice) part took less than half the time that the non-API version used to take.

Searches of the web show a low incidence of PER_DATASET ALPHA and gdal ecw "band 1 has no color table".

Has anybody seen this phenomenon - where an RGBA raster is telling software that it's only RGB (despite advertising one of its FOUR layers as 'ALPHA') and the software is believing it?

Here's what gdalinfo says about the imagery (my thoughts are the coloured bits):

enter image description here

I could algebra a way through it - treat it as a numpy array and do some linear algebra - but I am actually really quite keen to improve this workflow so I can pass the task to someone else.

  • 3
    What do you mean by "something prevents GDAL from recognising the fourth layer" and "CV2 has the same problem with the image/s as GDAL"? Are the 4th bands not present? (gdalinfo says they are.) Is your problem with libraries interpreting the images as RGB and not RGBA by default? – Logan Byers Feb 24 '17 at 13:41