While improving process of working with OpenCV in Detecting water areas on satellite images with computer vision I discovered one interesting issue that could be important for many cases of processing. To my surprise this was not mentioned in official docs and even books that described this kind of issues.

One of the steps in this process assumed creating of combined raster from different independent bands that were downloaded from Copernicus hub.

In fact this is a very important step, but there are couple of other related operations that have similar impact on quality of input data. Bands reordering is one them.

Let’s say, we have a multiband raster. In most of the cases this raster will contain at least 3 bands: Red, Green, Blue (and some additional). Normally this order is the default one (RGB). This order and this combination gives us natural colored image where forrest is green, water is blue. At some moment there could be a need of changing the default one. Since working with images in natural color is a rare case in Remote Sensing this is absolutely fine.

How could this issue be addressed?

99% of the people will use gdal_translate utility. 99% of commercial and open-source software will use this utility behind the scenes while performing this issues from GUI.

There is a great tutorial from Robert Simmon how to do this.

Let’s have a look what is the process.

Since we got Sentinel images and converted it for browsing we can run gdalinfo and check what we see for particular 3 bands image.

This is what we have as output

A lot of info, but for now we’re interested in bottom lines. Here we see bands, their names and order (Red is first and Blue is the last one).

Let’s run gdal_translate and see what we’ll get if you change Green and Red bands. In other words, wi’ll change RGB to GRB.

While running we notice first of all performance. On my machine (2.3 GHz Intel Core i5, 8 GB 1333 MHz DDR3) it took about 30 seconds to process 368 MB of input raster and to generate a new one!

Let’s run gdalinfo for newly created raster

GRB. Bands are reordered!

Visually it looks like this

Image processed by gdal_translate

Nice! But performance is poor. In addition, we call this from command line. This is absolutely fine, but what if we’re writing software for this and want to avoid calling this utility from command line? Let’s look if we can find alternative for this in OpenCV! If we can read georeferenced image with OpenCV there has to be a way to reorder bands with it.

I was struggling a lot with OpenCV to get the same image that I got with gdal_translate, but for me it was impossible. I used this snippet:

destinationFile = 'reorderdWithOpenCV.jpeg' inpuDataset = gdal . Open( 'T18TWL_20170423T155131_TCI..gtiff' ) format = 'JPEG' driver = gdal . GetDriverByName(format) outputDataset = driver . CreateCopy(destinationFile, inpuDataset, 0 ) inpuDataset = None outputDataset = None image = cv2 . imread(destinationFile, cv2 . IMREAD_LOAD_GDAL) print image . shape reordered = image[ ... , [ 1 , 0 , 2 ]] cv2 . imwrite(destinationFile, reordered)

I experimented a lot with indexes (direct and reverse order), but I was not able to get the same picture as I got above. What I got looked like this

Image processed by OpenCV

Size of that outcome was also strange – about 27Mb.

Most important that visually this raster is definitely different even from original (non-reordered) image (see below). I suspect that this could be caused by compression/different format of input and output images (stiff and jpeg), but for now this is not so important. This different look and feel means that at least some modifications were made for sure.

Original (non-reordered) image with RGB order of bands

What was even more strange that outcome of gdal_translate showed that bands were not reordered:

Order is RGB.

So, even though OpenCV could correctly change order of bands (internally is uses NumPy arrays) metadata is not updated. And that’s a problem. We can have ‘correct’ look and feel from image visually, but all our modifications should also update image’s metadata. Which OpenCV obviously can’t do.

Flag

cv2 . IMREAD_LOAD_GDAL

supports only ‘reading’, but for writing OpenCV has nothing.

Ok, let’s leave OpenCV for now and check what other people do in this cases.

In his book Joel Lawhead gives this snippet for reordering bands of the raster which for our case would look like this

src = 'T18TWL_20170423T155131_TCI..gtiff' arr = gdal_array . LoadFile(src) output = gdal_array . SaveArray(arr[[ 1 , 0 , 2 ], :], 'reorderedWithArray.gtiff' , format = "GTiff" , prototype = src) output = None

Here gdal_array is used. Good thing that this snippet produces exactly what we want visually:

But what’s with metadata?

We get this

RGB! Though values of the bands are reordered, metadata is not updated.

So, only partially this workaround could be used.

After a long research I found one ready to use solution that could solve both issues of bands reordering: actual reordering and metadata updating.

Here there is a detailed explanation how to implement this in Python and how to use VRT driver.

Adopted version of this snippet looks like this

def read_vsimem (fn): '''Read GDAL vsimem files''' vsifile = gdal . VSIFOpenL(fn, 'r' ) gdal . VSIFSeekL(vsifile, 0 , 2 ) vsileng = gdal . VSIFTellL(vsifile) gdal . VSIFSeekL(vsifile, 0 , 0 ) return gdal . VSIFReadL( 1 , vsileng, vsifile) def write_vsimem (fn,data): '''Write GDAL vsimem files''' vsifile = gdal . VSIFOpenL(fn, 'w' ) size = len (data) gdal . VSIFWriteL(data, 1 , size, vsifile) return gdal . VSIFCloseL(vsifile) def reorderWithPureGDAL (): infn = 'T18TWL_20170423T155131_TCI..gtiff' outfn = 'reorderedWithVRT.gtiff' #List of bands to retain, output bands will be reordered to match bands = [ 2 , 1 , 3 ] ds = gdal . Open(os . path . abspath(infn)) #Ensure path is absolute not relative path vrtdrv = gdal . GetDriverByName( 'VRT' ) tifdrv = gdal . GetDriverByName( 'GTIFF' ) #Create an in memory VRT copy of the input raster vfn = 'test.vrt' vrtds = vrtdrv . CreateCopy(vfn,ds) #Read the XML from memory and parse it #Could also write the VRT to a temp file on disk instead of /vsimem #and used etree.parse(vrtfilepath) instead #i.e. #vfn='some/path/on/disk/blah.vrt' #vrtds=vrtdrv.CreateCopy(vfn,ds) #tree=etree.parse(vfn) #os.unlink(vfn) vrtds = vrtdrv . CreateCopy(vfn,ds) vrtxml = read_vsimem(vfn) tree = etree . fromstring(vrtxml) #Ignore bands not in band list #And put bands we want to keep in the correct order for band in tree . findall( 'VRTRasterBand' ): try : i = bands . index( int (band . attrib[ 'band' ])) band . attrib[ 'band' ] = str (i + 1 ) bands[i] = band except ValueError : pass finally : tree . remove(band) for band in bands: tree . append(band) #Get the XML string from the tree vrtxml = etree . tostring(tree, pretty_print = True ) #Write it to the in memory file #Could also just write to a file on disk write_vsimem(vfn,vrtxml) #Open the VRT containing only the selected bands vrtds = gdal . Open(vfn) #Write to a new raster tifds = tifdrv . CreateCopy(outfn,vrtds, 0 , [ 'COMPRESS=JPEG' , 'PHOTOMETRIC=YCBCR' , 'TFW=YES' ]) #Close dataset to ensure cache is properly flushed del tifds return

Here the trick is in constructing xml tree with metadata and the parsing it. Visually this gives exactly what we want. But what’s with metadata?

GRB! Metadata is updated! So, this approach with VRT could be used for doing this things correctly.

Anyway, bands reordering is not so easy as it might look when we just run gdal_translate. Correct programmatic implementation of this operation includes a lot of additional aspects that we should take under consideration when doing this by ourselves.