😀

深い深い穴を探せ

2023/01/13に公開

postgis

GIS

tech

はじめに

基盤地図情報DEM配信サービスを公開した記念の記事 DEMのCOGを始めましたで「-134mってどこかはともかく」と申しました。

これ、しれっと言っていますが、データ変換とかやってるときにおかしくしたんじゃないかと思い、これはやばいぞと思いました。

とにかく原因究明と、自分のせいでおかしくしてたのならデータの作り直しになります。一刻も早く究明せねばなりません。

さっそく見ていきましょう。

最小値を出す関数を作ってみる

まずはファイルを特定しないといけません。方針としては、各行（基盤地図情報DEMのGMLファイルの1ファイルに相当）ごとに最小値を出して、-100mより下になるような行を抽出すると、それで概ね位置が特定できます。

そのうえで、ファイル内の最小値を取るセル位置も特定できるようにしておくと、さらにいいだろうと考えました。

で、作ってみたのが次の関数です。

CREATE TYPE rst_minval_result AS (
  x INTEGER,
  y INTEGER,
  v DOUBLE PRECISION
);


CREATE OR REPLACE FUNCTION rst_minvalue(
  rast RASTER,
  band INTEGER DEFAULT 1) RETURNS rst_minval_result AS $$
DECLARE
  x INTEGER;
  y INTEGER;
  nw INTEGER;
  nh INTEGER;
  arr DOUBLE PRECISION[][];
  ret rst_minval_result;
  v DOUBLE PRECISION;
BEGIN
  nw := ST_Width(rast);
  nh := ST_Width(rast);
  arr := ST_DumpValues(rast, band);
  FOR y IN 1..nh LOOP
    FOR x IN 1..nw LOOP
      v := arr[y][x];
      IF NOT v is NULL THEN
        IF ret IS NULL THEN
          ret := (x, y, v);
        ELSE
          IF v < ret.v THEN
            ret := (x, y, v);
          END IF;
        END IF;
      END IF;
    END LOOP;
  END LOOP;
  RETURN ret;
END;
$$ LANGUAGE plpgsql;

ST_DumpValues()でダンプしたものを二重ループで最小値を取るだけのものです。

なお、ST_Value()で値を取ろうとすると、とんでもなく遅くなりました。2次配列にダンプさせた方がいいです。

探索開始

上記の関数を dem10b の入ったテーブルで実行してみました。

db=# SELECT filename FROM dem10b WHERE (rst_minvalue(rast)).v < -100;
              filename              
------------------------------------
 FG-GML-6041-54-dem10b-20161001.xml
(1 □s)

どうもファイルはひとつだけのよう。

db=# SELECT rst_minvalue(rast) FROM dem10b WHERE filename = 'FG-GML-6041-54-dem10b-20161001.xml';
         rst_minvalue          
-------------------------------
 (334,432,-134.60000610351562)
(1 □s)

最小値を取るセルは X=334, Y=432 にあるものらしいです。

経度緯度はというと、

db=# SELECT ST_AsEWKT(ST_PixelAsPoint(rast, 334, 432))
  FROM dem10b
  WHERE filename = 'FG-GML-6041-54-dem10b-20161001.xml';
            st_asewkt
----------------------------------
 POINT(141.537 40.45211111130267)
(1 row)

東経141.537度、北緯40.45217度のようです。