线性插值

来自 PostgreSQL wiki
跳转到导航跳转到搜索

片段

线性插值

适用于 PostgreSQL

9.0

用

SQL

依赖于

无

PostgreSQL 中的线性插值

以下代码在 PostgreSQL 中基于包含要插值曲线的 x 和 y 值的数组（或 PostGIS 轨迹插值的 LINESTRING 几何图形）执行简单的线性插值。特别是，我用它来线性插值在特定时间移动对象的轨迹。代码很简单，但效率相当高且纯 PL/pgSQL，因此它应该在大多数安装上运行（执行地理空间计算的函数需要 PostGIS 扩展版本 >=2.0）。在 PostgreSQL 的较新版本中，某些代码是冗余的（例如，现在可以使用基本 generate_series 函数生成一系列浮点数或时间戳），但它仍然有效，并且如果需要，应该很容易理解如何现代化代码。

--- Simple function to do matlab style Xmin:DX:Xmax
CREATE OR REPLACE FUNCTION generate_float_series(Xmin float,Xmax float,DX float)
       RETURNS SETOF float AS
$BODY$
SELECT $1 + $3 * generate_series(0,floor(($2-$1)/$3)::int);
$BODY$ LANGUAGE 'sql' IMMUTABLE;

-------------------------------------------------------------------
--- Main linear interpolation function
-------------------------------------------------------------------
--- Last argument (na) determines if values outside limits of x values
--- generates an error or not.  If na is TRUE, then NULL will be
--- returned for values outside range.  If na is FALSE, then an error
--- will be generated.

--- EXAMPLES:
--- SELECT LinearInterpolation( ARRAY[1.,2,3,4], ARRAY[1.,4,9,20.5], ARRAY[-1.,2.5,2,25], FALSE );
--- SELECT LinearInterpolation( ARRAY[1.,2,3,4], ARRAY[1.,4,9,20.5], ARRAY[-1.,2.5,2,25], TRUE );

CREATE OR REPLACE FUNCTION LinearInterpolation(x float[], y float[], x2 float[], na boolean) 
       RETURNS float[] AS
$BODY$
DECLARE
	c1 SCROLL CURSOR FOR SELECT unnest(x) AS xi, unnest(y) AS yi 
	   	  	     	    ORDER BY xi ASC;
	c2 CURSOR FOR SELECT generate_subscripts(x2,1) AS i2, unnest(x2) AS xi2 
	   	      	     ORDER BY xi2 ASC;
	xi float;
	yi float;
	xj float;
	yj float;
	xi2 float;
	p float;
	i2 int;
	y2 float[] := array_fill(NULL::float, ARRAY[array_length(x2,1)],ARRAY[array_lower(x2,1)]);
	myex text=NULL::text;
BEGIN
	IF array_length(x,1) != array_length(y,1) THEN
	   RAISE EXCEPTION 'x,y input arrays not of same lengths';
	END IF;

	OPEN c1;
	OPEN c2;

	<<curs_block>>
	BEGIN

	FETCH c1 INTO xj,yj;
	FETCH c2 INTO i2,xi2;

	IF xj IS NULL THEN
	   myex := 'No input data';
	   EXIT curs_block;
	END IF;

	IF xi2 IS NULL THEN
	   myex := 'No interpolate data';
	   EXIT curs_block;
	END IF;

	WHILE xi2 IS NOT NULL AND xi2 < xj LOOP
	      FETCH c2 INTO i2,xi2;
	END LOOP;

	IF xi2 IS NULL THEN
	   EXIT curs_block;	   
	END IF;

	IF xi2 < xj AND NOT na THEN
	   myex := 'Interpolate data outside input data range';
	   EXIT curs_block;
	END IF;

	xi := xj;
	yi := yj;
	FETCH c1 INTO xj,yj;
	
	<<loop1>>
	WHILE xj IS NOT NULL LOOP

	      <<loop2>>
	      WHILE xi2 < xj LOOP      
	      	    p := (xi2-xi)/(xj-xi);
	      	    y2[i2] := (1.0-p) * yi + p * yj;

		    FETCH c2 INTO i2,xi2;

		    IF xi2 IS NULL THEN
		       EXIT curs_block;
		    END IF;
	      END LOOP;

	      xi := xj;
	      yi := yj;
	      FETCH c1 INTO xj,yj;
	END LOOP;

	WHILE xi2 IS NOT NULL AND xi2 = xi LOOP
	      y2[i2] := yi;
	      FETCH c2 INTO i2,xi2;
	END LOOP;

	IF xi2 IS NOT NULL AND NOT na THEN
	   myex := 'Interpolate data outside input data range';
	   EXIT curs_block;
	END IF;

	END; --curs_block

	CLOSE c1;
	CLOSE c2;

	IF myex IS NOT NULL THEN
	   RAISE EXCEPTION '%', myex;
	ELSE
	   RETURN y2;
	END IF;
END;	   
$BODY$ LANGUAGE 'plpgsql' IMMUTABLE;

--- Same as above, but defines the default value for na to TRUE for float arrays
CREATE OR REPLACE FUNCTION LinearInterpolation(
       x float[], y float[], x2 float[]) 
       RETURNS float[] AS
$BODY$
SELECT LinearInterpolation( $1, $2, $3, TRUE );
$BODY$ LANGUAGE 'sql' IMMUTABLE;

-------------------------------------------------------------------
--- Linear interpolation function for PostGIS LINESTRING geometries
-------------------------------------------------------------------

--- EXAMPLES:
--- SELECT LinearInterpolation( ARRAY[1.,2,3,4], ST_MakeLine(ARRAY[ST_MakePoint(1,1),ST_MakePoint(2,4),ST_MakePoint(3,9),ST_MakePoint(4,16)]), ARRAY[-1.,2.5,2,25], FALSE );
--- SELECT LinearInterpolation( ARRAY[1.,2,3,4], ST_MakeLine(ARRAY[ST_MakePoint(1,1),ST_MakePoint(2,4),ST_MakePoint(3,9),ST_MakePoint(4,16)]), ARRAY[-1.,2.5,2,25], TRUE );
--- SELECT LinearInterpolation( ARRAY[1.,2,3,4], ST_MakeLine(ARRAY[ST_MakePoint(1,1),ST_MakePoint(2,4),ST_MakePoint(3,9),ST_MakePoint(4,16)]), ARRAY[1.,2.5,2,4], TRUE );

CREATE OR REPLACE FUNCTION LinearInterpolation(x float[], g geometry, x2 float[],
       na boolean) 
       RETURNS geometry[] AS
$BODY$
DECLARE
	c1 SCROLL CURSOR FOR SELECT unnest(x) AS xi, 
	   	  	     	    (ST_DumpPoints(g)).geom AS gi;
	c2 CURSOR FOR SELECT generate_subscripts(x2,1) AS i2, unnest(x2) AS xi2 
	   	      	     ORDER BY xi2 ASC;
	xi float;
	gi geometry;
	xj float;
	gj geometry;
	xi2 float;
	p float;
	i2 int;
	g2 geometry[] := array_fill(ST_SetSRID('POINT EMPTY'::geometry,ST_SRID(g)), ARRAY[array_length(x2,1)],ARRAY[array_lower(x2,1)]);
	myex text=NULL::text;
	myexout text='Interpolate data outside input data range';
BEGIN
	IF array_length(x,1) != ST_NumPoints(g) THEN
	   RAISE EXCEPTION 'Number x not equal to number of points in g';
	END IF;

	OPEN c1;
	OPEN c2;

	<<curs_block>>
	BEGIN

	FETCH c1 INTO xj,gj;

	IF xj IS NULL THEN
	   myex := 'No input data';
	   EXIT curs_block;
	END IF;

	FETCH c2 INTO i2,xi2;
	IF xi2 IS NULL THEN
	   myex := 'No interpolate data';
	   EXIT curs_block;
	END IF;

	-- Look for points outside data range
	WHILE xi2 < xj LOOP
	   myex := myexout;
	   FETCH c2 INTO i2,xi2;
	END LOOP;

	xi := xj;
	gi := gj;
	FETCH c1 INTO xj,gj;
	
	<<loop1>>
	WHILE xj IS NOT NULL LOOP
	      IF xj < xi THEN
	      	 myex := 'Input x not increasing';
	   	 EXIT curs_block;
	      END IF;
	
	      <<loop2>>
	      WHILE xi2 < xj LOOP      
	      	    p := (xi2-xi)/(xj-xi);
	      	    g2[i2] := ST_SetSRID( 
		    	   ST_MakePoint((1.0-p) * ST_X(gi) + p * ST_X(gj),
		    	             	(1.0-p) * ST_Y(gi) + p * ST_Y(gj)),
				ST_SRID(g));

		    FETCH c2 INTO i2,xi2;

		    IF xi2 IS NULL THEN
		       EXIT curs_block;
		    END IF;
	      END LOOP;

	      xi := xj;
	      gi := gj;
	      FETCH c1 INTO xj,gj;
	END LOOP;

	WHILE xi2 IS NOT NULL AND xi2 = xi LOOP
	      g2[i2] := gi;

	      FETCH c2 INTO i2,xi2;
	END LOOP;

	-- Look for points outside data range
	IF xi2 IS NOT NULL THEN
	   myex := myexout;
	   EXIT curs_block;
	END IF;

	END; -- curs_block

	CLOSE c1;
	CLOSE c2;

	IF myex = myexout AND na THEN
	   -- If you data outside range and na=TRUE, return geo collection
	   RAISE WARNING '%', myex;
	   RETURN g2;
	ELSIF myex IS NOT NULL THEN
	   -- Else if any error string, raise exception
	   RAISE EXCEPTION '%', myex;
	ELSE
	   RETURN g2;
	END IF;
END;	   
$BODY$ LANGUAGE 'plpgsql' IMMUTABLE;

--- Same as above, but defines the default value for na to FALSE for geometries
--- Also converts geometry array to linestring

--- EXAMPLES:
--- SELECT LinearInterpolation( ARRAY[1.,2,3,4], ST_MakeLine(ARRAY[ST_MakePoint(1,1),ST_MakePoint(2,4),ST_MakePoint(3,9),ST_MakePoint(4,16)]), ARRAY[1.,2.5,2,4]);
CREATE OR REPLACE FUNCTION LinearInterpolation(
       x float[], y geometry, x2 float[]) 
       RETURNS geometry AS
$BODY$
SELECT ST_MakeLine(LinearInterpolation( $1, $2, $3, FALSE ));
$BODY$ LANGUAGE 'sql' IMMUTABLE;

--- A function that will generate series of timestamps
CREATE OR REPLACE FUNCTION generate_timestamp_series(
       Tmin TIMESTAMP WITHOUT TIME ZONE,
       Tmax TIMESTAMP WITHOUT TIME ZONE,
       DT INTERVAL)
       RETURNS SETOF TIMESTAMP WITHOUT TIME ZONE AS
$BODY$
SELECT $1 + 
       generate_float_series(0,extract(epoch FROM $2-$1),extract(epoch FROM $3))
	* INTERVAL '1 second';
$BODY$ LANGUAGE 'sql' IMMUTABLE;

CREATE OR REPLACE FUNCTION generate_timestamp_series(
       Tmin TIMESTAMP WITH TIME ZONE,
       Tmax TIMESTAMP WITH TIME ZONE,
       DT INTERVAL)
       RETURNS SETOF TIMESTAMP WITH TIME ZONE AS
$BODY$
SELECT $1 + 
       generate_float_series(0,extract(epoch FROM $2-$1),extract(epoch FROM $3))
	* INTERVAL '1 second';
$BODY$ LANGUAGE 'sql' IMMUTABLE;

--- Some functions to convert arrays of timestamps to arrays of epoch numbers
--- USE THE TIME ZONE STUFF CAREFULLY!!!
CREATE OR REPLACE FUNCTION extract_epoch_from_array(
       ts_array TIMESTAMP WITHOUT TIME ZONE[],tz text)
       RETURNS float[] AS
$BODY$
SELECT array_agg(a) FROM
       (SELECT extract( epoch FROM unnest($1) AT TIME ZONE $2 ) AS a) t;
$BODY$ LANGUAGE 'sql' IMMUTABLE;

CREATE OR REPLACE FUNCTION extract_epoch_from_array(
       ts_array TIMESTAMP WITH TIME ZONE[])
       RETURNS float[] AS
$BODY$
SELECT array_agg(a) FROM
       (SELECT extract( epoch FROM unnest($1) ) AS a) t;
$BODY$ LANGUAGE 'sql' IMMUTABLE;

-------------------------------------------------------------------
--- Now linear interpolation for timestamp arrays
-------------------------------------------------------------------
CREATE OR REPLACE FUNCTION LinearInterpolation(
       x TIMESTAMP WITHOUT TIME ZONE[], g geometry, 
       x2 TIMESTAMP WITHOUT TIME ZONE[], tz text, na boolean) 
       RETURNS geometry[] AS
$BODY$
SELECT LinearInterpolation( 
       extract_epoch_from_array($1,$4),
       $2,
       extract_epoch_from_array($3,$4),$5);
$BODY$ LANGUAGE 'sql' IMMUTABLE;

CREATE OR REPLACE FUNCTION LinearInterpolation(
       x TIMESTAMP WITHOUT TIME ZONE[], g geometry, 
       x2 TIMESTAMP WITHOUT TIME ZONE[], tz text) 
       RETURNS geometry AS
$BODY$
SELECT LinearInterpolation( 
       extract_epoch_from_array($1,$4),
       $2,
       extract_epoch_from_array($3,$4));
$BODY$ LANGUAGE 'sql' IMMUTABLE;

CREATE OR REPLACE FUNCTION LinearInterpolation(
       x TIMESTAMP WITH TIME ZONE[], g geometry, 
       x2 TIMESTAMP WITH TIME ZONE[], na boolean) 
       RETURNS geometry[] AS
$BODY$
SELECT LinearInterpolation( 
       extract_epoch_from_array($1),
       $2,
       extract_epoch_from_array($3),$4);
$BODY$ LANGUAGE 'sql' IMMUTABLE;

CREATE OR REPLACE FUNCTION LinearInterpolation(
       x TIMESTAMP WITH TIME ZONE[], g geometry, 
       x2 TIMESTAMP WITH TIME ZONE[]) 
       RETURNS geometry AS
$BODY$
SELECT LinearInterpolation( 
       extract_epoch_from_array($1),
       $2,
       extract_epoch_from_array($3));
$BODY$ LANGUAGE 'sql' IMMUTABLE;

--- Now linear interpolation for timestamp arrays on float arrays
--- instead of geometries
CREATE OR REPLACE FUNCTION LinearInterpolation(
       x TIMESTAMP WITHOUT TIME ZONE[], y float[], 
       x2 TIMESTAMP WITHOUT TIME ZONE[], tz text, na boolean) 
       RETURNS float[] AS
$BODY$
SELECT LinearInterpolation( 
       extract_epoch_from_array($1,$4),
       $2,
       extract_epoch_from_array($3,$4),$5);
$BODY$ LANGUAGE 'sql' IMMUTABLE;

CREATE OR REPLACE FUNCTION LinearInterpolation(
       x TIMESTAMP WITHOUT TIME ZONE[], y float[], 
       x2 TIMESTAMP WITHOUT TIME ZONE[], tz text) 
       RETURNS float[] AS
$BODY$
SELECT LinearInterpolation( 
       extract_epoch_from_array($1,$4),
       $2,
       extract_epoch_from_array($3,$4));
$BODY$ LANGUAGE 'sql' IMMUTABLE;

CREATE OR REPLACE FUNCTION LinearInterpolation(
       x TIMESTAMP WITH TIME ZONE[], y float[], 
       x2 TIMESTAMP WITH TIME ZONE[], na boolean) 
       RETURNS float[] AS
$BODY$
SELECT LinearInterpolation( 
       extract_epoch_from_array($1),
       $2,
       extract_epoch_from_array($3),$4);
$BODY$ LANGUAGE 'sql' IMMUTABLE;

CREATE OR REPLACE FUNCTION LinearInterpolation(
       x TIMESTAMP WITH TIME ZONE[], y float[], 
       x2 TIMESTAMP WITH TIME ZONE[]) 
       RETURNS float[] AS
$BODY$
SELECT LinearInterpolation( 
       extract_epoch_from_array($1),
       $2,
       extract_epoch_from_array($3));
$BODY$ LANGUAGE 'sql' IMMUTABLE;
检索自 "https://wiki.postgresql.ac.cn/index.php?title=Linear_Interpolation&oldid=32989"