# Linear Interpolation

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