2

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:

  1. 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

  2. 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

  3. 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)
  4. 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.

Improve this question
0

1 Answer 1

Reset to default
3

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
Improve this answer
8
  • Thanks @kyle! I'll take a look at this. I don't necessarily need a tif of the stack, but I need to be able to translate the stack into a 3-d numpy array (specifically of shape #images x #rows x #columns) so I can then add together the pixel values of each image (producing a 2-d array of shape #rows x #cols) and convert that into a 1-band tiff. Any ideas on how to do that? I'm not very familiar with gdal, but would I be able to convert the stack.vrt into a 3-d numpy array like specified above? Thanks!!!
    – user20408
    Commented Dec 30, 2013 at 22:48
  • 1
    Winging it a little here, but if you just need to do addition, try gdal_calc.py at: trac.osgeo.org/gdal/browser/trunk/gdal/swig/python/scripts/… It should be on the utilities page, as I think it's distributed, but it's not right now. Otherwise, just open the vrt, read each band using ReadAsArray(), and I bet numpy's '+' operator will to element-wise addition. So: ds = gdal.Open('my.vrt') b1 = ds.GetRasterBand(1).ReadAsArray() b2 = ds.GetRasterBand(2).ReadAsArray() b3 = ds.GetRasterBand(3).ReadAsArray() all = b1+b2+b3
    – user10353
    Commented Dec 30, 2013 at 22:53
  • awesome! Do you know if gdal_calc takes into account the fact that the images are different sizes? i'm just worried about summing the images while they are not aligned according to their coordinate information. If not, the vrt method makes sense to me, if I can implement it correctly. Another concern of mine is making sure the total tif has the correct metadata/geotag information. Sorry for all the questions, I'm new to using python/gdal for image manipulation (and fairly new to python for that matter). Thanks much!
    – user20408
    Commented Dec 30, 2013 at 23:07
  • No clue. The VRT should work, although you may have to use a numpy masked array to handle no data issues. I've not had to use that much.
    – user10353
    Commented Dec 30, 2013 at 23:26
  • do you know what i'm doing wrong here? i keep getting a syntax error (tried several different arrangements including using variables, pathnames, using a list in a text file, etc.): <br/> gdalbuildvrt -separate -input_file_list inputlist.txt stack.vrt <br/> #gdalbuildvrt -separate stack.vrt inclassfiles <br/> And other ways. What am I missing? Thanks
    – user20408
    Commented Dec 31, 2013 at 20:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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