PostGIS导入导出栅格数据

刺骨的言语ヽ痛彻心扉 2022-07-12 05:40 501阅读 0赞

上一篇博文[PostGIS导入导出ESRI Shapefile数据][PostGIS_ESRI Shapefile]介绍了如何导入空间矢量数据到PostgreSQL中，紧接上一篇，本文将介绍如何使用PostGIS导入导出空间栅格数据。  
实验环境和上一篇的博文一样： Ubuntu 16.04 LTS 。  
数据为全球影像（ESRI ArcGIS提供的示例数据），下载连接：[全球影像百度网盘下载][Link 1]  
可以使用GDAL的gdalinfo命令查看其详细信息：  
![gdalinfo][]  
在QGIS中查看如下：  
![wsiearth.tif][]

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

![PostGIS插入栅格数据][PostGIS]  
可以使用psql进入数据库查看：  
![PostGIS插入栅格数据][PostGIS 1]  
还可以使用`\d staging.wsiearth`查看关系表的结构：  
![PostGIS中的栅格关系表][PostGIS 2]  
可以看到有一个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_postgis_gdal_enabled_drivers.html]  
[http://postgis.net/docs/manual-2.2/postgis\_enable\_outdb\_rasters.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][http_initd.org_psycopg_docs_index.html]。  
我这里直接贴出代码，因为代码很容易理解。ST\_AsTIFF函数参见：[http://postgis.net/docs/RT\_ST\_AsTIFF.html][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中查看导出结果如下：  
![PostGIS导出结果][PostGIS 3]  
可以发现导出的结果，丢失了图像的colormap，但是每个像素值都是正确的。

[PostGIS_ESRI Shapefile]: http://blog.csdn.net/theonegis/article/details/55046419
[Link 1]: https://pan.baidu.com/s/1gfxE5OJ
[gdalinfo]: /images/20220711/1faa29dea8b848d6a46a5cc5351eabff.png
[wsiearth.tif]: /images/20220711/7428df20c8674532993f62c940c73570.png
[PostGIS]: /images/20220711/82e2dd772b2c48f5ae80b7075f325779.png
[PostGIS 1]: /images/20220711/309e46b4e88944778980c789bd09e96c.png
[PostGIS 2]: /images/20220711/a72839899b6f48699ca661447496cee1.png
[http_postgis.net_docs_postgis_gdal_enabled_drivers.html]: http://postgis.net/docs/postgis_gdal_enabled_drivers.html
[http_postgis.net_docs_manual-2.2_postgis_enable_outdb_rasters.html]: http://postgis.net/docs/manual-2.2/postgis_enable_outdb_rasters.html
[http_initd.org_psycopg_docs_index.html]: http://initd.org/psycopg/docs/index.html
[http_postgis.net_docs_RT_ST_AsTIFF.html]: http://postgis.net/docs/RT_ST_AsTIFF.html
[PostGIS 3]: /images/20220711/67096b6b3415419bb0fc961ffd608634.png