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函数用于构造一个矩形范围,其参数分别是最小X值,最小Y值,最大X值,最大Y值和坐标系代码;
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,37.28,111.55,50.52,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中查看结果,并和原图进行对比:
转载自:https://blog.csdn.net/theonegis/article/details/55211846