diff --git a/doc/08_interpolation.md b/doc/08_interpolation.md
index 22fc1bc..87b5124 100644
--- a/doc/08_interpolation.md
+++ b/doc/08_interpolation.md
@@ -2,7 +2,7 @@
 
 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)
+* [Nearest neighbor(s)](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)
 
@@ -15,7 +15,7 @@ Function to interpolate a numeric attribute of a point in a scatter dataset of p
 | 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|
+| p1   | integer     | limit the number of neighbors, IDW: 0->no limit, NN: 0-> closest one|
 | 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)
@@ -28,7 +28,7 @@ Function to interpolate a numeric attribute of a point in a scatter dataset of p
 | 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|
+| p1   | integer     | limit the number of neighbors, IDW: 0->no limit, NN: 0-> closest one|
 | p2   | integer     | IDW: order of distance decay, 0-> order 1|
 
 ### Returns
diff --git a/src/pg/sql/08_interpolation.sql b/src/pg/sql/08_interpolation.sql
index 0937f09..608583d 100644
--- a/src/pg/sql/08_interpolation.sql
+++ b/src/pg/sql/08_interpolation.sql
@@ -1,4 +1,4 @@
--- 0: nearest neighbor
+-- 0: nearest neighbor (s)
 -- 1: barymetric
 -- 2: IDW
 
@@ -51,11 +51,17 @@ DECLARE
     output numeric;
 BEGIN
     output :=  -999.999;
-    -- nearest
+
+    -- nearest neighbors
+    -- p1: limit the number of neighbors, 0-> closest one
     IF method = 0 THEN
 
+        IF p1 = 0 THEN
+            p1 := 1;
+        END IF;
+
         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;
+        SELECT avg(a.v) INTO output FROM a ORDER BY point<->a.g LIMIT p1::integer;
         RETURN output;
 
     -- barymetric
@@ -121,6 +127,11 @@ BEGIN
         SELECT sum(b.f)/sum(b.k) INTO output FROM b;
         RETURN output;
 
+    -- krigin
+    ELSIF method = 3 THEN
+
+
+
     END IF;
 
     RETURN -777.777;