Linear Interpolation

From PostgreSQL wiki
Jump to navigationJump to search

Snippets

Linear Intperpolation

Works with PostgreSQL

9.0

Written in

SQL

Depends on

Nothing

Linear interpolation in PostgreSQL

The code below does simple linear interpolation in PostgreSQL based on arrays containing x and y values of the curve to be interpolated (or LINESTRING geometries for PostGIS trajectory interpolation). In particular, I use it to linearly interpolate trajectories of moving objects at specific times. The code is simple, but reasonably efficient and pure PL/pgSQL, so it should work on most installations (functions doing geospatial calculations require the PostGIS extension version >=2.0). Some of the code is redundant in newer versions of PostgreSQL (e.g., generating series of floats or timestamps can now be generated with the base generate_series function), but it still works and it should be reasonably obvious how to modernize the code if need be.

--- 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;