Build a table with columns of raster values from multiple raster datasets, using Python, GDAL, or PyQGIS?

Question

I have a series of raster datasets, clipped to the same extent, with matching cell resolutions, each of which represent land surface temperatures (LST) captured by the MODIS sensor. Each raster stores LST values for sequential time intervals (daily or 8-day intervals) for a year.

I need to build a table with the following specifications:

1. each row represents a grid cell
2. one column with a unique grid cell ID
3. each subsequent column stores LST values from rasters for each time interval (so up to 365 columns max.)

I realize I could do this easily with ArcPy by creating grid cell centroids from one of the rasters, then use the Multi Value to Point tool to build an attribute table for the resulting point dataset with a field for each LST raster. However, I am building this tool as a PyQGIS plugin, and so need to do this with FOSS libraries and resources (i.e. GDAL/OGR, PyQGIS, etc.)

At this point, I'm a bit unsure as to the best method for doing this. Here are the alternative methods I've considered:

1. Convert all of the rasters to xyz files, then use Python's built-in csv library to merge the xyz files together. I like this option, but I'm still a Python novice, and I haven't been able to figure out how to iteratively merge csv columns together.
2. Create a point shapefile from grid cell centroids, then try to incorprate Borys Jurgiel's Point sampling tool plugin to extract the values from all of the remaining rasters. This tool is written as a plugin however, and would take quite a bit of work to fit into my existing code.
3. Turn the XYZ value into a 1-D Numpy array, and do some kind of Numpy-related magic? (I'm not very experienced with Numpy at this point...)
4. Use a RDBMS and SQL to do a relational join between csv files. I'm familiar with PostGIS and SpatiaLite/SQLite, but this will eventually be part of a QGIS plugin, so I'd raster not force users to have to install an RDBMS just to fulfill this one task in my workflow.

Are there any other suggestions as to how I should attack this? The first option (joining csv files) seems like the best, but maybe I'm over thinking this problem.

UPDATE Just to follow up on this issue, I decided to go with option 1. I convert all of the rasters to xyz file, using gdal2xyz.py (included with standard GDAL installations), then I add a header row and unique id column for each xyz file (stored as a CSV). The final solution to merge these files together into a single CSV table can be found at StackOverflow, here.

Answer 1

Your first idea seems to me good. I would use cut and paste Linux command and a bash script on Linux and gdal to convert raster to sorted ASCII gridded xyz file.

If you would like to use Python use Pandas to handle csv files.