Linear Interpolation

== 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;
Category:PL/pgSQL