GanosBase低代碼實現免切片遙感影像瀏覽（一）：金字塔 - PolarDB

使用PolarDB PostgreSQL版資料庫，配合GanosBase時空資料庫引擎，在不藉助第三方工具的情況下，可僅使用SQL語句快速管理與展示遙感影像資料。GanosBase共提供兩種影像免切瀏覽的方法，一種藉助GanosBase Raster外掛程式，使用視窗範圍擷取影像資料展示，另一種通過固定瓦片範圍擷取影像資料展示，本文詳細介紹第一種方法並提供前後端實操代碼協助您快速理解GanosBase Raster的使用細節。

背景

當管理著日益增長的遙感資料時，想要將任一影像展現在地圖上往往要思考如下問題：

對於可能只使用幾次的影像進行傳統的切片、發布流程是否興師動眾？瓦片儲存在何處？瓦片資料如何進行後續的管理？
不進行預切片而採用即時切片的方案能否滿足前端的響應需求？遇到特別巨大的影像時，即時切片會受到多大影響？

不過，使用PolarDB PostgreSQL版管理影像時，藉助GanosBase 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'));
```

建立金字塔

金字塔是快速探索影像資料的基礎。對於新匯入的資料，建議首先建立金字塔。GanosBase提供ST_BuildPyramid函數以實現金字塔的建立。詳細的文法介紹請參考ST_BuildPyramid。

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

語句解析

建立金字塔成功後，可以使用GanosBase的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：使用的影像金字塔層級。
使用的金字塔層級越高，映像的清晰度越高，但其體積也隨之增加。因此，需要選擇最合適的金字塔層級，以確保傳輸效率。

擷取影像外包框範圍

使用GanosBase的ST_Envelope函數擷取影像的外包框範圍，然後使用ST_Transform函數將影像轉換到常用座標系（此處使用WGS 84座標系），最終將其轉換為前端方便使用的格式。

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

擷取最適宜金字塔層級

使用GanosBase的ST_BestPyramidLevelFunction Compute特定影像範圍內最適宜的金字塔層級。以下為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發布服務的方式進行，GanosBase同樣支援基於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中，即可快速探索庫內的影像資料，評估資料概況，增強資料管理的使用體驗。

總結

利用GanosBase Raster的相關函數可以實現提取庫內影像資料的功能，最終通過少量代碼實現一個可互動瀏覽庫內遙感影像的地圖應用。使用GanosBase管理日益增長的遙感影像，不僅可以降低管理成本，還可以通過Ganos Raster擴充，在不藉助第三方複雜工具的情況下，使用少量代碼即可直接瀏覽庫內影像資料，大大增強資料管理的體驗。