[GIS] Using GDAL/Python to stack georeferenced images of different sizes

gdalgeoreferencingpythonstack

I am in the process of porting a code I wrote in IDL (interactive data language) to Python but am running into a bit of a problem that I am hoping someone can help me with.

The code goes like this:

take individual classified Landsat GeoTIFFs (say there are N individual 1-band files per scene, each representing a different day) and further reduce these images to three binary-themed 1-band images (water and not water, land and not land, water/land and not water/land). This will be done by reading the rasters as matrices and replacing values.
** I don't actually need to have these images, so I can save them as memory or just keep them as numpy ndarrays to move to the next step
stack these images/arrays to produce 3 different (1 for each 'element') N-band stacks (or a 3-dimensional array– (samples, lines, N)) for each scene
total the stacks to get a total number of water/land/water&land observations per pixel (produces one 1-band total image for each scene)
other stuff

The problem I am running into is when I get to the stacking, as the individual images for each scene vary in size, although they mostly overlap with each other. I originally used an ENVI layer-stacking routine that takes the N different-sized 1-band images for each scene and stacks them into an N-band image with an extent that encompasses all of the images' extents, and then reading the resulting rasters in as 3-D arrays to do the totals. I would like to do something similar with GDAL/Python but am not sure how to go about doing so. I was thinking I would implement GDAL capabilities of GeoTIFFs by using the geotransform info of the images to somehow find the inclusive extent, possibly padding the edges of the images with 0's so they are all the same size, stacking these images/3-d arrays so that they are correctly aligned, then computing the totals. Hopefully there is an easier way, as I'm not sure how to pull that off.

Does anyone have any suggestions or ideas as to what would be the most efficient way (or, any way really), to do what I need to do? I'm open to anything.

Best Answer

I'm not totally certain on what you need, but I think:

gdalbuildvrt -separate stack.vrt lsat1.tif lsat2.tif ...

Should give you a gdal dataset that is 'stacked and is covers the extent of all images. If you need a tif after that, use

gdal_translate stack.vrt stack.tif

Related Solutions

[GIS] raster algebra in python with rasters of different extents

I think you can do that pretty easily with GDAL and numpy. Mind you, I think that you will need to do more complex analysis afterwards (speckle reduction etc), but in principle, you could stack all your observations into a single multiband file using eg

gdalbuildvrt -separate -te xmin ymin xmax ymax -input_file_list my_filenames.txt output_file.vrt

output_file.vrt is a dataset that gdal should understand. xmin ymin xmax ymax are the maximum extent (and I'm assuing they all your observations share the same spatial resolution). In python, you load up the data, set a mask for your no data value (I'm assuming 0 here, but do check), and then average:

from osgeo import gdal
g = gdal.Open ( "output_file.vrt" )
data = g.ReadAsArray()
mdata = np.ma.array ( data, mask=( data == 0 ) )
mean_s0 = mdata.mean ( axis=0 )

[GIS] gdalbuildvrt error, when using in Python

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"])

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)

Best Answer

Related Solutions

[GIS] raster algebra in python with rasters of different extents

[GIS] gdalbuildvrt error, when using in Python

Related Question