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
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
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 format='GTiff', cutlineDSName='PG:' + cs, # cs is the PG connection string cutlineSQL=sql_01, # cut geometry via PG is faster cropToCutline=True, # <- obviously copyMetadata=True, metadataConflictValue=None, srcAlpha=True, dstAlpha=True, creationOptions=["PHOTOMETRIC=RGB"], setColorInterpretation=False, callback=None, callback_data=None)
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
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
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?
gdalinfo says about the imagery (my thoughts are the coloured bits):
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.