説明

2つの入力ラスタ(rast1, rast2)に対してexpressionで定義される妥当なPostgreSQL代数演算で形成されラスタを返します。 band1とband2が指定されていない場合は1番と仮定します。結果ラスタは、ひとつめのラスタのアラインメント(スケール、スキュー、ピクセル角位置)にあわされます。範囲は"extenttype"引数で定義されます。

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

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

初出: 2.0.0

例: 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
				

INTERSECTIONの場合

UNIONの場合

例: 別個のバンドとしてのキャンバス上へのラスタのオーバレイ

-- ST_AsPNGを使って1バンドで灰色に見える画像を作成します --
WITH mygeoms 
    AS ( SELECT 2 As bnum, ST_Buffer(ST_Point(1,5),10) As geom
            UNION ALL
            SELECT 3 AS bnum, 
                ST_Buffer(ST_GeomFromText('LINESTRING(50 50,150 150,150 50)'), 10,'join=bevel') As geom
            UNION ALL
            SELECT 1 As bnum, 
                ST_Buffer(ST_GeomFromText('LINESTRING(60 50,150 150,150 50)'), 5,'join=bevel') As geom
            ),
   -- define our canvas to be 1 to 1 pixel to geometry
   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,
                    Max(ST_SRID(geom)) As srid 
                    from mygeoms 
                    ) As foo
            ),
   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]') As rast
                FROM mygeoms AS m CROSS JOIN canvas
                ORDER BY m.bnum) As rasts
                )
          SELECT rasts[1] As rast1 , rasts[2] As rast2, rasts[3] As rast3, ST_AddBand(
                    ST_AddBand(rasts[1],rasts[2]), rasts[3]) As final_rast
            FROM rbands;

rast1	rast2
rast3	final_rast

例: 選択した区画の2メートル幅境界線の空中写真へのオーバレイ

-- 前2つの切り取ったバンドと
--   第3バンドにジオメトリをオーバレイしたバンドの
--   3バンドのラスタを生成します
-- このクエリはWindows 64ビット版の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
-- あとでインタセクトできるように対象区画を結合して1つのジオメトリにします
		(SELECT ST_Union(ST_Transform(the_geom,26986)) AS geom 
		  FROM landparcels WHERE pid IN('0303890000', '0303900000')) As g
		ON ST_Intersects(rast::geometry, ST_Expand(g.geom,50))
),
-- ラスタの破片を結合します
-- ラスタのST_Unionは非常に遅いですがラスタを小さくすれば速度が上がります
-- 先に切り出しておいてから結合しているのはこのためです
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)
-- ラスタの破片が結合され、区画境界線のオーバレイも持った
-- 最終のラスタを返します
-- 最初の2バンドと3番バンド+ジオメトリの地図代数関数とを追加します
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;

The blue lines are the boundaries of select parcels

関連情報

ST_MapAlgebraExpr, ST_AddBand, ST_AsPNG, ST_AsRaster, ST_MapAlgebraFct, ST_BandPixelType, ST_GeoReference, ST_Value , ST_Union, ST_Union