1. 引言

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

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

不过，如果你使用云原生关系型数据库PolarDB（兼容PostgreSQL版本）管理影像，可以利用Ganos Raster扩展，无需借助第三方工具，仅利用SQL语句，即可快速将库内影像高效的展现在地图上。

本文将按照如下顺序展开，介绍如何利用Ganos实现影像浏览功能：

前期准备：将影像从OSS中导入PolarDB，并创建金字塔；
语句解析：了解如何从数据库中以图片格式提取遥感影像；
代码实践：使用Python语言编写一个简易的地图服务，配合Mapbox框架将数据库内的遥感影像展示在地图上。

2. 导入影像数据

在数据库内创建raster扩展。

CREATE EXTENSION ganos_raster CASCADE;

创建一张表raster_table，作为我们的测试表。

create table raster_table (ID INT PRIMARY KEY NOT NULL,name text,rast raster);

接下来我们从OSS中导入一幅影像作为测试数据。

我们使用Ganos的ST_ImportFrom方法导入影像到分块表chunk_table，更详细的参数描述可以参考官方文档。

insert into raster_table values (1, 'xxxx影像', ST_ImportFrom('chunk_table','OSS://ABCDEFG:1234567890@oss-cn.aliyuncs.com/mybucket/data/4.tif'));

3. 创建金字塔

金字塔是快速浏览影像数据的前提，对于新导入的数据，我们建议先创建金字塔。Ganos为我们提供了ST_BuildPyramid方法来创建金字塔。

方法原型如下，更详细的参数描述可以参考官方文档。

raster ST_BuildPyramid(raster source,  [cstring chunkTableName])

其中

source：需要创建金字塔的影像。
chunkTableName：金字塔所存储的分块表名称。
该方法返回更新后的raster对象。

实际调用如下

update raster_table set rast = st_buildpyramid(raster_table,'chunk_table') where name = 'xxxx影像';

我们为导入的测试数据创建了金字塔，存储于分块表chunk_table中，并使用Update方法更新了raster对象。

4. 从数据库内提取影像

创建好金字塔，我们就可以使用Ganos的ST_AsImage方法从数据库中获取指定范围的影像。

方法原型如下，更详细的参数描述可以参考官方文档。

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

4.1 静态参数解析

静态参数指一般不随用户操作而变化的参数，我们可以在代码中固定它们：

bands：指定要显示的波段。
待显示的波段列表，可以使用0-3或0,1,2,3等形式，不能超过影像本身的波段。
format：格式。
默认为PNG格式，可选JPEG格式。由于PNG格式对数据的压缩效果不如JPEG格式的有损压缩，因此在传输时将耗费更多的时间，在没有透明度需求的情况下应尽量选择JPEG格式。
option：其他参数。
JSON字符串格式。可以定义额外的渲染参数。

4.2 动态参数解析

动态参数随用户操作而变化，需要动态生成它们：

extent：需要获取的影像范围。
在条件一样的情况下，显示范围越大，数据库处理时间越长，返回图片的体积也越大，总体响应时间越长。
因此最好只获取用户视域范围内的影像，以保证传输的效率。
pyramidLevel：使用的金字塔级别。
使用的金字塔级别越高，图片越清晰，体积也就越大。
因此需要选择最适宜的金字塔级别，以保证传输的效率。

4.2.1 获取影像外包框范围

使用Ganos的ST_Envelope方法可以获取到影像的外包框范围，并且使用ST_Transform方法将影像转换到常用的WGS 84坐标系下，最终将其转换为前端方便使用的格式。

实际调用如下

select replace((box2d(st_transform(st_envelope(geom),4326)))::text,'BOX','') from rat_clip where name = 'xxxx影像'

4.2.2 获取最适宜金字塔级别

使用Ganos的ST_BestPyramidLevel方法可以计算特定影像范围内最适宜的金字塔级别。

方法原型如下，更详细的参数描述可以参考官方文档。

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

其中

extent：可视区域所要获取的影像范围。
与ST_AsImage方法中使用的范围一致。
width/height：可视区域的像素宽高。
一般来说，这就是我们使用的前端地图框的尺寸。

值得一提的是，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影像';
-- 返回 (180.0,58.077876),(-180.0,-58.077876)

注：Ganos在6.0版本（计划2023年11月下旬上线发布）会提供Raster Box类型与Geometry Box2d类型之间的转换函数，上述嵌套的replace操作可以通过调用::box进行类型转换后简化

5. 基于Ganos栅格能力的最佳实践

栅格数据的浏览一般可采用Geoserver发布服务的方式，Ganos同样也支持基于Geoserver的栅格/矢量服务发布，详情可参考官方文档。本文介绍了一种更为简单的低代码方法，不需要依赖任何GIS Server工具，用最少的代码快速搭建一个能随用户对地图的拖拽与缩放，动态更新显示我们指定影像数据的地图应用，这种方式更便于业务系统集成，也方便用户理解Ganos栅格的相关能力。

5.1 整体结构

5.2 后端代码

为了代码简洁，更侧重于逻辑的描述，我们选择了Python作为后端语言，Web框架使用了基于Python的Flask框架，数据库连接框架使用了基于Python的Psycopg2（使用pip install psycopg2-binary进行安装）。

我们在后端实现了一个简易的地图服务，自动建立金字塔，响应前端请求返回指定范围内的影像数据。

# -*- coding: utf-8 -*-
# @File : Raster.py
import json
from flask import Flask, request, Response, send_from_directory
import binascii
import psycopg2
# 连接参数
CONNECTION = "dbname=postgres user=postgres password=postgres host=YOUR_HOST port=5432"
# 影像地址
OSS_RASTER = "OSS://ABCDEFG:1234567890@oss-cn.aliyuncs.com/mybucket/data/4.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)

将以上代码保存为Raster.py文件，执行python Raster.py命令即可启动服务。

5.3 后端代码解析

可以看出，针对影像数据，后端代码主要实现了3个功能：

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

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

5.4 前端代码

前端我们采用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>

5.5 前端代码解析

由于我们的影像获取接口不是一个标准地图协议接口，因此需要手动实现与影像更新相关的功能。

我们要解决的最大问题是：如何确定用户视域范围内的影像范围。

解决的基本逻辑为：

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

为了在前端实现预定的空间判断功能，我们引入了Turf框架，将简单的计算在前端实现，减少不必要的请求次数。

为了方便空间判断与各种格式间转换，我们还实现了一个辅助类Extent，包括如下功能：

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

6. 效果

6.1 数据总览效果

6.2 在PGAdmin中集成影像浏览功能

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

7. 总结

本文介绍了如何利用Ganos Raster的相关方法实现提取库内影像数据的功能，并最终通过少量代码实现了一个可交互浏览库内遥感影像的地图应用。

从本实践中我们可以看到，使用Ganos管理日益增长的遥感影像，不仅可以降低管理成本，还可以通过Ganos Raster扩展，不借助第三方复杂工具，使用少量代码即可直接浏览库内影像数据，大大增强了数据管理的体验。

用Ganos低代码实现免切片遥感影像浏览（一）

1. 引言

2. 导入影像数据

3. 创建金字塔

4. 从数据库内提取影像

4.1 静态参数解析

4.2 动态参数解析

4.2.1 获取影像外包框范围

4.2.2 获取最适宜金字塔级别

5. 基于Ganos栅格能力的最佳实践

5.1 整体结构

5.2 后端代码

5.3 后端代码解析

5.4 前端代码

5.5 前端代码解析

6. 效果

6.1 数据总览效果

6.2 在PGAdmin中集成影像浏览功能

7. 总结