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?
I would first 'stack' them:
gdalbuildvrt -separate -input_file_list my_files.txt my.vrt
where the my_files contains a list of the files you want to use, in the specific order. Then crop with gdal_translate:
gdal_translate -projwin ulx uly lrx lry my.vrt my.tif
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 ; done
Use 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.tif
I 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.