[gdal-dev] VRT derived band pixel functions written in Python

Even Rouault even.rouault at spatialys.com
Mon Sep 12 05:31:43 PDT 2016


I wanted to mention a new (and I think pretty cool ) feature I've added in 
trunk: the possibility to define pixel functions in Python in a VRT derived 

Let's start with a simple example to multiply a raster by 1.5 :

<VRTDataset rasterXSize="20" rasterYSize="20">
  <VRTRasterBand dataType="Byte" band="1" subClass="VRTDerivedRasterBand">
    <PixelFunctionArguments factor="1.5"/>
import numpy as np
def multiply(in_ar, out_ar, xoff, yoff, xsize, ysize, raster_xsize,
                   raster_ysize, buf_radius, gt, **kwargs):
    factor = float(kwargs['factor'])
    out_ar[:] = np.round_(np.clip(in_ar[0] * factor,0,255))
      <SourceFilename relativeToVRT="1">byte.tif</SourceFilename>

or to add 2 (or more) rasters:

<VRTDataset rasterXSize="20" rasterYSize="20">
  <VRTRasterBand dataType="Byte" band="1" subClass="VRTDerivedRasterBand">
import numpy as np
def add(in_ar, out_ar, xoff, yoff, xsize, ysize, raster_xsize,
                   raster_ysize, buf_radius, gt, **kwargs):
    np.round_(np.clip(np.sum(in_ar, axis = 0, dtype = 'uint16'),0,255),
              out = out_ar)
      <SourceFilename relativeToVRT="1">byte.tif</SourceFilename>
      <SourceFilename relativeToVRT="1">byte2.tif</SourceFilename>

You can put any python module code inside PixelFunctionCode, with at least one 
function with the following arguments :
- in_ar: list of input numpy arrays (one numpy array for each source)
- out_ar: output numpy array to fill (initialized at the right dimensions and 
with the VRTRasterBand.dataType)
- xoff: pixel offset to the top left corner of the accessed region of the band
- yoff: line offset to the top left corner of the accessed region of the band
- xsize: width of the region of the accessed region of the band
- ysize: height of the region of the accessed region of the band
- raster_xsize: total with of the raster band
- raster_ysize: total height of the raster band
- buf_radius: radius of the buffer (in pixels) added to the left, right, top 
and bottom of in_ar / out_ar
- gt: geotransform
- kwargs: dictionnary with user arguments defined in PixelFunctionArguments

For basic operations, you just need to care about in_ar and out_ar.

With all that, you can do interesting stuff like hillshading (code ported from 

<VRTDataset rasterXSize="121" rasterYSize="121">
  <VRTRasterBand dataType="Byte" band="1" subClass="VRTDerivedRasterBand">
      <SourceFilename relativeToVRT="1">n43.dt0</SourceFilename>
    <PixelFunctionArguments scale="111120" z_factor="30" />
# Licence: X/MIT
# Copyright 2016, Even Rouault
import math

def hillshade_int(in_ar, out_ar, xoff, yoff, xsize, ysize, raster_xsize,
                         raster_ysize, radius, gt, z, scale):
    ovr_scale_x = float(out_ar.shape[1] - 2 * radius) / xsize
    ovr_scale_y = float(out_ar.shape[0] - 2 * radius) / ysize
    ewres = gt[1] / ovr_scale_x
    nsres = gt[5] / ovr_scale_y
    inv_nsres = 1.0 / nsres
    inv_ewres = 1.0 / ewres

    az = 315
    alt = 45
    degreesToRadians = math.pi / 180

    sin_alt = math.sin(alt * degreesToRadians)
    azRadians = az * degreesToRadians
    z_scale_factor = z / (8 * scale)
    cos_alt_mul_z_scale_factor = \
              math.cos(alt * degreesToRadians) * z_scale_factor
    cos_az_mul_cos_alt_mul_z_scale_factor_mul_254 = \
                254 * math.cos(azRadians) * cos_alt_mul_z_scale_factor
    sin_az_mul_cos_alt_mul_z_scale_factor_mul_254 = \
                254 * math.sin(azRadians) * cos_alt_mul_z_scale_factor
    square_z_scale_factor = z_scale_factor * z_scale_factor
    sin_alt_mul_254 = 254.0 * sin_alt

    for j in range(radius, out_ar.shape[0]-radius):
        win_line = in_ar[0][j-radius:j+radius+1,:]
        for i in range(radius, out_ar.shape[1]-radius):
            win = win_line[:,i-radius:i+radius+1].tolist()
            x = inv_ewres * ((win[0][0] + win[1][0] + win[1][0] + win[2][0])-\
                             (win[0][2] + win[1][2] + win[1][2] + win[2][2])) 
            y = inv_nsres * ((win[2][0] + win[2][1] + win[2][1] + win[2][2])-\
                             (win[0][0] + win[0][1] + win[0][1] + win[0][2]))
            xx_plus_yy = x * x + y * y
            cang_mul_254 = (sin_alt_mul_254 - \
                (y * cos_az_mul_cos_alt_mul_z_scale_factor_mul_254 - \
                    x * sin_az_mul_cos_alt_mul_z_scale_factor_mul_254)) / \
                math.sqrt(1 + square_z_scale_factor * xx_plus_yy)
            if cang_mul_254 < 0:
                out_ar[j,i] = 1
                out_ar[j,i] = 1 + round(cang_mul_254)

def hillshade(in_ar, out_ar, xoff, yoff, xsize, ysize, raster_xsize, 
              raster_ysize, radius, gt, **kwargs):
    z = float(kwargs['z_factor'])
    scale= float(kwargs['scale'])
    hillshade_int(in_ar, out_ar, xoff, yoff, xsize, ysize, raster_xsize,
                  raster_ysize, radius, gt, z, scale)


You can completely offload the python code itself into a proper my_lib.py file 
and just specify

Technically, the interfacing with Python is done at run-time by dynamically 
loading Python symbols with dlopen()/GetProcAddress(), when they are already 
available in the process. For example if libgdal is loaded from a Python 
interpreter, or from a binary like QGIS which has already loaded the Python 
lib. Otherwise a few shared objects ("libpython2.7.so", "python27.dll",  
"libpython3.4m.so", etc.) are tried, unless the PYTHONSO config option is 
specified to point to a precise filename. The advantage of this approach is that 
the same GDAL library binary is compatible of all Python 2.X (tested: 2.6, 
2.7) and 3.X (tested: 3.1, 3.4) versions, and any numpy version that comes in 
the Python environment used (at compilation time you don't need any 
python/numpy development package). The numpy dependency is not a critical one: 
one could imagine a fallback mode where Python arrays would be used instead, 
but this has likely little interest.

Successfully tested on Linux, MacOSX, FreeBSD and Windows. Some extra tweaking 
of the predefined set of shared object names - that are probed when no already 
loaded Python environment is found - might be needed.
The implementation should be thread-safe regarding use of the Python Global 
Interpreter Lock (GIL).

I've also tested with the numba Just-In-Time compiler 
(http://numba.pydata.org/) and it provides major performance improvements for 
highly computational code (example given below).
With numba enabled, I found that my above Python hillshading on a 10Kx10K 
float32 raster was even faster than gdaldem (the reason is that the Python 
version is a bit simplified as it doesn't take into account input nodata 
values), whereas the non-jit'ed one is 100x slower. When removing the nodata 
flag, gdaldem is only twice faster as the jit'ed python code. So this is a good 
sign that such approach isn't only a toy or just for prototyping.
Speaking of JIT, there's no provision (yet?) for interfacing with PyPy. Would 
require a new backend as PyPy C API has nothing to do with the CPython one.

There are obvious security concerns in allowing Python code to be run when 
getting the content of a vrt file. The GDAL_VRT_ENABLE_PYTHON config option = 
IF_SAFE / NO / YES can be set to control the behaviour. The default is IF_SAFE 
(can be change at compilation time by defining 
-DGDAL_VRT_ENABLE_PYTHON_DEFAULT="NO" e.g. And Python code execution can be 
completely disabled with -DGDAL_VRT_DISABLE_PYTHON). Safe must be understood 
as: the code will not read, write, delete... files, spawn external code, do 
network activity, etc. Said otherwise, the code will only do "maths". But 
infinite looping is something definitely possible in the safe mode. The 
heuristics of the IF_SAFE mode is rather basic and I'd be grateful if people 
could point ways of breaking it. If any of the following strings - "import" 
(unless it is "import numpy" / "from numpy import ...", "import math" / "from 
math import ..." or "from numba import jit"), "eval", "compile", "open", 
"load", "file", "input", "save", "memmap", "DataSource", "genfromtxt", 
"getattr", "ctypeslib", "testing", "dump", "fromregex" - is found anywhere in 
the code, then the code is considered unsafe (there are interestingly a lot of 
methods in numpy to do file I/O. Hopefully I've captured them with the previous 
filters). Another 'unsafe' pattern is when the pixel function references an 
external module like my above my_lib.hillshade example (who knows if there 
will not be some day a shutil.reformat_your_hard_drive function with the right 

This new capability isn't yet documented in the VRT doc, although this message 
will be a start.

I'm interested in feedback you may have.

And to conclude with a fun example: a raster with a Mandelbrot fractal. Just a 
grey-level version. Let to the reader as an exercice: add a color table. To be 
opened for example in QGIS and enjoy the almost infinite zoom feeling. Make 
sure to disable contrast enhancement. It uses numba when available, and when 
this is the case, it's really fast when paning/zooming.

<VRTDataset rasterXSize="100000000" rasterYSize="100000000">
  <VRTRasterBand dataType="Byte" band="1" subClass="VRTDerivedRasterBand">
    from numba import jit
    #print('Using numba')
    g_max_iterations = 100
    class jit(object):
        def __init__(self, nopython = True, nogil = True):

        def __call__(self, f):
            return f

    #print('Using non-JIT version')
    g_max_iterations = 25

# Use a wrapper since the VRT code cannot access the jit decorated function
def mandelbrot(in_ar, out_ar, xoff, yoff, xsize, ysize, raster_xsize,
                        raster_ysize, r, gt, **kwargs):
    mandelbrot_jit(out_ar, xoff, yoff, xsize, ysize, raster_xsize, raster_ysize, 

@jit(nopython=True, nogil=True)
def mandelbrot_jit(out_ar, xoff, yoff, xsize, ysize, raster_xsize,
                        raster_ysize, max_iterations):
    ovr_factor_y = float(out_ar.shape[0]) / ysize
    ovr_factor_x = float(out_ar.shape[1]) / xsize
    for j in range( out_ar.shape[0]):
        y0 = 2.0 * (yoff + j / ovr_factor_y) / raster_ysize - 1
        for i in range(out_ar.shape[1]):
            x0 = 3.5 * (xoff + i / ovr_factor_x) / raster_xsize - 2.5
            x = 0.0
            y = 0.0
            x2 = 0.0
            y2 = 0.0
            iteration = 0
            while x2 + y2 < 4 and iteration < max_iterations:
                y = 2*x*y + y0
                x = x2 - y2 + x0
                x2 = x * x
                y2 = y * y
                iteration += 1

            out_ar[j][i] = iteration * 255 / max_iterations
      <MDI key="STATISTICS_MEAN">127</MDI>
      <MDI key="STATISTICS_STDDEV">127</MDI>

Statistics have been added just to make QGIS open the file a bit quicker.


Spatialys - Geospatial professional services

More information about the gdal-dev mailing list