From 5b6d2c568b9799b7ebc899c79b7a46999219d47a Mon Sep 17 00:00:00 2001 From: abelvm Date: Mon, 15 Aug 2016 17:15:36 -0400 Subject: [PATCH] first commit --- doc/19_contour.md | 47 ++++++ src/pg/sql/19_contour.sql | 191 +++++++++++++++++++++++ src/pg/test/expected/19_contour_test.out | 2 + src/pg/test/sql/19_contour_test.sql | 14 ++ 4 files changed, 254 insertions(+) create mode 100644 doc/19_contour.md create mode 100644 src/pg/sql/19_contour.sql create mode 100644 src/pg/test/expected/19_contour_test.out create mode 100644 src/pg/test/sql/19_contour_test.sql diff --git a/doc/19_contour.md b/doc/19_contour.md new file mode 100644 index 0000000..f34f7c0 --- /dev/null +++ b/doc/19_contour.md @@ -0,0 +1,47 @@ +## Contour maps + +Function to generate a contour map from an 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_Contour (geom geometry[], values numeric[], resolution integer, buffer numeric, method, classmethod integer, steps integer) + +#### Arguments + +| Name | Type | Description | +|------|------|-------------| +| geom | geometry[] | Array of points's geometries | +| values | numeric[] | Array of points' values for the param under study| +| resolution | integer | Size in meters of the cells in the grid| +| buffer | numeric | Value between 0 and 1 for spatial buffer of the set of points +| method | integer | 0:nearest neighbor, 1: barycentric, 2: IDW| +| classmethod | integer | 0:equals, 1: heads&tails, 2:jenks, 3:quantiles | +| steps | integer | Number of steps in the classification| + +### Returns +Returns a table object + +| Name | Type | Description | +|------|------|-------------| +| the_geom | geometry | Geometries of the classified contour map| +| avg_value | numeric | Avg value of the area| +| min_value | numeric | Min value of the area| +| max_value | numeric | Max value of the areal| +| bin | integer | Index of the class of the area| + +#### Example Usage + +```sql +WITH a AS ( + SELECT + ARRAY[800, 700, 600, 500, 400, 300, 200, 100]::numeric[] 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 + foo.* +FROM + a, + CDB_contour(a.g, a.vals, 500, 0.0, 1,3,5) foo; +``` diff --git a/src/pg/sql/19_contour.sql b/src/pg/sql/19_contour.sql new file mode 100644 index 0000000..5118972 --- /dev/null +++ b/src/pg/sql/19_contour.sql @@ -0,0 +1,191 @@ +CREATE OR REPLACE FUNCTION CDB_Contour( + IN geomin geometry[], + IN colin numeric[], + IN resolution integer, + IN buffer numeric, + IN intmethod integer, + IN classmethod integer, + IN steps integer + ) +RETURNS TABLE( + the_geom geometry, + bin integer, + min_value numeric, + max_value numeric, + avg_value numeric +) AS $$ +DECLARE + cell numeric; + tin geometry[]; +BEGIN + -- calc the cell size in web mercator units + WITH center as ( + SELECT ST_centroid(ST_Collect(geomin)) as c + ) + SELECT + round(resolution / cos(ST_y(c) * pi()/180)) + INTO cell + FROM center; + -- raise notice 'Resol: %', cell; + + -- we don't have iterative barycentric interpolation in CDB_interpolation, + -- and it's a costy function, so let's make a custom one here till + -- we update the code + -- tin := ARRAY[]::geometry[]; + IF intmethod=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) + SELECT array_agg(v) INTO tin FROM c; + END IF; + -- Delaunay stuff performed just ONCE!! + + -- magic + RETURN QUERY + WITH + convexhull as ( + SELECT + ST_ConvexHull(ST_Collect(geomin)) as g, + buffer * |/ st_area(ST_ConvexHull(ST_Collect(geomin)))/PI() as r + ), + envelope as ( + SELECT + st_expand(a.g, a.r) as e + FROM convexhull a + ), + envelope3857 as( + SELECT + ST_Transform(e, 3857) as geom + FROM envelope + ), + grid as( + SELECT + ST_Transform(CDB_RectangleGrid(geom, cell, cell), 4326) as geom + FROM envelope3857 + ), + interp as( + SELECT + geom, + CASE + WHEN intmethod=1 THEN cdb_crankshaft._interp_in_tin(geomin, colin, tin, ST_Centroid(geom)) + ELSE db_crankshaft.CDB_SpatialInterpolation(geomin, colin, ST_Centroid(geom), intmethod) + END as val + FROM grid + ), + classes as( + SELECT CASE + WHEN classmethod = 0 THEN + CDB_EqualIntervalBins(array_agg(val), steps) + WHEN classmethod = 1 THEN + CDB_HeadsTailsBins(array_agg(val), steps) + WHEN classmethod = 2 THEN + CDB_JenksBins(array_agg(val), steps) + ELSE + CDB_QuantileBins(array_agg(val), steps) + END as b + FROM interp + where val is not null + ), + classified as( + SELECT + i.*, + width_bucket(i.val, c.b) as bucket + FROM interp i left join classes c + ON 1=1 + ), + classified2 as( + SELECT + geom, + val, + CASE + WHEN bucket = steps THEN bucket - 1 + ELSE bucket + END as b + FROM classified + ), + final as( + SELECT + st_union(geom) as the_geom, + b as bin, + min(val) as min_value, + max(val) as max_value, + avg(val) as avg_value + FROM classified2 + GROUP BY bin + ) + SELECT + * + FROM final + where final.bin is not null + ; +END; +$$ language plpgsql; + + + +-- ===================================================================== +-- Interp in grid, so we can use barycentric with a precalculated tin (NNI) +-- ===================================================================== +CREATE OR REPLACE FUNCTION cdb_crankshaft._interp_in_tin( + IN geomin geometry[], + IN colin numeric[], + IN tin geometry[], + IN point geometry + ) +RETURNS numeric AS +$$ +DECLARE + g geometry; + vertex geometry[]; + sg numeric; + sa numeric; + sb numeric; + sc numeric; + va numeric; + vb numeric; + vc numeric; + output numeric; +BEGIN + -- get the cell the point is within + WITH + a as (SELECT unnest(tin) as v), + b as (SELECT v FROM a WHERE ST_Within(point, v)) + SELECT v INTO g FROM b; + + -- if we're out of the data realm, + -- return null + IF g is null THEN + RETURN null; + 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(geomin) as geo, unnest(colin) as c) + SELECT c INTO va FROM a WHERE ST_Equals(geo, vertex[1]); + + WITH a AS(SELECT unnest(geomin) as geo, unnest(colin) as c) + SELECT c INTO vb FROM a WHERE ST_Equals(geo, vertex[2]); + + WITH a AS(SELECT unnest(geomin) as geo, unnest(colin) as c) + SELECT c INTO vc FROM a WHERE ST_Equals(geo, vertex[3]); + + -- calc the areas + 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,1); + RETURN output; +END; +$$ +language plpgsql; + + diff --git a/src/pg/test/expected/19_contour_test.out b/src/pg/test/expected/19_contour_test.out new file mode 100644 index 0000000..da0e94b --- /dev/null +++ b/src/pg/test/expected/19_contour_test.out @@ -0,0 +1,2 @@ +SET client_min_messages TO WARNING; +\set ECHO none diff --git a/src/pg/test/sql/19_contour_test.sql b/src/pg/test/sql/19_contour_test.sql new file mode 100644 index 0000000..fa32609 --- /dev/null +++ b/src/pg/test/sql/19_contour_test.sql @@ -0,0 +1,14 @@ +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]::numeric[] 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 + foo.* +FROM + a, + CDB_contour(a.g, a.vals, 500, 0.0, 1,3,5) foo;