cartodb4/crankshaft
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

working versiongit branch

Browse Source
This commit is contained in:
abelvm
2017-02-16 16:51:07 +01:00
parent bd03f3483f
commit fb5744fccd
1 changed files with 94 additions and 36 deletions
Download Patch File Download Diff File Expand all files Collapse all files

130
src/pg/sql/12_contour2.sql
View File

@@ -1,13 +1,3 @@
-- =============================================================================================
--
-- CDB_Contour v2
--
-- Meander triangles
--
-- author: Abel Vázquez
--
-- =============================================================================================
CREATE OR REPLACE FUNCTION CDB_contour2(
IN geomin geometry[],
IN colin numeric[],
@@ -30,8 +20,10 @@ DECLARE
interp31 geometry[];
segment geometry[];
running_merge geometry[];
-- interlines geometry;
geomout geometry;
i integer;
cell_breaks integer[];
running_breaks integer[];
BEGIN
-- generate the breaks
@@ -56,6 +48,9 @@ BEGIN
c as (SELECT (ST_Dump(t)).geom AS v FROM b)
SELECT array_agg(v) INTO gs FROM c;
-- RAISE NOTICE 'TIN size: %',array_length(gs,1);
i:= 0;
-- ======================================================================================
-- loop in the TIN
@@ -74,89 +69,133 @@ BEGIN
WITH a AS(SELECT unnest(geomin) as geo, unnest(colin) as c)
SELECT c INTO v3 FROM a WHERE ST_Equals(geo, vertex[3]);
-- RAISE NOTICE 'Cell: %', i;
-- i := i +1;
-- retrieve the bucket index of each vertex
bu := ARRAY[width_bucket(v1, breaks), width_bucket(v2, breaks), width_bucket(v3, breaks)];
-- RAISE NOTICE ' - BU: %',bu;
-- continue when there is no contour line crossing the current cell
CONTINUE WHEN bu[1] = bu[2] and bu[1] = bu[3];
-- we have contour lines in this cell, let's find their intersections with triangle sides
interp12 := array_fill(null::geometry, ARRAY[steps]);
IF bu[1] <> bu[2] THEN
SELECT array_agg((t.x-v1)/(v2-v1)) INTO inter FROM unnest(breaks) as t(x);
-- RAISE NOTICE ' - Inter 12: %', inter;
WITH a as(
SELECT
ST_Line_Interpolate_Point(ST_MakeLine(vertex[1], vertex[2]),t.x)
CASE WHEN t.x BETWEEN 0 AND 1 THEN
ST_LineInterpolatePoint(ST_MakeLine(vertex[1], vertex[2]),t.x)
ELSE null::geometry END as point
FROM unnest(inter) as t(x)
WHERE t.x BETWEEN 0 AND 1
)
SELECT
array_agg(g) INTO interp12
array_agg(point) INTO interp12
FROM a;
END IF;
interp23 := array_fill(null::geometry, ARRAY[steps]);
IF bu[2] <> bu[3] THEN
SELECT array_agg((t.x-v2)/(v3-v2)) INTO inter FROM unnest(breaks) as t(x);
-- RAISE NOTICE ' - Inter 23: %', inter;
WITH a as(
SELECT
ST_Line_Interpolate_Point(ST_MakeLine(vertex[2], vertex[3]),t.x)
CASE WHEN t.x BETWEEN 0 AND 1 THEN
ST_LineInterpolatePoint(ST_MakeLine(vertex[2], vertex[3]),t.x)
ELSE null::geometry END as point
FROM unnest(inter) as t(x)
WHERE t.x BETWEEN 0 AND 1
)
SELECT
array_agg(g) INTO interp23
array_agg(point) INTO interp23
FROM a;
END IF;
interp31 := array_fill(null::geometry, ARRAY[steps]);
IF bu[3] <> bu[1] THEN
SELECT array_agg((t.x-v3)/(v1-v3)) INTO inter FROM unnest(breaks) as t(x);
-- RAISE NOTICE ' - Inter 31: %', inter;
WITH a as(
SELECT
ST_Line_Interpolate_Point(ST_MakeLine(vertex[3], vertex[1]),t.x)
CASE WHEN t.x BETWEEN 0 AND 1 THEN
ST_LineInterpolatePoint(ST_MakeLine(vertex[3], vertex[1]),t.x)
ELSE null::geometry END as point
FROM unnest(inter) as t(x)
WHERE t.x BETWEEN 0 AND 1
)
SELECT
array_agg(g) INTO interp31
array_agg(point) INTO interp31
FROM a;
END IF;
-- create segments per bucket
-- create segments crossing the cell per bucket
WITH
a as(
SELECT
generate_series(1,steps) as break,
unnest(interp12) as p12,
unnest(interp23) as p23,
unnest(interp31) as p31
),
b as(
SELECT
break,
ST_MakeLine(ARRAY[p12, p23, p31]::geometry[]) as segm
FROM a
WHERE
(p12 is not null and p23 is not null and ST_equals(p12, p23)=false) OR
(p23 is not null and p31 is not null and ST_equals(p23, p31)=false) OR
(p31 is not null and p12 is not null and ST_equals(p31, p12)=false)
)
SELECT
array_agg(segm) INTO segment
array_agg(break), array_agg(segm) INTO cell_breaks, segment
FROM b
WHERE
segm IS NOT null AND
ST_IsEmpty(segm) = false;
-- concat the segments
SELECT
CASE
WHEN array_length(running_merge,1)=0 THEN segment
ELSE array_cat(running_merge, segment) END
INTO running_merge;
-- RAISE NOTICE ' - Segments: %',array_length(segment,1);
-- concat the segments and breaks
IF array_length(running_merge,1)=0 THEN
running_merge := segment;
running_breaks := cell_breaks;
ELSE
running_merge := running_merge || segment ;
running_breaks := running_breaks || cell_breaks;
END IF;
-- loop end
END LOOP;
segment := unnest(running_merge);
RETURN ST_Collect(segment);
-- RAISE NOTICE 'TOTAL segment: %',array_length(running_merge,1);
-- ======================================================================================
-- merge and polygonize the results
running_merge := running_merge || ST_ExteriorRing(ST_Convexhull(ST_Collect(geomin))) ;
WITH
a as(
SELECT unnest(running_merge) as geo
)
SELECT ST_collect(geo) into geomout from a;
-- return some stuff
RETURN geomout;
END;
$$ language plpgsql IMMUTABLE;
-- ////////////////////////////////////////////////////////////
-- merge and polygonize the results
WITH
partials as(
SELECT
@@ -197,8 +236,27 @@ BEGIN
ST_GeometryType(geo) = 'ST_Polygon';
-- return some stuff
RETURN geomout;
-- ////////////////////////////////////////////////////////////
END;
$$ language plpgsql IMMUTABLE;
with
a as(
SELECT
array_agg(the_geom) as geomin,
array_agg(real_capacity) as colin
FROM uk_fcc
),
b as(
SELECT
CDB_contour2(
geomin,
colin,
2,
7
) as the_geom
from a
)
SELECT
*,
st_transform(the_geom::geometry, 3857) as the_geom_webmercator
from b;