gdal读取postgis数据库中的栅格数据（python）

---

前言

gdal读取tif格式的栅格影像非常方便，它还支持读取postgis的栅格影像。不过，一开始我以为读取方式与tif会有较大差异，所以就上网搜索gdal如何读取postgis栅格。很不幸，内网几乎没有，外网没查过，一是不会魔法上网，二是硬看英文太累。功夫不负有心人，最后还是找到了利用gdal读取postgis栅格的JAVA实现（https://blog.youkuaiyun.com/qq_24309981/article/details/84184772）。由于不同语言下，gdal的函数（方法）名称几乎是一致的，所以这也算解决了我的问题。

---

一、向postgis导入栅格

读取之前的让你的数据库有栅格影像。因为gdal暂不支持向postgis导入栅格，因此得用其他方法------最简单的就是用postgis自带的程序raster2pgsql。具体使用在这里不赘述，本人学习的是这篇文章：https://blog.youkuaiyun.com/theonegis/article/details/55050943
这篇文章在内网被转发了好多次，我也不知道谁是原创了。

二、gdal读取postgis栅格

其实很简单，就是文件路径跟读取tif时设置不一样，其他都一样。

 \\tif的情形
 path=r'D:abc\def\xxx.tif'
 \\postgis的情形
 path="PG:host=localhost port=5432 dbname='mydb' user='postgres' password='secret' schema='public' table=mytable"
 \\实际上就是表示文件路径的字符串改成了访问数据库的所需的字符串。上述字符串的具体含义参考gdal官方手册的postgisraster章节:
 \\https://www.osgeo.cn/gdal/drivers/raster/postgisraster.html#raster-postgisraster
 \\然后像打开tif一样读取postgis栅格即可
 from osgeo import gdal
 gdal.AllRegister()
 dataset=gdal.Open(path,gdal.GA_ReadOnly)
 print(dataset)

直接打印将返回一个栅格对象：

<osgeo.gdal.Dataset; proxy of <Swig Object of type 'GDALDatasetShadow *' at 0x00000274FAB24360> >

三、遇到问题

上述方法仍有bug。当栅格数据被手动分块（瓦片）或自动分块（常见于数据较大的情况）存进postgis数据库时，栅格在数据库中表现为一张多行的表。
1.栅格数据较小，或没手动分块时，数据库中的栅格表是这样的：

rid	物理地址（猜的，没看具体说明文档）
1	00011S021FSDFFS…（一长串数字及英文字符）

2.栅格较大或手动分块时，数据库中的栅格表是这样的：

rid	。。。。
1	xxxxx
2	xxxxxx
。。。	。。。。
n	xxxxxx

假如是第一种情况，上述读取方式能正确读取数据库中的栅格，但是第二种情况下，上述方法不能正确读取。具体情况为：栅格的投影，行列数，波段数等都为0或者空，读取栅格矩阵时，返回的也是None。
查阅gdal官方文档后发现，上述读取方式中还有一个“mode”的参数，该参数有两种情况：1.这是默认设置，认为数据表中一行表示一张影像；2.认为一张表的所有行同属一张影像。
显然，出现bug的原因是“mode”参数没指定，gdal默认使用mode=1的情形，然而我们的数据表其实是一张影像，因此在一张影像分多行存储时，gdal没能读取完整，从而导致bug出现。
因此，在影像分块后，记得在上述path字符串中加上“mode=2”，gdal就能正常读取了。