2

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.

6

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)
  • Okay, that makes more sense! I actually do want an N-band file like you described though. My only concern is having the pixels aligned correctly according to their coordinate information, because the rasters vary in size. I just want the smaller rasters to be padded with 0's so that each layer in the "stack" is the same size, and aligned. Do you know if using gdalbuildvrt (or anything) can/will reproject the layers to be of the same size (such that the extent is the union of all N images)? I really need to figure out a way, so if you could give me any suggestions, I'd be so grateful. Thanks! – user20408 Jan 3 '14 at 2:09
  • hmm, okay. that sounds like it can do what I need it to. Is there a way I can call the gdalbuildvrt process, but while using a variable as one of the inputs? I am going to be creating several hundred stacks using a loop, so the inputlist is going to change after every iteration. How can I define a variable such that gdalbuildvrt as a subprocess will understand? Or can I just simply define it in my python script and insert it as the input in subprocess.call([... ? Sorry for all the questions. Thanks so much. – user20408 Jan 3 '14 at 16:06
  • of course, the command you launch with subprocess can be handled like any list in Python (see update). – radouxju Jan 3 '14 at 19:03
  • that's super helpful. i was able to get a good result from gdalbuildvrt using your help. could you give me a hint on how to go about doing the sum w/i the vrt? I couldn't find anything helpful (that I can tell at least) on the link. I thought I could readAsArray the N bands of the vrt & insert each into the 1st dim of a 3d numpy array (of shape (Nbands, Nsamples, Nlines)), but it kept giving me "Error: Array is too big." After defining the dtype = byte, I now recieve MemoryError. So if I can figure out how to just do the sum in the vrt (or get around memory issues), i could be done. thank you! – user20408 Jan 3 '14 at 22:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.