PostGIS批量导入栅格数据

绝地灬酷狼 2022-07-12 06:20 618阅读 0赞

如果我们有一批以文件存储的影像数据如何利用PostGIS批量的导出到PostgreSQL数据库中进行管理呢？  
（单个数据的导入参见我上篇博文：[PostGIS导入导出栅格数据][PostGIS]）

我的实验环境如下：  
OS: Ubuntu 16.04 LTS  
PostgreSQL:9.5.5  
PostGIS: 2.2

我的实验数据是全球气象数据（不过，我只下载了中国部分19区的平均气温数据（TIF格式）做示例演示）  
下载网址：[http://www.worldclim.org/tiles.php][http_www.worldclim.org_tiles.php]  
下载好的数据进行解压可以看到包含了12个月的平均气温数据：  
![数据准备][SouthEast]  
使用GDAL的命令行工具gdalinfo查看其元数据，可以看到其采用的是WGS84坐标系，大小是3600×3600。  
![gdalinfo查看元数据][gdalinfo]  
我们可以在raster2pgsql命令中使用通配符批量的导入数据到PostgreSQL数据库，命令如下：  
（具体参数可自行查看命令的帮助文档）

raster2pgsql -s 4326 -I -C -M ./*.tif -F -t 256x256 staging.tmean_19 | psql -h localhost -p 5432 -U postgres -d postgis_in_action -W

需要说明的是：  
`-t`参数指定了数据的分块大小，原始的3600×3600以256×256的Tile进行划分存储。  
`-F`参数在关系表中增加名为filename的列用于保存原始文件的名称，这在我们进行Tile的合并的时候很有用。  
![raster2pgsql导入数据][raster2pgsql]  
进入数据库可以查看我们导入的数据表  
![栅格数据导入关系表][SouthEast 1]  
通过SQL进行查看，发现每个文件被分割成了225块。（但是(3600\*3600)/(256\*256)应该是197，不是225，我也不知道为什么被分割成了225块？？？）  
![导入影像分块数目][SouthEast 2]

--------------------

最后，说说怎么导出数据到文件。和上一篇博文一样，我们还是使用Python脚本进行导出。这里主要用到了两个函数ST\_AsTIFF和ST\_Union。ST\_AsTIFF上一篇中我们已经使用过，用于将PostgreSQL中的raster类型导出为TIFF格式对应的bytea类型，而ST\_Union对于栅格数据来说，用于将分块的影像重新拼接为完整的一副影像。  
详细代码如下：

# -*- 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;")
    cur.execute("SELECT ST_AsTIFF(ST_Union(rast), 'LZW') AS rasttiff FROM staging.tmean_19 WHERE filename='tmean1_19.tif';")
    
    # 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中查看和原来的影像一模一样：  
![PostGIS导出结果][PostGIS 1]

[PostGIS]: http://blog.csdn.net/theonegis/article/details/55050943
[http_www.worldclim.org_tiles.php]: http://www.worldclim.org/tiles.php
[SouthEast]: /images/20220711/337908ff81cc49759e3cd2a1c7133313.png
[gdalinfo]: /images/20220711/19e60b46c15c48d6ad7e713c4cd16d51.png
[raster2pgsql]: /images/20220711/d3d3141cead0411095c35e9ae67f16a0.png
[SouthEast 1]: /images/20220711/0e1fb70da5bf4d19a195e97846b85750.png
[SouthEast 2]: /images/20220711/af1d46f8b1c14ce08f2aa851da9fc6e8.png
[PostGIS 1]: /images/20220711/abc17ccf50464cee89e99fa44c02ac92.png