Although it would require another library (which currently is only available on OS X and Linux) you could use RSGISLib (http://rsgislib.org), which is built on top of GDAL to do this. There is a function to stack bands, as an example:
#!/usr/bin/env python
import rsgislib
from rsgislib import imageutils
# Create list of images
imageList = ['LC80430342013291LGN00_B1.TIF',
'LC80430342013291LGN00_B2.TIF',
'LC80430342013291LGN00_B3.TIF',
'LC80430342013291LGN00_B4.TIF',
'LC80430342013291LGN00_B5.TIF',
'LC80430342013291LGN00_B6.TIF',
'LC80430342013291LGN00_B7.TIF']
# Create list of band names
bandNamesList = ['Coastal','Blue','Green','Red','NIR','SWIR1','SWIR2']
# Set output image
outputImage = 'LC80430342013291LGN00_stack.tif'
# Set format and type
gdalFormat = 'GTiff'
dataType = rsgislib.TYPE_16UINT
# Stack
imageutils.stackImageBands(imageList, bandNamesList, outputImage, None, 0, gdalFormat, dataType)
It also allows passing in a list of band names.
However, if you were going to install RSGISLib to stack the Landsat bands you could just install ARCSI, which would stack the bands and convert to radiance or reflectance (if you were looking to do this). More details are here: http://spectraldifferences.wordpress.com/2014/05/27/arcsi/
Unfortunately, your problem description leaves some desires. You really should specify more clearly what does not work and which error messages you receive and so forth. And usually it is best in case of errors to try one step after the other.
In your case I see some problems that may cause your script to "not work". First, you should specify the correct path to your programs and to your data. Working with relative paths is problematic, since your script may not work under the same shell environment as you do. So, your relative path newfolder/band4.tif
may be the reason why the command does not work.
Second, you are fiddling too much with python commands you obviously do not understand well enough. If you want to continue to use the subprocess.Popen
it's best to read the manual page and understand it.
Anyway, a few tipps now that should get you running.
- use the full path to gdalwarp, ie. /usr/local/bin/gdalwapr
- do not use subprocess.Popen, but use subprocess.check_call() instead 1)
- omit all those proc.communicate() stuff, you probably don't need it
- '+proj=longlat +ellps=WGS84' is equal to 'EPSG:4326'
1) subprocess.check_call()
waits for the completion of the command before it continues in your script. Nice for seeing where your problem lies and vital if you want to treat the same data sequentially. You may substitute subprocess.check_call()
with other subprocess commands after your program runs
A little example script:
def warp(args):
"""with a def you can easily change your subprocess call"""
# command construction with binary and options
options = ['/usr/bin/gdalwarp']
options.extend(args)
# call gdalwarp
subprocess.check_call(options)
# just pass the list of commands without the gdalwarp command at the beginning, since it is added automatically
# this should run two commands of gdalwarp sequentially
warp(['-t_srs', 'EPSG:4326', '/home/user/maps/source.tif', '/home/user/maps/target.tif'])
warp(['-t_srs', 'EPSG:3857', '-of', 'AAIGrid', '/home/user/maps/source.tif', '/home/user/maps/target.xyz']))
I have a little script which you can examine. It creates a combined hillshade from SRTM data. It does not run well at the moment, due to other reasons though. There you can have a look how I've organized my workflow for lots of GDAL commands, and Imagemagick too.
Best Answer
Using
rasterio
you could doIt of course assumes that your input layers already share the same extent, resolution and data type