名前

ST_MapAlgebraExpr — 2バンド版: 二つの入力バンドに対する妥当な PostgreSQL 代数演算で形成された、指定したピクセルタイプとなる 1バンドラスターを生成します。バンドを指定しない場合には、どちらも 1番と仮定します。結果ラスターは、一つ目のラスターのアラインメント (スケール、スキュー、ピクセル角位置) にあわされます。範囲は"extenttype"引数で定義されます。取りうる"extenttype"の値は INTERSECTION, UNION, FIRST, SECOND です。

概要

raster ST_MapAlgebraExpr(raster rast1, raster rast2, text expression, text pixeltype=same_as_rast1_band, text extenttype=INTERSECTION, text nodata1expr=NULL, text nodata2expr=NULL, double precision nodatanodataval=NULL);

raster ST_MapAlgebraExpr(raster rast1, integer band1, raster rast2, integer band2, text expression, text pixeltype=same_as_rast1_band, text extenttype=INTERSECTION, text nodata1expr=NULL, text nodata2expr=NULL, double precision nodatanodataval=NULL);

説明

[警告]

ST_MapAlgebraExpr は 2.1.0で非推奨になりました。代わりにST_MapAlgebra (expression version) を使います。

expressionで定義された妥当な二つのバンドへの PostgreSQL 代数演算を入力ラスター rast1, (rast2) に適用して、一つのバンドを持つラスターを生成します。band1, band2が指定されない場合には、1番バンドと仮定します。新しいラスターは、一つ目のラスターと同じアラインメント (スケール、スキュー、ピクセル隅) を持ちます。新しいラスターは、extenttype引数で定義される範囲になります。

expression

二つのラスターと PostgreSQL 定義済み関数/演算子を含む PostgreSQL 代数式です。関数と演算子は、二つのピクセルがインタセクトするピクセルの値を定めます。たとえば(([rast1] + [rast2])/2.0)::integer といったふうになります。

pixeltype

出力ラスターのピクセルタイプです。必ずST_BandPixelTypeに挙げられたものの一つになるか、省略されるか、NULL に設定されます。引数として渡されないか NULL が渡された場合には、一つ目のラスターのピクセルタイプになります。

extenttype

新しいラスターの範囲を制御します。

  1. INTERSECTION - 新しいラスターの範囲は二つのラスターのインタセクトした領域です。これがデフォルトです。

  2. UNION - 新しいラスターの範囲は二つのラスターの結合です。

  3. FIRST - 新しいラスターの範囲は一つ目のラスターと同じです。

  4. SECOND - 新しいラスターの範囲は二つ目のラスターと同じです。

nodata1expr

rast1が NODATA 値で、特に rast2ピクセルに値がある時に、rast2だけを返すか返すべき値を定義する定数を含む代数式です。

nodata2expr

rast2が NODATA 値で、特に rast2ピクセルに値がある時に、rast1だけを返すか返すべき値を定義する定数を含む代数式です。

nodatanodataval

rast1と rast2のピクセルの両方が NOADTA 値になる場合に返すべき定数です。

pixeltypeが渡された場合には、新しいラスターは、指定されたピクセルタイプのバンドを持ちます。pixeltype として NULL が渡されたりピクセルタイプを指定しない場合には、新しいラスターはrast1と同じピクセルタイプになります。

数式の中で使える語は、元バンドのピクセル値を参照する [rast1.val], [rast2.val]、1始まりの列/行インデックスを参照する[rast1.x], [rast1.y]などです。

Availability: 2.0.0

2バンドのインターセクトしている部分と結合。

元のラスターから 1バンドラスターを生成します。元のラスターバンドの値について 2で割った余りが入ります。

ラスターのいかした集合を生成します。

メートル単位マサチューセッツ州平面でボストン周囲にイかした形状を挿入して、インターセクト部分と結合部分とを示します。

コード
DROP TABLE IF EXISTS fun_shapes;
CREATE TABLE fun_shapes(rid serial PRIMARY KEY, fun_name text, rast raster);

INSERT INTO fun_shapes(fun_name, rast)
VALUES ('ref', ST_AsRaster(ST_MakeEnvelope(235229, 899970, 237229, 901930, 26986), 200, 200, '8BUI', 0, 0));

INSERT INTO fun_shapes(fun_name, rast)
WITH ref(rast) AS (SELECT rast FROM fun_shapes WHERE fun_name = 'ref' )
SELECT 'area' AS fun_name, ST_AsRaster(ST_Buffer(ST_SetSRID(ST_Point(236229, 900930), 26986), 1000),
            ref.rast, '8BUI', 10, 0) As rast
FROM ref
UNION ALL
SELECT 'rand bubbles',
       ST_AsRaster(
           (
               SELECT ST_Collect(geom)
               FROM (
                   SELECT ST_Buffer(
                       ST_SetSRID(ST_Point(236229 + i*random()*100, 900930 + j*random()*100), 26986),
                       random()*20
                   ) AS geom
                   FROM generate_series(1, 10) AS i,
                        generate_series(1, 10) AS j
               ) AS foo
           ),
           ref.rast, '8BUI', 200, 0
       )
FROM ref;

SELECT  ST_MapAlgebraExpr(area.rast, bub.rast, '[rast2.val]', '8BUI', 'INTERSECTION', '[rast2.val]', '[rast1.val]') As interrast,
        ST_MapAlgebraExpr(area.rast, bub.rast, '[rast2.val]', '8BUI', 'UNION', '[rast2.val]', '[rast1.val]') As unionrast
FROM
  (SELECT rast FROM fun_shapes WHERE fun_name = 'area') As area
CROSS JOIN
  (SELECT rast FROM fun_shapes WHERE fun_name = 'rand bubbles') As bub
                    

二つの決定的かつ自己完結型のラスターを使た同じインターセクト部分と結合部分の計算。

コード
WITH canvas AS (
    SELECT ST_AddBand(
        ST_MakeEmptyRaster(80, 80, 0, 80, 1, -1, 0, 0, 0),
        1, '8BUI', 40, 0
    ) AS rast
), inputs AS (
    SELECT
        rast AS raster_1,
        ST_AsRaster(
            ST_Buffer(ST_Point(40, 40), 28),
            rast,
            '8BUI', 220, 0
        ) AS raster_2
    FROM canvas
), variants AS (
    SELECT 'raster 1' AS title, raster_1 AS rendered FROM inputs
    UNION ALL
    SELECT 'raster 2', raster_2 FROM inputs
    UNION ALL
    SELECT 'intersection', ST_MapAlgebraExpr(
        raster_1, raster_2,
        '[rast2.val]', '8BUI', 'INTERSECTION',
        '[rast2.val]', '[rast1.val]'
    ) FROM inputs
    UNION ALL
    SELECT 'union', ST_MapAlgebraExpr(
        raster_1, raster_2,
        '[rast2.val]', '8BUI', 'UNION',
        '[rast2.val]', '[rast1.val]'
    ) FROM inputs
)
SELECT title, ST_AsPNG(rendered) AS image
FROM variants
ORDER BY CASE title
    WHEN 'raster 1' THEN 1
    WHEN 'raster 2' THEN 2
    WHEN 'intersection' THEN 3
    ELSE 4
END;
出力
title        | image
--------------+---------------------------
 raster 1     | PNG image, 80 x 80 pixels
 raster 2     | PNG image, 56 x 56 pixels
 intersection | PNG image, 56 x 56 pixels
 union        | PNG image, 80 x 80 pixels
Geometry figure for visual-rt-st-mapalgebraexpr2-01

別のバンドとして 1単位 1ピクセルのキャンバスの上にオーバーレイしているラスター。ST_AsPNG を使って画像を描画しています。単一バンドのラスターは全ての灰色に見えます。

コード
WITH mygeoms AS (
    SELECT 2 AS bnum, ST_Buffer(ST_Point(1, 5), 10) AS geom
    UNION ALL
    SELECT 3,
        ST_Buffer(
            ST_GeomFromText('LINESTRING(50 50,150 150,150 50)'),
            10,
            'join=bevel'
        )
    UNION ALL
    SELECT 1,
        ST_Buffer(
            ST_GeomFromText('LINESTRING(60 50,150 150,150 50)'),
            5,
            'join=bevel'
        )
), canvas AS (
    SELECT ST_AddBand(
        ST_MakeEmptyRaster(
            200,
            200,
            ST_XMin(e)::integer,
            ST_YMax(e)::integer,
            1, -1, 0, 0),
        '8BUI'::text,
        0
    ) AS rast
    FROM (
      SELECT ST_Extent(geom) AS e
      FROM mygeoms
    ) AS bounds
), rbands AS (
    SELECT ARRAY(
        SELECT ST_MapAlgebraExpr(
            canvas.rast,
            ST_AsRaster(m.geom, canvas.rast, '8BUI', 100),
            '[rast2.val]', '8BUI', 'FIRST',
            '[rast2.val]', '[rast1.val]'
        )
        FROM mygeoms AS m
        CROSS JOIN canvas
        ORDER BY m.bnum
    ) AS rasts
)
SELECT title, ST_AsPNG(rendered) AS image
FROM rbands
CROSS JOIN LATERAL (VALUES
    ('band 1', rasts[1]),
    ('band 2', rasts[2]),
    ('band 3', rasts[3]),
    ('RGB', ST_AddBand(ST_AddBand(rasts[1], rasts[2]), rasts[3]))
) AS variants(title, rendered);
出力
title  | image
--------+-----------------------------
 band 1 | PNG image, 200 x 200 pixels
 band 2 | PNG image, 200 x 200 pixels
 band 3 | PNG image, 200 x 200 pixels
 RGB    | PNG image, 200 x 200 pixels
Geometry figure for visual-rt-st-mapalgebraexpr2-02

空中写真上に、選択した区画の 2メートル幅の境界線を重ねます。クエリーでは最初に区画ジオメトリーを結合して、元ラスターの破片を領域で切り抜き、小さい破片を結合し、最後に、バンドの先頭 2バンドと、地図代数で描いた第 3バンドを追加します。

新しい 3バンドのラスターを生成します。最初の二つは切り抜いたバンドで、三つ目のバンドはジオメトリーから生成します。クエリーは Windows 用 PostGIS で 3.6秒かかりました。

コード
WITH pr AS (
    SELECT ST_Clip(rast, ST_Expand(geom, 50)) AS rast, g.geom
    FROM aerials.o_2_boston AS r
    INNER JOIN (
        SELECT ST_Union(ST_Transform(geom, 26986)) AS geom
        FROM landparcels
        WHERE pid IN ('0303890000', '0303900000')
    ) AS g
        ON ST_Intersects(rast::geometry, ST_Expand(g.geom, 50))
), prunion AS (
    SELECT ST_AddBand(
        NULL,
        ARRAY[
            ST_Union(rast, 1),
            ST_Union(rast, 2),
            ST_Union(rast, 3)
        ]
    ) AS clipped,
    geom
    FROM pr
    GROUP BY geom
)
SELECT ST_AddBand(
    ST_Band(clipped, ARRAY[1, 2]),
    ST_MapAlgebraExpr(
     ST_Band(clipped, 3),
     ST_AsRaster(ST_Buffer(ST_Boundary(geom), 2), clipped, '8BUI', 250),
     '[rast2.val]',
     '8BUI',
     'FIRST',
     '[rast2.val]',
     '[rast1.val]') ) As rast
FROM prunion;
                    

区画境界による上書き。

青線が選択した区画の境界です。

上の生成された 4個パネルのキャンバスの例では、外部空中写真や区画テーブルを求めることなしに同じバンド重ね合わせ技術を示しています。