From 9db4b7f5192c3f946bc6c6dab4a956e3e9a39d16 Mon Sep 17 00:00:00 2001 From: abelvm Date: Tue, 14 Jun 2016 17:55:45 +0200 Subject: [PATCH] first commit --- doc/08_interpolation.md | 51 +++++++ src/pg/sql/08_interpolation.sql | 127 ++++++++++++++++++ .../test/expected/08_interpolation_test.out | 4 + src/pg/test/sql/08_interpolation_test.sql | 6 + 4 files changed, 188 insertions(+) create mode 100644 doc/08_interpolation.md create mode 100644 src/pg/sql/08_interpolation.sql create mode 100644 src/pg/test/expected/08_interpolation_test.out create mode 100644 src/pg/test/sql/08_interpolation_test.sql 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;