PostGIS查询指定范围的数据

对于上一篇PostGIS批量导入栅格数据中导入的气温数据，如何查询指定范围的气温呢？
比如，给定了经纬度范围，如何取出给定月份的数据？
下面的SQL代码给出了查询方法：

SELECT ST_Union(ST_Clip(rast,geom)) AS rast FROM staging.tmean_19 CROSS JOIN ST_MakeEnvelope(3.87,73.67,53.55,135.05,4326) As geom WHERE ST_Intersects(rast,geom) AND month=1;

其中，
ST_MakeEnvelope函数用于构造一个矩形范围，其参数分别是最小Ｘ值，最小Ｙ值，最大Ｘ值，最大Ｙ值和坐标系代码；
ST_Intersects函数用于选择出与geom矩形相交的栅格Tiles；
ST_Clip函数用于将选择出来的Tiles进行裁剪，得到geom范围的数据；
ST_Union函数用于聚合选择出来的数据为一个整体；
上述的SQL返回的结果是raster类型的数据，如果想要将结果导出为TIFF格式的数据，SQL代码如下：

SELECT ST_AsTIFF(rast,'LZW') FROM ( SELECT ST_Union(ST_Clip(rast,geom)) AS rast FROM staging.tmean_19 CROSS JOIN ST_MakeEnvelope(97.51,37.28,111.55,50.52,4326) As geom WHERE month=1 AND ST_Intersects(rast,geom) ) AS rasttiff;

完整的Python代码如下：

# -*- coding: utf-8 -*-

import psycopg2

# Connect to an existing database
conn = psycopg2.connect('host=localhost port=5432 user=postgres password=post1231 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;")
# cur.execute("SELECT ST_AsTIFF(ST_Union(rast),'LZW') AS rasttiff FROM staging.tmean_19 WHERE filename='tmean1_19.tif';")
strsql = "SELECT ST_AsTIFF(rast,'LZW') "  "FROM ("  "SELECT ST_Union(ST_Clip(rast,geom)) AS rast "  "FROM "  "staging.tmean_19 "  "CROSS JOIN "  "ST_MakeEnvelope(97.51,4326) As geom "  "WHERE month=1 AND 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('/home/theone/Desktop/tmean1.tif','wb').write(str(rasttiff[0])) # Close communication with the database cur.close() conn.close()

我们可以在QGIS中查看结果，并和原图进行对比：

（编辑：李大同）

【声明】本站内容均来自网络，其相关言论仅代表作者个人观点，不代表本站立场。若无意侵犯到您的权利，请及时与联系站长删除相关内容!