Linear Interpolation
From PostgreSQL wiki
Jump to navigationJump to searchLinear 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;
