smart guessing optional
This commit is contained in:
abelvm
123
src/pg/sql/13_PIA.sql
Normal file
123
src/pg/sql/13_PIA.sql
Normal file
@@ -0,0 +1,123 @@
|
||||
-- Based on:
|
||||
-- https://github.com/mapbox/polylabel/blob/master/index.js
|
||||
-- https://sites.google.com/site/polesofinaccessibility/
|
||||
-- Requires: https://github.com/CartoDB/cartodb-postgresql
|
||||
|
||||
-- Based on:
|
||||
-- https://github.com/mapbox/polylabel/blob/master/index.js
|
||||
-- https://sites.google.com/site/polesofinaccessibility/
|
||||
-- Requires: https://github.com/CartoDB/cartodb-postgresql
|
||||
|
||||
CREATE OR REPLACE FUNCTION CDB_PIA(
|
||||
IN polygon geometry,
|
||||
IN tolerance numeric DEFAULT 1.0
|
||||
)
|
||||
RETURNS geometry AS $$
|
||||
DECLARE
|
||||
env geometry[];
|
||||
cells geometry[];
|
||||
cell geometry;
|
||||
best_c geometry;
|
||||
best_d numeric;
|
||||
test_d numeric;
|
||||
test_mx numeric;
|
||||
test_h numeric;
|
||||
test_cells geometry[];
|
||||
width numeric;
|
||||
height numeric;
|
||||
h numeric;
|
||||
i integer;
|
||||
n integer;
|
||||
sqr numeric;
|
||||
p geometry;
|
||||
BEGIN
|
||||
sqr := |/2;
|
||||
polygon := ST_Transform(polygon, 3857);
|
||||
|
||||
-- grid #0 cell size
|
||||
height := ST_YMax(polygon) - ST_YMin(polygon);
|
||||
width := ST_XMax(polygon) - ST_XMin(polygon);
|
||||
h := 0.5*LEAST(height, width);
|
||||
|
||||
-- grid #0
|
||||
with c1 as(
|
||||
SELECT cdb_crankshaft.CDB_RectangleGrid(polygon, h, h) as c
|
||||
)
|
||||
SELECT array_agg(c) INTO cells FROM c1;
|
||||
|
||||
-- 1st guess: centroid
|
||||
best_d := cdb_crankshaft._Signed_Dist(polygon, ST_Centroid(Polygon));
|
||||
|
||||
-- looping the loop
|
||||
n := array_length(cells,1);
|
||||
i := 1;
|
||||
LOOP
|
||||
|
||||
EXIT WHEN i > n;
|
||||
|
||||
cell := cells[i];
|
||||
i := i+1;
|
||||
|
||||
-- cell side size, it's square
|
||||
test_h := ST_XMax(cell) - ST_XMin(cell) ;
|
||||
|
||||
-- check distance
|
||||
test_d := cdb_crankshaft._Signed_Dist(polygon, ST_Centroid(cell));
|
||||
IF test_d > best_d THEN
|
||||
best_d := test_d;
|
||||
best_c := cells[i];
|
||||
END IF;
|
||||
|
||||
-- longest distance within the cell
|
||||
test_mx := test_d + (test_h/2 * sqr);
|
||||
|
||||
-- if the cell has no chance to contains the desired point, continue
|
||||
CONTINUE WHEN test_mx - best_d <= tolerance;
|
||||
|
||||
-- resample the cell
|
||||
with c1 as(
|
||||
SELECT cdb_crankshaft.CDB_RectangleGrid(cell, test_h/2, test_h/2) as c
|
||||
)
|
||||
SELECT array_agg(c) INTO test_cells FROM c1;
|
||||
|
||||
-- concat the new cells to the former array
|
||||
cells := cells || test_cells;
|
||||
|
||||
-- prepare next iteration
|
||||
n := array_length(cells,1);
|
||||
|
||||
END LOOP;
|
||||
|
||||
RETURN ST_transform(ST_Centroid(best_c), 4326);
|
||||
|
||||
END;
|
||||
$$ language plpgsql IMMUTABLE;
|
||||
|
||||
|
||||
-- signed distance point to polygon with holes
|
||||
-- negative is the point is out the polygon
|
||||
CREATE OR REPLACE FUNCTION _Signed_Dist(
|
||||
IN polygon geometry,
|
||||
IN point geometry
|
||||
)
|
||||
RETURNS numeric AS $$
|
||||
DECLARE
|
||||
i integer;
|
||||
within integer;
|
||||
holes integer;
|
||||
dist numeric;
|
||||
BEGIN
|
||||
dist := 1e999;
|
||||
SELECT LEAST(dist, ST_distance(point, ST_ExteriorRing(polygon))::numeric) INTO dist;
|
||||
SELECT CASE WHEN ST_Within(point,polygon) THEN 1 ELSE -1 END INTO within;
|
||||
SELECT ST_NumInteriorRings(polygon) INTO holes;
|
||||
IF holes > 0 THEN
|
||||
FOR i IN 1..holes
|
||||
LOOP
|
||||
SELECT LEAST(dist, ST_distance(point, ST_InteriorRingN(polygon, i))::numeric) INTO dist;
|
||||
END LOOP;
|
||||
END IF;
|
||||
dist := dist * within::numeric;
|
||||
RETURN dist;
|
||||
END;
|
||||
$$ language plpgsql IMMUTABLE;
|
||||
67
src/pg/sql/14_densify.sql
Normal file
67
src/pg/sql/14_densify.sql
Normal file
@@ -0,0 +1,67 @@
|
||||
--
|
||||
-- Iterative densification of a set of points using Delaunay triangulation
|
||||
-- the new points have as assigned value the average value of the 3 vertex (centroid)
|
||||
--
|
||||
-- @param geomin - array of geometries (points)
|
||||
--
|
||||
-- @param colin - array of numeric values in that points
|
||||
--
|
||||
-- @param iterations - integer, number of iterations
|
||||
--
|
||||
--
|
||||
-- Returns: TABLE(geomout geometry, colout numeric)
|
||||
--
|
||||
--
|
||||
CREATE OR REPLACE FUNCTION CDB_Densify(
|
||||
IN geomin geometry[],
|
||||
IN colin numeric[],
|
||||
IN iterations integer
|
||||
)
|
||||
RETURNS TABLE(geomout geometry, colout numeric) AS $$
|
||||
DECLARE
|
||||
geotemp geometry[];
|
||||
coltemp numeric[];
|
||||
i integer;
|
||||
gs geometry[];
|
||||
g geometry;
|
||||
vertex geometry[];
|
||||
va numeric;
|
||||
vb numeric;
|
||||
vc numeric;
|
||||
center geometry;
|
||||
centerval numeric;
|
||||
tmp integer;
|
||||
BEGIN
|
||||
geotemp := geomin;
|
||||
coltemp := colin;
|
||||
FOR i IN 1..iterations
|
||||
LOOP
|
||||
-- generate TIN
|
||||
WITH a as (SELECT unnest(geotemp) 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 gs FROM c;
|
||||
-- loop cells
|
||||
FOREACH g IN ARRAY gs
|
||||
LOOP
|
||||
-- append centroid
|
||||
SELECT ST_Centroid(g) INTO center;
|
||||
geotemp := array_append(geotemp, center);
|
||||
-- retrieve the value of each vertex
|
||||
WITH a AS (SELECT (ST_DumpPoints(g)).geom AS v)
|
||||
SELECT array_agg(v) INTO vertex FROM a;
|
||||
WITH a AS(SELECT unnest(geotemp) as geo, unnest(coltemp) as c)
|
||||
SELECT c INTO va FROM a WHERE ST_Equals(geo, vertex[1]);
|
||||
WITH a AS(SELECT unnest(geotemp) as geo, unnest(coltemp) as c)
|
||||
SELECT c INTO vb FROM a WHERE ST_Equals(geo, vertex[2]);
|
||||
WITH a AS(SELECT unnest(geotemp) as geo, unnest(coltemp) as c)
|
||||
SELECT c INTO vc FROM a WHERE ST_Equals(geo, vertex[3]);
|
||||
-- calc the value at the center
|
||||
centerval := (va + vb + vc) / 3;
|
||||
-- append the value
|
||||
coltemp := array_append(coltemp, centerval);
|
||||
END LOOP;
|
||||
END LOOP;
|
||||
RETURN QUERY SELECT unnest(geotemp ) as geomout, unnest(coltemp ) as colout;
|
||||
END;
|
||||
$$ language plpgsql IMMUTABLE;
|
||||
43
src/pg/sql/15_tinmap.sql
Normal file
43
src/pg/sql/15_tinmap.sql
Normal file
@@ -0,0 +1,43 @@
|
||||
CREATE OR REPLACE FUNCTION CDB_TINmap(
|
||||
IN geomin geometry[],
|
||||
IN colin numeric[],
|
||||
IN iterations integer
|
||||
)
|
||||
RETURNS TABLE(geomout geometry, colout numeric) AS $$
|
||||
DECLARE
|
||||
p geometry[];
|
||||
vals numeric[];
|
||||
gs geometry[];
|
||||
g geometry;
|
||||
vertex geometry[];
|
||||
centerval numeric;
|
||||
va numeric;
|
||||
vb numeric;
|
||||
vc numeric;
|
||||
coltemp numeric[];
|
||||
BEGIN
|
||||
SELECT array_agg(dens.geomout), array_agg(dens.colout) INTO p, vals FROM cdb_crankshaft.CDB_Densify(geomin, colin, iterations) dens;
|
||||
WITH a as (SELECT unnest(p) 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 gs FROM c;
|
||||
FOREACH g IN ARRAY gs
|
||||
LOOP
|
||||
-- retrieve the vertex of each triangle
|
||||
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(p) as geo, unnest(vals) as c)
|
||||
SELECT c INTO va FROM a WHERE ST_Equals(geo, vertex[1]);
|
||||
WITH a AS(SELECT unnest(p) as geo, unnest(vals) as c)
|
||||
SELECT c INTO vb FROM a WHERE ST_Equals(geo, vertex[2]);
|
||||
WITH a AS(SELECT unnest(p) as geo, unnest(vals) as c)
|
||||
SELECT c INTO vc FROM a WHERE ST_Equals(geo, vertex[3]);
|
||||
-- calc the value at the center
|
||||
centerval := (va + vb + vc) / 3;
|
||||
-- append the value
|
||||
coltemp := array_append(coltemp, centerval);
|
||||
END LOOP;
|
||||
RETURN QUERY SELECT unnest(gs) as geomout, unnest(coltemp ) as colout;
|
||||
END;
|
||||
$$ language plpgsql IMMUTABLE;
|
||||
Reference in New Issue
Block a user