Ganos低代码实现免切片遥感影像浏览（一）：金字塔_云原生数据库 PolarDB(PolarDB)-阿里云帮助中心

使用PolarDB PostgreSQL版数据库，配合Ganos时空数据库引擎，在不借助第三方工具的情况下，可仅使用SQL语句快速管理与展示遥感影像数据。Ganos共提供两种影像免切浏览的方法，一种借助Ganos Raster插件，使用窗口范围获取影像数据展示，另一种通过固定瓦片范围获取影像数据展示，本文详细介绍第一种方法并提供前后端实操代码帮助您快速理解Ganos Raster的使用细节。

背景

当管理着日益增长的遥感数据时，想要将任一影像展现在地图上往往要思考如下问题：

对于可能只使用几次的影像进行传统的切片、发布流程是否兴师动众？瓦片存储在何处？瓦片数据如何进行后续的管理？
不进行预切片而采用实时切片的方案能否满足前端的响应需求？遇到特别巨大的影像时，实时切片会受到多大影响？

不过，使用PolarDB PostgreSQL版管理影像时，借助Ganos Raster扩展，无需依赖第三方工具，仅通过SQL语句即可高效快速地将数据库中的影像展示在地图上。

最佳实践

前期准备

导入影像数据

安装ganos_raster插件。
```
CREATE EXTENSION ganos_raster CASCADE;
```

创建包含Raster属性列的测试表raster_table。

CREATE TABLE raster_table (ID INT PRIMARY KEY NOT NULL,name text,rast raster);

从OSS中导入一幅影像作为测试数据。使用ST_ImportFrom函数导入影像数据到分块表chunk_table，详细的语法介绍请参考ST_ImportFrom。
```
INSERT INTO raster_table VALUES (1, 'xxxx影像', ST_ImportFrom('chunk_table','oss://<access_id>:<secrect_key>@<Endpoint>/<bucket>/path_to/file.tif'));
```

创建金字塔

金字塔是快速浏览影像数据的基础。对于新导入的数据，建议首先创建金字塔。Ganos提供ST_BuildPyramid函数以实现金字塔的创建。详细的语法介绍请参考ST_BuildPyramid。

UPDATE raster_table SET rast = st_buildpyramid(raster_table,'chunk_table') WHERE name = 'xxxx影像';

语句解析

创建金字塔成功后，可以使用Ganos的ST_AsImage函数从数据库中获取指定范围的影像。以下为ST_AsImage函数基础语法介绍，详细介绍请参考ST_AsImage。

bytea ST_AsImage(raster raster_obj,
        box extent,
        integer pyramidLevel default 0,
        cstring bands default '',
        cstring format default 'PNG',
        cstring option default '');

ST_AsImage函数的参数可以分为两类：静态参数与动态参数。

静态参数

一般不随操作而变化的参数，可以在代码中固定，减少重复工作。

bands：需要获取的波段列表。
可以使用'0-3'或'0,1,2,3'等形式，不能超过影像本身的波段。
format：输出的影像格式。
默认为PNG格式，可选JPEG格式。由于PNG格式在数据压缩效果上不及JPEG格式的有损压缩，因此在传输时将消耗更多的时间。在没有透明度需求的情况下，建议优先选择JPEG格式。
option：JSON字符串类型的转换选项。可定义额外的渲染参数。

动态参数

随操作而变化，需要动态生成。

extent：需要获取的影像范围。
在相同条件下，显示范围越大，数据库的处理时间越长，返回的图片体积也相应增大，从而导致总体响应时间延长。因此，建议仅获取用户视域范围内的影像，以确保传输效率。
pyramidLevel：使用的影像金字塔层级。
使用的金字塔级别越高，图像的清晰度越高，但其体积也随之增加。因此，需要选择最合适的金字塔级别，以确保传输效率。

获取影像外包框范围

使用Ganos的ST_Envelope函数获取影像的外包框范围，然后使用ST_Transform函数将影像转换到常用坐标系（此处使用WGS 84坐标系），最终将其转换为前端方便使用的格式。

SELECT replace((box2d(st_transform(st_envelope(geom),4326)))::text,'BOX','') FROM rat_clip WHERE name = 'xxxx影像';

获取最适宜金字塔级别

使用Ganos的ST_BestPyramidLevel函数计算特定影像范围内最适宜的金字塔级别。以下为ST_BestPyramidLevel函数基础语法介绍，详细介绍请参考ST_BestPyramidLevel。

integer ST_BestPyramidLevel(raster rast, 
                            Box extent, 
                            integer width, 
                            integer height)

其中：

extent：可视区域所要获取的影像范围。与ST_AsImage函数中使用的范围一致。
width/heigth：可视区域的像素宽高。一般情况下，即使用的前端地图框的尺寸。

需要注意的是，ST_BestPyramidLevel和ST_AsImage函数中使用的是原生的Box类型，而非Geometry兼容类型，因此需要将其转换。前端传回的是BBox数组，要将其转换为原生Box类型需要如下步骤：

在SQL语句中将前端传回的字符串构建为Geometry类型对象。
将Geometry类型对象转换为Text类型对象。
使用文本替换方式将其转换为Box类型兼容格式的Text类型对象。
将Text类型对象转换为原生Box类型对象。

SELECT  Replace(Replace(Replace(box2d(st_transform(st_setsrid(ST_Multi(ST_Union(st_point(-180,-58.077876),st_point(180,58.077876))),4326),st_srid(rast)))::text, 'BOX', '') , ',', '),('),' ',',')::box FROM rat_clip WHERE name = 'xxxx影像';

说明

Ganos的6.0版本提供Raster Box类型与Geometry Box2d类型之间的转换函数，上述嵌套的replace操作可以通过调用::box进行类型转换后简化。

SELECT st_extent(rast)::box2d::box FROM rat_clip WHERE name = 'xxxx影像';

代码实践

栅格数据的浏览通常可以通过Geoserver发布服务的方式进行，Ganos同样支持基于Geoserver的栅格和矢量服务发布，详情可参考地图服务。此处介绍一种更为简单的低代码方法，不需要依赖任何GIS Server工具，使用最少的Python语言快速搭建一个能随用户对地图的拖拽与缩放，动态更新显示我们指定影像数据的地图应用，更便于业务系统集成。

整体结构

后端代码

为实现代码的简洁性并更注重逻辑描述，选用Python作为后端语言，Web框架采用了基于Python的Flask框架，数据库连接框架则使用了基于Python的Psycopg2（可使用pip install psycopg2安装）。本案例在后端实现了一个简易的地图服务，自动建立金字塔，响应前端请求返回指定范围内的影像数据。将以下代码保存为Raster.py文件，执行python Raster.py命令即可启动服务。

## -*- coding: utf-8 -*-
## @File : Raster.py

import json
from flask import Flask, request, Response, send_from_directory
import binascii
import psycopg2

## 连接参数
CONNECTION =  "dbname=<database_name> user=<user_name> password=<user_password> host=<host> port=<port>"

## 影像地址
OSS_RASTER = "oss://<access_id>:<secrect_key>@<Endpoint>/<bucket>/path_to/file.tif"

## 分块表表名
RASTER_NAME = "xxxx影像"

## 分块表表名
CHUNK_TABLE = "chunk_table"

## 主表名
RASTER_TABLE = "raster_table"

## 字段名
RASTER_COLUMN = "rast"

## 默认渲染参数
DEFAULT_CONFIG = {
    "strength": "ratio",
    "quality": 70
}


class RasterViewer:
    def __init__(self):
        self.pg_connection = psycopg2.connect(CONNECTION)
        self.column_name = RASTER_COLUMN
        self.table_name = RASTER_TABLE
        self._make_table()
        self._import_raster(OSS_RASTER)

    def poll_query(self, query: str):
        pg_cursor = self.pg_connection.cursor()
        pg_cursor.execute(query)
        record = pg_cursor.fetchone()
        self.pg_connection.commit()
        pg_cursor.close()
        if record is not None:
            return record[0]

    def poll_command(self, query: str):
        pg_cursor = self.pg_connection.cursor()
        pg_cursor.execute(query)
        self.pg_connection.commit()
        pg_cursor.close()

    def _make_table(self):
        sql = f"create table if not exists {self.table_name} (ID INT PRIMARY KEY NOT NULL,name text, {self.column_name} raster);"
        self.poll_command(sql)

    def _import_raster(self, raster):
        sql = f"insert into {self.table_name} values (1, '{RASTER_NAME}', ST_ComputeStatistics(st_buildpyramid(ST_ImportFrom('{CHUNK_TABLE}','{raster}'),'{CHUNK_TABLE}'))) on conflict (id) do nothing;;"
        self.poll_command(sql)
        self.identify = f" name= '{RASTER_NAME}'"

    def get_extent(self) -> list:
        """获取影像范围"""
        import re
        sql = f"select replace((box2d(st_transform(st_envelope({self.column_name}),4326)))::text,'BOX','') from {self.table_name} where {self.identify}"
        result = self.poll_query(sql)

        # 转换为前端方便识别的形式
        bbox = [float(x) for x in re.split(
                '\(|,|\s|\)', result) if x != '']
        return bbox

    def get_jpeg(self, bbox: list, width: int, height: int) -> bytes:
        """
        获取指定位置的影像
        :param bbox: 指定位置的外包框
        :param width: 视区控件宽度
        :param height: 视区控件高度
        """

        # 指定波段和渲染参数
        bands = "0-2"
        options = json.dumps(DEFAULT_CONFIG)

        # 获取范围
        boxSQl = f"Replace(Replace(Replace(box2d(st_transform(st_setsrid(ST_Multi(ST_Union(st_point({bbox[0]},{bbox[1]}),st_point({bbox[2]},{bbox[3]}))),4326),st_srid({self.column_name})))::text, 'BOX', ''), ',', '),('),' ',',')::box"
        sql = f"select encode(ST_AsImage({self.column_name},{boxSQl} ,ST_BestPyramidLevel({self.column_name},{boxSQl},{width},{height}),'{bands}','jpeg','{options}'),'hex')  from {self.table_name} where {self.identify}"
        result = self.poll_query(sql)
        result = binascii.a2b_hex(result)
        return result


rasterViewer = RasterViewer()
app = Flask(__name__)


@app.route('/raster/image')
def raster_image():
    bbox = request.args['bbox'].split(',')
    width = int(request.args['width'])
    height = int(request.args['height'])
    return Response(
        response=rasterViewer.get_jpeg(bbox, width, height),
        mimetype="image/jpeg"
    )


@app.route('/raster/extent')
def raster_extent():
    return Response(
        response=json.dumps(rasterViewer.get_extent()),
        mimetype="application/json",
    )


@app.route('/raster')
def raster_demo():
    """代理前端页面"""
    return send_from_directory("./", "Raster.html")


if __name__ == "__main__":
    app.run(port=5000, threaded=True)

针对影响数据，后端代码主要实现以下三个功能：

初始化数据：包括创建测试表，导入数据，建立金字塔并统计。
获取影像位置：辅助前端地图快速定位影像，跳转到影像所在的位置。
提取影像为图片格式：
当前针对该影像，已经提前获取其元数据，了解该影像为3波段数据，对应也可以通过ST_NumBands函数动态获取其波段数。
根据前端传回的范围信息和可视区域控件的像素宽高确定其范围。
在使用psycopg2库时，以16进制传递数据的效率更高，如果使用其他框架，或者使用其他语言时，直接以二进制形式获取即可，无需编码转换。

在此使用Python作为示例代码，若使用其他语言，同时实现相同的逻辑，也能达到同样的效果。

前端代码

本案例选择Mapbox作为前端地图框架，同时引入Turf前端空间库，计算用户进行拖拽、缩放等地图操作后，当前视域与原始影像的交集，并向服务端请求该区域更清晰的图像，实现地图随操作更新影像的效果。在后端代码的同一文件目录下，新建一个名为Raster.html的文件，并写入以下代码。在后端服务启动后，就可以通过http://localhost:5000/raster访问。

<!DOCTYPE html>
<html>

<head>
  <meta charset="UTF-8" />
  <title></title>
  <link href="https://cdn.bootcdn.net/ajax/libs/mapbox-gl/1.13.0/mapbox-gl.min.css" rel="stylesheet" />
</head>
<script src="https://cdn.bootcdn.net/ajax/libs/mapbox-gl/1.13.0/mapbox-gl.min.js"></script>
<script src="https://cdn.bootcdn.net/ajax/libs/axios/0.21.0/axios.min.js"></script>
<script src="https://cdn.bootcdn.net/ajax/libs/lodash.js/4.17.20/lodash.min.js"></script>
<script src="https://cdn.bootcdn.net/ajax/libs/Turf.js/5.1.6/turf.min.js"></script>

<body>
  <div id="map" style="height: 100vh" />
  <script>

    // 初始化地图控件
    const map = new mapboxgl.Map({
      container: "map",
      style: { version: 8, layers: [], sources: {} },
    });

    class Extent {
      constructor(geojson) {
        // 默认使用Geojson格式
        this._extent = geojson;
      }

      static FromBBox(bbox) {
        // 从BBOx格式生成Extent对象
        return new Extent(turf.bboxPolygon(bbox));
      }

      static FromBounds(bounds) {
        // 从Mapbox的Bounds生成Extent对象
        const bbox = [
          bounds._sw.lng,
          bounds._sw.lat,
          bounds._ne.lng,
          bounds._ne.lat,
        ];
        return Extent.FromBBox(bbox);
      }

      intersect(another) {
        // 判断相交区域
        const inrersect = turf.intersect(this._extent, another._extent);
        return inrersect ? new Extent(inrersect) : null;
      }

      toQuery() {
        // 转换为query格式
        return turf.bbox(this._extent).join(",");
      }

      toBBox() {
        // 转换为BBox格式
        return turf.bbox(this._extent);
      }

      toMapboxCoordinates() {
        // 转换为Mapbox Coordinate格式
        const bbox = this.toBBox();
        const coordinates = [
          [bbox[0], bbox[3]],
          [bbox[2], bbox[3]],
          [bbox[2], bbox[1]],
          [bbox[0], bbox[1]],
        ];
        return coordinates;
      }
    }

    map.on("load", async () => {
      map.resize();

      const location = window.location.href;

      // 拼接查询语句
      const getUrl = (extent) => {
        const params = {
          bbox: extent.toQuery(),
          height: map.getCanvas().height,
          width: map.getCanvas().width,
        };
        // 拼接请求体
        const url = `${location}/image?${Object.keys(params)
          .map((key) => `${key}=${params[key]}`)
          .join("&")}`;
        return url;
      };

      // 查询影像范围
      const result = await axios.get(`${location}/extent`);
      const extent = Extent.FromBBox(result.data);
      const coordinates = extent.toMapboxCoordinates();

      // 添加数据源
      map.addSource("raster_source", {
        type: "image",
        url: getUrl(extent),
        coordinates,
      });

      //添加图层
      //使用了Mapbox的image类型图层，采用将图片附着在指定位置上的方法，让影像在地图上展示
      map.addLayer({
        id: "raster_layer",
        paint: { "raster-fade-duration": 300 },
        type: "raster",
        layout: { visibility: "visible" },
        source: "raster_source",
      });

      // 跳转到影像位置
      map.fitBounds(extent.toBBox());

      // 绑定刷新方法
      map.once("moveend", () => {
        const updateRaster = () => {
          const _extent = Extent.FromBounds(map.getBounds());
          let intersect;
          // 当可视区域没有图形时不重新请求
          if (!_extent || !(intersect = extent.intersect(_extent))) return;

          // 更新图形
          map.getSource("raster_source").updateImage({
            url: getUrl(intersect),
            coordinates: intersect.toMapboxCoordinates(),
          });
        };

        // 添加防抖,减少无效请求
        const _updateRaster = _.debounce(updateRaster, 200);
        map.on("zoomend", _updateRaster);
        map.on("moveend", _updateRaster);
      });
    });
  </script>
</body>

</html>

由于当前的影像获取接口不是一个标准地图协议接口，因此需要手动实现与影像更新相关的功能。需要解决的最大问题是：如何确定用户视域范围内的影像范围。解决的基本逻辑为：

获取影像的空间范围。
获取用户当前视域的空间范围。
获取两者的交集，即为预期的影像范围。

为在前端实现预定的空间判断功能，引入了Turf框架，将简单的计算在前端实现，减少不必要的请求次数。为了方便空间判断与各种格式间转换，我们还实现了一个辅助类Extent，包括如下功能：

格式互转：
- BBox <=> Geojson
- Mapbox Coordinate <=> Geojson
- Geojson => Query
- Bounds => Geojson
封装了Turf的空间判断方法。

效果展示

数据总览效果

在pgAdmin中集成影像浏览功能

支持将影像浏览功能集成在兼容PolarDB的数据库客户端软件pgAdmin中，即可快速浏览库内的影像数据，评估数据概况，增强数据管理的使用体验。

总结

利用Ganos Raster的相关函数可以实现提取库内影像数据的功能，最终通过少量代码实现一个可交互浏览库内遥感影像的地图应用。使用Ganos管理日益增长的遥感影像，不仅可以降低管理成本，还可以通过Ganos Raster扩展，在不借助第三方复杂工具的情况下，使用少量代码即可直接浏览库内影像数据，大大增强数据管理的体验。