Given a time series of an area (single band), how you subset them and combine the subsets into a single multiband file (eg. vrt) using gdal or gdals' bindings for python?
python – Creating Subsetted Multiband Image from Multiple Files Using GDAL or Python
boundless-suitegdalpython
Related Solutions
Update:
As of GDAL 2.1, you can use the gdal.Translate() function and pass the bands in the bandList argument (list of arguments).
from osgeo import gdal
infn='in.tif'
outfn='out.tif'
gdal.Translate(outfn, infn, bandList=[1,2,3], **other_kwargs)
Original answer:
gdal_translate does this by building a VRT from scratch, copying all the metadata from the source dataset (unless everything matches, then it just uses CreateCopy).
You can do this in python but it's a bit fiddly. An easier way is use the VRT driver to create a VRT copy, then modify the VRT XML to remove the bands you aren't interested in before writing the modified VRT out to new file. See example below:
from osgeo import gdal
from lxml import etree
import os,sys
def read_vsimem(fn):
'''Read GDAL vsimem files'''
vsifile = gdal.VSIFOpenL(fn,'r')
gdal.VSIFSeekL(vsifile, 0, 2)
vsileng = gdal.VSIFTellL(vsifile)
gdal.VSIFSeekL(vsifile, 0, 0)
return gdal.VSIFReadL(1, vsileng, vsifile)
def write_vsimem(fn,data):
'''Write GDAL vsimem files'''
vsifile = gdal.VSIFOpenL(fn,'w')
size = len(data)
gdal.VSIFWriteL(data, 1, size, vsifile)
return gdal.VSIFCloseL(vsifile)
infn='test.tif'
outfn='test2.tif'
#List of bands to retain, output bands will be reordered to match
bands=[4,3,2]
ds = gdal.Open(os.path.abspath(infn)) #Ensure path is absolute not relative path
vrtdrv=gdal.GetDriverByName('VRT')
tifdrv=gdal.GetDriverByName('GTIFF')
#Create an in memory VRT copy of the input raster
vfn='/vsimem/test.vrt'
vrtds=vrtdrv.CreateCopy(vfn,ds)
#Read the XML from memory and parse it
#Could also write the VRT to a temp file on disk instead of /vsimem
#and used etree.parse(vrtfilepath) instead
#i.e.
#vfn='some/path/on/disk/blah.vrt'
#vrtds=vrtdrv.CreateCopy(vfn,ds)
#tree=etree.parse(vfn)
#os.unlink(vfn)
vrtds=vrtdrv.CreateCopy(vfn,ds)
vrtxml=read_vsimem(vfn)
tree=etree.fromstring(vrtxml)
#Ignore bands not in band list
#And put bands we want to keep in the correct order
for band in tree.findall('VRTRasterBand'):
try:
i=bands.index(int(band.attrib['band']))
band.attrib['band']=str(i+1)
bands[i]=band
except ValueError:pass
finally:
tree.remove(band)
for band in bands:
tree.append(band)
#Get the XML string from the tree
vrtxml=etree.tostring(tree, pretty_print=True)
#Write it to the in memory file
#Could also just write to a file on disk
write_vsimem(vfn,vrtxml)
#Open the VRT containing only the selected bands
vrtds=gdal.Open(vfn)
#Write to a new raster
tifds=tifdrv.CreateCopy(outfn,vrtds, 0, ['COMPRESS=JPEG','PHOTOMETRIC=YCBCR','TFW=YES'])
#Close dataset to ensure cache is properly flushed
del tifds
I complemented the code with a few lines, to set spatial reference and same geotransform parameters of original rasters, and it works perfectly in a Linux system (python 2.7.x and GDAL 1.11.x). In a Windows 8.x system (python 2.7.x and GDAL 1.11.x), it too worked but it printed these console errors (first one is similar to your error but I got a resulting raster):
with next command:
python ndvi_github.py b3.ND.tif b4.ND.tif output-ndvi.tif
Apparently, resulting raster hasn't any problem.
Complete used code was:
#!/usr/bin/env/python
# ndvi.py red.tif nir.tif output-ndvi.tif
# Calculate NDVI (see Wikipedia). Assumes atmospheric correction.
# (Although I use it without all the time for quick experiments.)
import numpy as np
from sys import argv
from osgeo import gdal, gdalconst, osr
#type for internal calculations
t = np.float32
red, nir = map(gdal.Open, argv[1:3])
print argv[1:3]
transform = red.GetGeoTransform()
geotiff = gdal.GetDriverByName('GTiff')
output = geotiff.CreateCopy(argv[3], red, 0)
output = geotiff.Create(argv[3],
red.RasterXSize,
red.RasterYSize,
1,
gdal.GDT_Float32)
# Ugly syntax, but fast:
r = red.GetRasterBand(1).ReadAsArray(0,0,red.RasterXSize,red.RasterYSize)
n = nir.GetRasterBand(1).ReadAsArray(0,0,nir.RasterXSize,nir.RasterYSize)
# Convert the 16-bit Landsat 8 values to floats for the division operation:
r = r.astype(t)
n = n.astype(t)
# Tell numpy not to complain about division by 0:
np.seterr(invalid='ignore')
# Here's the meat of this whole thing, the actual NDVI formula:
ndvi = (n - r)/(n + r)
# The ndvi value is in the range -1..1, but we want it to be displayable, so:
# Make the value positive and scale it back up to the 16-bit range:
ndvi = (ndvi + 1) * (2**15 - 1)
# And do the type conversion back:
ndvi = ndvi.astype(np.uint16)
output.GetRasterBand(1).WriteArray(ndvi)
output.SetGeoTransform(transform)
wkt = red.GetProjection()
#setting spatial reference of output raster
srs = osr.SpatialReference()
srs.ImportFromWkt(wkt)
output.SetProjection( srs.ExportToWkt() )
red = None
nir = None
output = None
and next image was resulting raster, obtained in Windows 8 system, loaded in QGIS.
Best Answer
I would use gdal_translate and gdal_merge.py:
Translate the images to crop (subset them) using gdal_translate. You could use a bash script to automate. Something along the lines below.
for f in *.tif;do gdal_translate -projwin ulx uly lrx lry "$f" "$f".cropped.tif ; doneUse gdal_merge.py to 'stack' the images. Here we are not explicitly controlling stack order. I believe these will stack alphanumerically in ascending order (0-9, then a-z).
gdal_merge.py -separate -o myoutput.tif *.cropped.tifI used -o myoutput.tif because I do not know if myoutput.vrt will work with gdal_merge.py. I assume it would, as it is a GDAL supported format, but I have never tested it.