PostGIS导出指定范围内的栅格数据

针对已经通过PostGIS导入到PostgreSQL中的栅格数据，给定经纬度范围，实现栅格数据的导出。
查询SQL语句如下：

SELECT ST_Union(ST_Clip(rast,geom)) AS rast FROM public.alt CROSS JOIN
ST_MakeEnvelope(90,60,120,45,4326) As geom WHERE ST_Intersects(rast,geom);

其中PostGIS函数的含义如下：

ST_MakeEnvelope：函数用于构造一个矩形范围，其参数分别是最小Ｘ值，最小Ｙ值，最大Ｘ值，最大Ｙ值和坐标系代码 
ST_Intersects：函数用于选择出与geom矩形相交的栅格Tiles
ST_Clip：函数用于将选择出来的Tiles进行裁剪，得到geom范围的数据
ST_Union：函数用于聚合选择出来的数据为一个整体

导出TIF格式的SQL语句如下：

SELECT ST_AsTIFF(rast, 'LZW') FROM (SELECT ST_Union(ST_Clip(rast,geom)) AS rast FROM public.alt
CROSS JOIN ST_MakeEnvelope(90,60,120,45,4326) As geom WHERE ST_Intersects(rast,geom)) AS rasttiff;

完整的Python代码如下：

import psycopg2

# Connect to an existing database
conn = psycopg2.connect('host=localhost port=5432 user=postgres password=123456 dbname=rastertest')

# Open a cursor to perform database operations
cur = conn.cursor()
strsql = "SELECT ST_AsTIFF(rast, 'LZW') " 
         "FROM (" 
         "SELECT ST_Union(ST_Clip(rast,geom)) AS rast " 
         "FROM " 
         "alt " 
         "CROSS JOIN " 
         "ST_MakeEnvelope(90,60,120,45,4326) As geom " 
         "WHERE ST_Intersects(rast,geom)" 
         ") AS rasttiff"

cur.execute(strsql)

# Fetch data as Python objects
rasttiff = cur.fetchone()

# Write data to file
if rasttiff is not None:
    open('E:/personalfile/cs1/alt_ex1.tif', 'wb').write(str(rasttiff[0]))

# Close communication with the database
cur.close()
conn.close()

具体裁切导出效果如下图所示：

影像导出效果图.png

参考原文链接：
http://blog.csdn.net/theonegis/article/details/55211846