How can I modify a VRTRasterBand sub-class,etc from Python?

Question

I am using GDAL to merge a series of overlapping rasters using Pixel Functions to combine the color bands. As a proof-of-concept, my approach has been to create a VRT, and then to modify it in a text editor:

• Add the VRTDerivedRasterBand subclass to the VRTRasterBand tags
• Add <PixelFunctionCode> (and related) tags to implement the functions.
• Cut&paste the <ComplexSource> tags between bands. Specifically, the new Green (2) band required input from both Red (1) and Green(2) bands. The order here matters, but as long as it is consistent then my custom function can be modified accordingly.

This proof-of-concept had three input rasters, but I wish to work with a variable number.

I also need to automate the entire sequence with one script or program. I'm currently working with Python, but could switch to C/C++.

Can I use the GDAL API to perform these manipulations of the VRT? Can I change the sub-class?

Looks like SetMetadataItem() might be usable to add new XML tags(?) but how can I change the sub-class of an existing VRTRasterBand? Or to remove and replace with a new one?

I can't create a new blank VRT and build it from the ground up, because I do not know the grid size and transform for the final image.

Or is the best approach to directly modify the VRT's XML using something like XSLT?

Answer 1

In the end I directly modified the VRT's XML. Rather than using XSLT, I used Python's lxml library with a sprinkling of XPath. I chose lxml rather than ElementTree because of the need to insert CDATA sections.

The following code adds a PixelFunction to the first band:

from copy import deepcopy
from lxml import etree

vrt_tree = etree.parse(tempfile)
vrt_root = vrt_tree.getroot()
vrtband1 = vrt_root.findall(".//VRTRasterBand[@band='1']")[0]

vrtband1.set("subClass","VRTDerivedRasterBand")
pixelFunctionType = etree.SubElement(vrtband1, 'PixelFunctionType')
pixelFunctionType.text = "find_min"
pixelFunctionLanguage = etree.SubElement(vrtband1, 'PixelFunctionLanguage')
pixelFunctionLanguage.text = "Python"
pixelFunctionCode = etree.SubElement(vrtband1, 'PixelFunctionCode')
pixelFunctionCode.text = etree.CDATA("""
import numpy as np

def find_min(in_ar, out_ar, xoff, yoff, xsize, ysize, raster_xsize, raster_ysize, buf_radius, gt, **kwargs):
    np.amin(in_ar, axis=0, initial=255, out=out_ar)
""")

lxml can also be used to copy the ComplexSource tags from Band 1 to Band 2. The only potential gotcha is that they need to be deep copied:

# Fetch the Red(1) sources - we copy these to Green(2)
vrtband2 = vrt_root.findall(".//VRTRasterBand[@band='2']")[0]
red_sourcelist = vrtband1.findall("ComplexSource")

# Copy the red(1) sources
idx = 0
for src_ele in red_sourcelist:
    vrtband2.insert(idx, deepcopy(src_ele) )
    idx = idx + 1