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..76fad01 --- /dev/null +++ b/src/pg/sql/08_interpolation.sql @@ -0,0 +1,130 @@ +-- 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[]; + output numeric; +BEGIN + EXECUTE 'WITH a AS('||query||') SELECT array_agg(the_geom), array_agg(attrib) FROM a' INTO gs, vs; + SELECT CDB_SpatialInterpolation(gs, vs, point, method, p1,p2) INTO output FROM a; + + RETURN output; +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..635ca2a --- /dev/null +++ b/src/pg/test/expected/08_interpolation_test.out @@ -0,0 +1,5 @@ +SET client_min_messages TO WARNING; +\set ECHO none +cdb_spatialinterpolation +t +(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..bd9c729 --- /dev/null +++ b/src/pg/test/sql/08_interpolation_test.sql @@ -0,0 +1,10 @@ +SET client_min_messages TO WARNING; +\set ECHO none +\pset format unaligned + +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_crankshaft.CDB_SpatialInterpolation(g, vals, ST_GeomFromText('POINT(2.154 41.37)'), 1) - 780.79470198683925288365) / 780.79470198683925288365 < 0.001 As cdb_spatialinterpolation FROM a; diff --git a/src/py/crankshaft/crankshaft/clustering/moran.py b/src/py/crankshaft/crankshaft/clustering/moran.py index 39b3ff6..4bced89 100644 --- a/src/py/crankshaft/crankshaft/clustering/moran.py +++ b/src/py/crankshaft/crankshaft/clustering/moran.py @@ -7,6 +7,7 @@ Moran's I geostatistics (global clustering & outliers presence) import pysal as ps import plpy +from collections import OrderedDict # crankshaft module import crankshaft.pysal_utils as pu @@ -21,11 +22,11 @@ def moran(subquery, attr_name, core clusters with PySAL. Andy Eschbacher """ - qvals = {"id_col": id_col, - "attr1": attr_name, - "geom_col": geom_col, - "subquery": subquery, - "num_ngbrs": num_ngbrs} + qvals = OrderedDict([("id_col", id_col), + ("attr1", attr_name), + ("geom_col", geom_col), + ("subquery", subquery), + ("num_ngbrs", num_ngbrs)]) query = pu.construct_neighbor_query(w_type, qvals) @@ -65,11 +66,11 @@ def moran_local(subquery, attr, # geometries with attributes that are null are ignored # resulting in a collection of not as near neighbors - qvals = {"id_col": id_col, - "attr1": attr, - "geom_col": geom_col, - "subquery": subquery, - "num_ngbrs": num_ngbrs} + qvals = OrderedDict([("id_col", id_col), + ("attr1", attr), + ("geom_col", geom_col), + ("subquery", subquery), + ("num_ngbrs", num_ngbrs)]) query = pu.construct_neighbor_query(w_type, qvals) @@ -101,12 +102,12 @@ def moran_rate(subquery, numerator, denominator, Moran's I Rate (global) Andy Eschbacher """ - qvals = {"id_col": id_col, - "attr1": numerator, - "attr2": denominator, - "geom_col": geom_col, - "subquery": subquery, - "num_ngbrs": num_ngbrs} + qvals = OrderedDict([("id_col", id_col), + ("attr1", numerator), + ("attr2", denominator) + ("geom_col", geom_col), + ("subquery", subquery), + ("num_ngbrs", num_ngbrs)]) query = pu.construct_neighbor_query(w_type, qvals) @@ -145,13 +146,14 @@ def moran_local_rate(subquery, numerator, denominator, # geometries with values that are null are ignored # resulting in a collection of not as near neighbors - query = pu.construct_neighbor_query(w_type, - {"id_col": id_col, - "numerator": numerator, - "denominator": denominator, - "geom_col": geom_col, - "subquery": subquery, - "num_ngbrs": num_ngbrs}) + qvals = OrderedDict([("id_col", id_col), + ("numerator", numerator), + ("denominator", denominator), + ("geom_col", geom_col), + ("subquery", subquery), + ("num_ngbrs", num_ngbrs)]) + + query = pu.construct_neighbor_query(w_type, qvals) try: result = plpy.execute(query) @@ -186,12 +188,12 @@ def moran_local_bv(subquery, attr1, attr2, """ plpy.notice('** Constructing query') - qvals = {"num_ngbrs": num_ngbrs, - "attr1": attr1, - "attr2": attr2, - "subquery": subquery, - "geom_col": geom_col, - "id_col": id_col} + qvals = OrderedDict([("id_col", id_col), + ("attr1", attr1), + ("attr2", attr2), + ("geom_col", geom_col), + ("subquery", subquery), + ("num_ngbrs", num_ngbrs)]) query = pu.construct_neighbor_query(w_type, qvals)