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:
where the my_files contains a list of the files you want to use, in the specific order. Then crop with gdal_translate:
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.
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).
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.