In general, it seems there are roughly five (5) ways to get "file data" (e.g. a GeoTIFF) out of a PostGIS geoprocessing workflow:
- Export just the raster field as an ASCII grid
- Connect to the database using a desktop client (e.g. QGIS) 
- Use a procedural language (like PLPGSQL or PLPYthon) 
- Use the COPY declaration to get a hex dump out and convert it to a binary file
- Fill a 2D NumPy array with a byte array and serialize it to a binary file using GDAL or psycopg2 [3, 4]
- Use ST_AsTiff() or the more general ST_AsGDALRaster() to get a byte array, which can be written to a binary file
Out of all of these, I prefer the fifth option, using a Python programming environment to execute the query, receive the byte array, and serialize it to a file, particularly because the byte array buffer can be delivered to a web framework as a file-attachment response. I will briefly review some of the other options first.
 Obe, Regina. Leo S. Hsu. "Using PostGIS in a desktop environment." Chapter 12. PostGIS in Action. 2011.
Exporting as an ASCII Grid (Text File)
This works great! Because an ASCII grid file (or "ESRI Grid" file, with the *.asc or *.grd extension, typically) is just plain text, you can directly export it from the database. The GDAL driver name is "AAIGrid" which should be the second argument to ST_AsGDALRaster(). Be sure to remove the column header from your export (see image below). However, what you get is a file that has no projection information that you may need to convert to another format. This can present problems for your workflow, especially if you're trying to automate the production of raster files, say, through a web API.
Connect Using QGIS Desktop Client
There is a plug-in for QGIS that promises to allow you to load raster data from PostGIS directly into a QGIS workspace. I used the Plugins Manager ("Plugins" > "Fetch Python Plugins...") in QGIS to get this plug-in package. The first time I selected the "Load PostGIS Raster to QGIS" plug-in and tried to install it, I found that I couldn't write to the plug-ins directory (this with a relatively fresh installation of QGIS). After creating and setting myself as the owner of the python/plugins directory, I was able to install the plug-in without any further trouble. Connecting to the database and viewing the available relations was also no trouble at all. One minor irritation is that you need to enter your password every time the plug-in interfaces with the database, which can be very often, at every time the list of available relations needs to be updated.
There are a few options available to you in displaying raster data from the database: "Read table's vector representation," "Read one table as a raster," "Read one row as a raster," or "Read the dataset as a raster." It's not clear what the second and last choices are, but "Reading the table as a raster" did not work for me where my table has one raster field and a couple of non-raster, non-geometry/geography fields; QGIS hung for a few seconds then said it "Could not load PG..." Reading one row worked, however, you have to select the row by its primary key (or row number in a random selection, not sure which it is returning). This may not be the easiest way for you to find a particular raster image in a table.
Using the COPY Declaration in SQL
pip install psycopg2 pygresql
The basic idea is to use the
COPY declaration in SQL to export the raster to a hexadecimal file, then to convert that file to a binary file using
import os, stat, pg # See: http://www.pygresql.org/install.html # pip install psycopg2, pygresql # Designate path to output file outfile = '/home/myself/temp.tiff' # Name of PostgreSQL table to export pg_table = 'geowepp_soil' # PostgreSQL connection parameters pg_server = 'my_server' pg_database = 'my_database' pg_user = 'myself' # Desginate a file to receive the hex data; make sure it exists with the right permissions pg_hex = '/home/myself/temp.hex' os.mknod(pg_hex, stat.S_IRUSR | stat.S_IWUSR | stat.S_IRGRP | stat.S_IWGRP) conn = pg.connect(pg_database, pg_server, 5432, None, None, pg_user) sql = "COPY (SELECT encode(ST_AsTIFF(ST_Union(" + pg_table + ".rast)), 'hex') FROM " + pg_table + ") TO '" + pg_hex + "'" # You can check it with: print sql conn.query(sql) cmd = 'xxd -p -r ' + pg_hex + ' > ' + outfile os.system(cmd)
This needs to be done on the file system of the database server, which is where PostgreSQL will write.
Using ST_AsTIFF() and Serializing from a Byte Array Buffer
Despite the seeming complexity of this option, I think it is the most straightforward and flexible approach. I'll provide two examples here, with code: using GeoDjango to execute a raw query and using GeoAlchemy2's object-relational model to execute the query. Finally, I'll show an example of writing the output to a file or to a Django
First, some setup. We'll define a
RasterQuery class to help with handling the details. While a new class isn't exactly an idiomatic example, I'm hoping it will succinctly illustrate the considerations involved in using performing raw SQL queries with Django.
class RasterQuery: ''' Assumes some global FORMATS dictionary describes the valid file formats, their file extensions and MIME type strings. ''' def __init__(self, qs, params=None, file_format='geotiff'): assert file_format in FORMATS.keys(), 'Not a valid file format' self.cursor = connection.cursor() self.params = params self.query_string = qs self.file_format = file_format self.file_extension = FORMATS[file_format]['file_extension'] def execute(self, params=None): '''Execute the stored query string with the given parameters''' self.params = params if self.params is not None: self.cursor.execute(self.query_string, params) else: self.cursor.execute(self.query_string) def fetch_all(self): '''Return all results in a List; a List of buffers is returned''' return [ row for row in self.cursor.fetchall() ] def write_all(self, path, name=None): '''For each raster in the query, writes it to a file on the given path''' name = name or 'raster_query' i = 0 results = self.fetch_all() for each in results: name = name + str(i + 1) + self.file_extension with open(os.path.join(path, name), 'wb') as stream: stream.write(results[i]) i += 1
RasterQuery class available to us, we can more cleanly execute our raw SQL queries and serialize the response to a file attachment in a Django view:
def clip_one_raster_by_another(request): # Our raw SQL query, with parameter strings query_string = ''' SELECT ST_AsGDALRaster(ST_Clip(landcover.rast, ST_Buffer(ST_Envelope(burnedarea.rast), %s)), %s) AS rast FROM landcover, burnedarea WHERE ST_Intersects(landcover.rast, burnedarea.rast) AND burnedarea.rid = %s''' # Create a RasterQuery instance; apply the parameters query = RasterQuery(query_string) query.execute([1000, 'GTiff', 2]) filename = 'blah.tiff' # Outputs: # [(<read-only buffer for 0x2613470, size 110173, offset 0 at 0x26a05b0>), # (<read-only buffer for 0x26134b0, size 142794, offset 0 at 0x26a01f0>)] # Return only the first item response = HttpResponse(query.fetch_all(), content_type=FORMATS[_format]['mime']) response['Content-Disposition'] = 'attachment; filename="%s"' % filename return response
Seem simple enough? To write to a file instead, see the
write_all() method of the
RasterQuery class. The
query.fetch_all() at the end is contrived. I'll show a better way of getting to a nested buffer in the next example.
GeoAlchemy2's object-relational model (ORM) allows tables to be represented as classes, just like in Django.
class LandCover(DeclarativeBase): __tablename__ = 'landcover' rid = Column(Integer, primary_key=True) rast = Column(ga2.types.Raster) filename = Column(String(255)) class BurnedArea(DeclarativeBase): __tablename__ = 'burnedarea' rid = Column(Integer, primary_key=True) rast = Column(ga2.types.Raster) filename = Column(String(255)) burndate = Column(Date) burnname = Column(String(255))
ENGINE global variables are available, the gist of this approach can be seen in this example:
def clip_fccs_by_mtbs_id2(request): query = SESSION.query(LandCover.rast\ .ST_AsGDALRaster(LandCover.rast\ .ST_Clip(LandCover.rast, BurnedArea.rast\ .ST_Envelope()\ .ST_Buffer(1000)), 'GTiff').label('rast'))\ .filter(LandCover.rast.ST_Intersects(BurnedArea.rast), BurnedArea.rid==2) filename = 'blah.tiff' # Outputs: # [(<read-only buffer for 0x2613470, size 110173, offset 0 at 0x26a05b0>), # (<read-only buffer for 0x26134b0, size 142794, offset 0 at 0x26a01f0>)] result = query.all() while type(result) != buffer: result = result # Unwrap until a buffer is found # Consequently, it returns only the first item response = HttpResponse(result, content_type=FORMATS[_format]['mime']) response['Content-Disposition'] = 'attachment; filename="%s"' % filename return response
Here we see a better way of getting at a nested
buffer. If we wanted all of the rasters that were returned (all of the buffers), we could call
ST_Union on our final raster selection before passing it to