diff --git a/doc/08_interpolation.md b/doc/08_interpolation.md
new file mode 100644
index 0000000..22fc1bc
--- /dev/null
+++ b/doc/08_interpolation.md
@@ -0,0 +1,51 @@
+## Spacial interpolation
+
+Function to interpolate a numeric attribute of a point in a scatter dataset of points, using one of three methos:
+
+* [Nearest neighbor](https://en.wikipedia.org/wiki/Nearest-neighbor_interpolation)
+* [Barycentric](https://en.wikipedia.org/wiki/Barycentric_coordinate_system)
+* [IDW](https://en.wikipedia.org/wiki/Inverse_distance_weighting)
+
+### CDB_SpatialInterpolation (query text, point geometry, method integer DEFAULT 1, p1 integer DEFAULT 0, ps integer DEFAULT 0)
+
+#### Arguments
+
+| Name | Type | Description |
+|------|------|-------------|
+| query   | text | query that returns at least `the_geom` and a numeric value as `attrib` |
+| point   | geometry | The target point to calc the value |
+| method   | integer     | 0:nearest neighbor, 1: barycentric, 2: IDW|
+| p1   | integer     | IDW: limit the number of neighbors, 0->no limit|
+| p2   | integer     | IDW: order of distance decay, 0-> order 1|
+
+### CDB_SpatialInterpolation (geom geometry[], values numeric[], point geometry, method integer DEFAULT 1, p1 integer DEFAULT 0, ps integer DEFAULT 0)
+
+#### Arguments
+
+| Name | Type | Description |
+|------|------|-------------|
+| geom   | geometry[]  | Array of points's geometries |
+| values | numeric[]   | Array of points' values for the param under study|
+| point   | geometry | The target point to calc the value |
+| method   | integer     | 0:nearest neighbor, 1: barycentric, 2: IDW|
+| p1   | integer     | IDW: limit the number of neighbors, 0->no limit|
+| p2   | integer     | IDW: order of distance decay, 0-> order 1|
+
+### Returns
+
+| Column Name | Type | Description |
+|-------------|------|-------------|
+| value  | numeric | Interpolated value at the given point, `-888.888` if the given point is out of the boundaries of the source points set |
+
+
+#### Example Usage
+
+```sql
+with a as (
+    select
+        array_agg(the_geom) as geomin,
+        array_agg(temp::numeric) as colin
+    from table_4804232032
+)
+SELECT CDB_SpatialInterpolation(geomin, colin, CDB_latlng(41.38, 2.15),1) FROM a;
+```
diff --git a/src/pg/sql/08_interpolation.sql b/src/pg/sql/08_interpolation.sql
new file mode 100644
index 0000000..04f1584
--- /dev/null
+++ b/src/pg/sql/08_interpolation.sql
@@ -0,0 +1,127 @@
+-- 0: nearest neighbor
+-- 1: barymetric
+-- 2: IDW
+
+CREATE OR REPLACE FUNCTION CDB_SpatialInterpolation(
+    IN query text,
+    IN point geometry,
+    IN method integer DEFAULT 1,
+    IN p1 numeric DEFAULT 0,
+    IN p2 numeric DEFAULT 0
+    )
+RETURNS numeric AS
+$$
+DECLARE
+    gs geometry[];
+    vs numeric[];
+BEGIN
+    EXECUTE 'WITH a AS('||query||') SELECT array_agg(the_geom), array_agg(attrib) FROM a' INTO gs, vs;
+    RETURN QUERY SELECT CDB_SpatialInterpolation(gs, vs, point, method, p1,p2) FROM a;
+END;
+$$
+language plpgsql IMMUTABLE;
+
+CREATE OR REPLACE FUNCTION CDB_SpatialInterpolation(
+    IN geomin geometry[],
+    IN colin numeric[],
+    IN point geometry,
+    IN method integer DEFAULT 1,
+    IN p1 numeric DEFAULT 0,
+    IN p2 numeric DEFAULT 0
+    )
+RETURNS numeric AS
+$$
+DECLARE
+    gs geometry[];
+    vs numeric[];
+    gs2 geometry[];
+    vs2 numeric[];
+    g geometry;
+    vertex geometry[];
+    sg numeric;
+    sa numeric;
+    sb numeric;
+    sc numeric;
+    va numeric;
+    vb numeric;
+    vc numeric;
+    output numeric;
+BEGIN
+    output :=  -999.999;
+    -- nearest
+    IF method = 0 THEN
+
+        WITH    a as (SELECT unnest(geomin) as g, unnest(colin) as v)
+        SELECT a.v INTO output FROM a ORDER BY point<->a.g LIMIT 1;
+        RETURN output;
+
+    -- barymetric
+    ELSIF method = 1 THEN
+        WITH    a as (SELECT unnest(geomin) AS e),
+                b as (SELECT ST_DelaunayTriangles(ST_Collect(a.e),0.001, 0) AS t FROM a),
+                c as (SELECT (ST_Dump(t)).geom as v FROM b),
+                d as (SELECT v FROM c WHERE ST_Within(point, v))
+            SELECT v INTO g FROM d;
+        IF g is null THEN
+            -- out of the realm of the input data
+            RETURN -888.888;
+        END IF;
+        -- vertex of the selected cell
+        WITH a AS (SELECT (ST_DumpPoints(g)).geom AS v)
+                SELECT array_agg(v) INTO vertex FROM a;
+
+            -- retrieve the value of each vertex
+        WITH a AS(SELECT unnest(vertex) as geo, unnest(colin) as c)
+            SELECT c INTO va FROM a WHERE ST_Equals(geo, vertex[1]);
+        WITH a AS(SELECT unnest(vertex) as geo, unnest(colin) as c)
+            SELECT c INTO vb FROM a WHERE ST_Equals(geo, vertex[2]);
+        WITH a AS(SELECT unnest(vertex) as geo, unnest(colin) as c)
+                SELECT c INTO vc FROM a WHERE ST_Equals(geo, vertex[3]);
+
+        SELECT ST_area(g), ST_area(ST_MakePolygon(ST_MakeLine(ARRAY[point, vertex[2], vertex[3], point]))), ST_area(ST_MakePolygon(ST_MakeLine(ARRAY[point, vertex[1], vertex[3], point]))), ST_area(ST_MakePolygon(ST_MakeLine(ARRAY[point,vertex[1],vertex[2], point]))) INTO sg, sa, sb, sc;
+
+        output := (coalesce(sa,0) * coalesce(va,0) + coalesce(sb,0) * coalesce(vb,0) + coalesce(sc,0) * coalesce(vc,0)) / coalesce(sg);
+        RETURN output;
+
+    -- IDW
+    -- p1: limit the number of neighbors, 0->no limit
+    -- p2: order of distance decay, 0-> order 1
+    ELSIF method = 2 THEN
+
+        IF p2 = 0 THEN
+            p2 := 1;
+        END IF;
+
+        WITH    a as (SELECT unnest(geomin) as g, unnest(colin) as v),
+                b as (SELECT a.g, a.v FROM a ORDER BY point<->a.g)
+        SELECT array_agg(b.g), array_agg(b.v) INTO gs, vs FROM b;
+        IF p1::integer>0 THEN
+            gs2:=gs;
+            vs2:=vs;
+            FOR i IN 1..p1
+            LOOP
+                gs2 := gs2 || gs[i];
+                vs2 := vs2 || vs[i];
+            END LOOP;
+        ELSE
+            gs2:=gs;
+            vs2:=vs;
+        END IF;
+
+        WITH    a as (SELECT unnest(gs2) as g, unnest(vs2) as v),
+                b as (
+                    SELECT
+                    (1/ST_distance(point, a.g)^p2::integer) as k,
+                    (a.v/ST_distance(point, a.g)^p2::integer) as f
+                    FROM a
+                )
+        SELECT sum(b.f)/sum(b.k) INTO output FROM b;
+        RETURN output;
+
+    END IF;
+
+    RETURN -777.777;
+
+END;
+$$
+language plpgsql IMMUTABLE;
diff --git a/src/pg/test/expected/08_interpolation_test.out b/src/pg/test/expected/08_interpolation_test.out
new file mode 100644
index 0000000..42d24cb
--- /dev/null
+++ b/src/pg/test/expected/08_interpolation_test.out
@@ -0,0 +1,4 @@
+ cdb_spatialinterpolation
+--------------------------
+ 780.79470198683925288365
+(1 row)
diff --git a/src/pg/test/sql/08_interpolation_test.sql b/src/pg/test/sql/08_interpolation_test.sql
new file mode 100644
index 0000000..c8db89d
--- /dev/null
+++ b/src/pg/test/sql/08_interpolation_test.sql
@@ -0,0 +1,6 @@
+WITH a AS (
+    SELECT
+    ARRAY[800, 700, 600, 500, 400, 300, 200, 100] AS vals,
+    ARRAY[ST_GeomFromText('POINT(2.1744 41.403)'),ST_GeomFromText('POINT(2.1228 41.380)'),ST_GeomFromText('POINT(2.1511 41.374)'),ST_GeomFromText('POINT(2.1528 41.413)'),ST_GeomFromText('POINT(2.165 41.391)'),ST_GeomFromText('POINT(2.1498 41.371)'),ST_GeomFromText('POINT(2.1533 41.368)'),ST_GeomFromText('POINT(2.131386 41.41399)')] AS g
+)
+SELECT CDB_SpatialInterpolation(g, vals, ST_GeomFromText('POINT(2.154 41.37)'),1) FROM a;