[GIS] gdalbuildvrt error, when using in Python

gdalgeoreferencingpythonvrt

I am new to python/GDAL and am running into perhaps a trivial issue. This may stem from the fact that I don't really understand how to use GDAL properly in python, or something careless, but even though I think I am following the help doc, I keep getting a syntax error when trying to use "gdalbuildvrt".

What I want to do is take several (amount varies for each set, call it N) geotagged 1-band binary rasters [all values are either 0 or 1] of different sizes (each raster in the set overlaps for the most part though), and "stack" them on top of each other so that they are aligned properly according to their coordinate information. I want this "stack" simply so I can sum the values and produce a 'total' tiff that has an extent to match the exclusive extent (meaning not just the overlap region) of all the original rasters. The resulting tiff would have values ranging from 0 to N, to represent the total number of "hits" the pixel in that location received over the course of the N rasters.

I was led to gdalbuildvrt [ http://www.gdal.org/gdalbuildvrt.html ] and after reading about it, it seemed that by using the keyword -separate, I would be able to achieve what I need. However, each time I try to run my program, I get a syntax error. The following shows two of the several different ways I tried calling gdalbuildvrt:

gdalbuildvrt -separate -input_file_list stack.vrt inputlist.txt
gdalbuildvrt -separate stack.vrt inclassfiles

Where inputlist.txt is a text file with a path to the tif on every line, just like the help doc specifies. And inclassfiles is a python list of the pathnames. Every single time, no matter which way I call it, I get a syntax error on the first word after the keywords (i.e. 'inputlist' in inputlist.txt, or 'stack' in stack.vrt).

Could someone please shed some light on what I might be doing wrong? Alternatively, does anyone know how else I could use python to get what I need?

Thanks so much.

Best Answer

First of all, I don't think that gdalbuildvrt will do exactly what you want with the "-separate" option : the stacked file will be a layer of N bands containing each individual image in one of those bands.

Concerning your syntax, I would write :

gdalbuildvrt -separate -input_file_list inputlist.txt stack.vrt

In python, I usually call gdal directly

import subprocess
subprocess.call(["gdalbuildvrt", "-separate", "-input_file_list", "inputlist.txt", "stack.vrt"])

or 

subprocess.call(["gdalbuildvrt", "-separate", "stack.vrt", "im1.tif", "im2.tif" , "im3.tif"])

The extent can be modified with the option "-te xmin ymin xmax ymax" . The values must be expressed in georeferenced units. If not specified, the extent of the VRT is the minimum bounding box of the set of source rasters, so you don't really need to use this option in your case, but I usually do.

If I understand your comment, you should also use -vrtnodata 0 in order to have a value of 0 where there is no image.

Eventually, for the sum of your pixel, it can done within the vrt but this is not straightforward. (see http://www.gdal.org/gdal_vrttut.html) . I prefer using OTB bandmathfilter (see http://orfeo-toolbox.org/otb/otb-applications.html ,it can be wrapped in Python).

Here is an example for building a command list automatically, note that all parameters must be strings:

for year in years:
   command=["gdalbuildvrt",  '-te', '-20015109.354',  '-10007554.677', '20015109.354','10007554.677', path_vrt+ "out_A"+str(year)+".vrt"]       

   list = glob.glob(path_or + "*/MCD64A1.A"+str(year)+ "*.tif")
   for myfile in list:
       command.append(myfile)
   subprocess.call(command)
Related Question