上一篇文章PostGIS导入导出ESRI Shapefile数据（http://www.linuxdiyf.com/linux/28395.html）介绍了如何导入空间矢量数据到PostgreSQL中，紧接上一篇，本文将介绍如何使用PostGIS导入导出空间栅格数据。
实验环境和上一篇的文章一样：Ubuntu 16.04 LTS。
数据为全球影像（ESRI ArcGIS提供的示例数据），下载连接：https://pan.baidu.com/s/1gfxE5OJ
可以使用GDAL的gdalinfo命令查看其详细信息：
在QGIS中查看如下：
PostGIS提供了raster2pgsql工具用于栅格数据的导入（可以使用man raster2pgsql命令查看帮助文档）。
使用如下命令进行影像数据的插入，具体参数的含义这里不在累赘，上篇文章有说明。关于raster2pgsql的参数可以使用man命令进行查看。
raster2pgsql -s 4326 -C ~/Downloads/gisdata/wsiearth.tif staging.wsiearth | psql -h localhost -p 5432 -U postgres -d postgis_in_action -W
可以使用psql进入数据库查看：
还可以使用\d staging.wsiearth查看关系表的结构：
可以看到有一个rid和rast的列，rid是以一个整形对插入的栅格数据进行的标示，rast列的类型是raster存储了具体数据。
如果想要导出数据，可以使用PostGIS提供的内置函数。其中ST_AsGDALRaster 是一个通用性的函数，可以导出为GDAL支持的任意格式。此外，还有ST_AsPNG，ST_AsJPEG和ST_AsTIFF等函数，提供了直接导出到指定格式，并且接受更少的参数。
注意：PostGIS的默认导出GDAL支持的格式数据是关闭的，启用的话需要修改PostgreSQL的配置文件。修改的方法比较多，而且不同的PostgreSQL版本的支持情况不一样。我使用如下命令进行修改：
sudo vim /etc/postgresql/9.5/main/postgresql.conf
在该配置文件后面添加：
postgis.gdal_enabled_drivers = 'ENABLE_ALL'
postgis.enable_outdb_rasters = True
重启PostgreSQL：
sudo service postgresql restart
参考链接：
http://postgis.net/docs/postgis_gdal_enabled_drivers.html
http://postgis.net/docs/manual-2.2/postgis_enable_outdb_rasters.html
第二个问题是：这些内置函数提供的导出结果是PostgreSQL的bytea (byte array)数据类型，我们需要自己写程序进行转换成实际的图像文件。
下面我使用Python的Psycopg库连接PostgreSQL数据库，进行查询并导出最终的结果。Psycopg库参见：http://initd.org/psycopg/docs/index.html。
我这里直接贴出代码，因为代码很容易理解。ST_AsTIFF函数参见：http://postgis.net/docs/RT_ST_AsTIFF.html。
# -*- coding: utf-8 -*-
import psycopg2
# Connect to an existing database
conn = psycopg2.connect('host=localhost port=5432 user=postgres password=password dbname=postgis_in_action')
# Open a cursor to perform database operations
cur = conn.cursor()
# Execute SQL query
cur.execute("SELECT ST_AsTIFF(rast, 'LZW') AS rasttiff FROM staging.wsiearth WHERE rid=1;")
# Fetch data as Python objects
rasttiff = cur.fetchone()
# Write data to file
if rasttiff is not None:
open('/home/theone/Desktop/wsiearth.tif', 'wb').write(str(rasttiff[0]))
# Close communication with the database
cur.close()
conn.close()
在QGIS中查看导出结果如下：
可以发现导出的结果，丢失了图像的colormap，但是每个像素值都是正确的。