针对已经通过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:函数用于构造一个矩形范围,其参数分别是最小X值,最小Y值,最大X值,最大Y值和坐标系代码
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