missed files

This commit is contained in:
abelvm 2016-08-30 05:50:20 -04:00
parent 605e6e88ae
commit 72ee1a6d48
4 changed files with 296 additions and 0 deletions

38
doc/09_voronoi.md Normal file
View File

@ -0,0 +1,38 @@
## Voronoi
Function to construct the [Voronoi Diagram](https://en.wikipedia.org/wiki/Voronoi_diagram) from a dataset of scatter points, clipped to the significant area
PostGIS wil include this in future versions ([doc for dev branch](http://postgis.net/docs/manual-dev/ST_Voronoi.html)) and will perform faster for sure, but in the meantime...
### CDB_Voronoi (geom geometry[], buffer numeric DEFAULT 0.5, tolerance numeric DEFAULT 1e-9)
#### Arguments
| Name | Type | Description |
|------|------|-------------|
| geom | geometry[] | Array of points's geometries |
| buffer | numeric | enlargment ratio for the envelope area used for the restraints|
| tolerance | numeric | Delaunay tolerance, optional |
### Returns
| Column Name | Type | Description |
|-------------|------|-------------|
| geom | geometry collection | Collection of polygons of the Voronoi cells|
#### Example Usage
```sql
WITH a AS (
SELECT
ARRAY[ST_GeomFromText('POINT(2.1744 41.403)', 4326),ST_GeomFromText('POINT(2.1228 41.380)', 4326),ST_GeomFromText('POINT(2.1511 41.374)', 4326),ST_GeomFromText('POINT(2.1528 41.413)', 4326),ST_GeomFromText('POINT(2.165 41.391)', 4326),ST_GeomFromText('POINT(2.1498 41.371)', 4326),ST_GeomFromText('POINT(2.1533 41.368)', 4326),ST_GeomFromText('POINT(2.131386 41.41399)', 4326)] AS geomin
)
SELECT
st_transform(
(st_dump(CDB_voronoi(geomin, 0.2, 1e-9)
)).geom
, 3857) as the_geom_webmercator
FROM a;
```

236
src/pg/sql/09_voronoi.sql Normal file
View File

@ -0,0 +1,236 @@
-- =============================================================================================
--
-- CDB_Voronoi
--
-- =============================================================================================
CREATE OR REPLACE FUNCTION CDB_voronoi(
IN geomin geometry[],
IN buffer numeric DEFAULT 0.5,
IN tolerance numeric DEFAULT 1e-9
)
RETURNS geometry AS $$
DECLARE
geomout geometry;
BEGIN
-- we need to make the geometry calculations in (pseudo)meters!!!
with a as (
SELECT unnest(geomin) as g1
),
b as(
SELECT st_transform(g1, 3857) g2 from a
)
SELECT array_agg(g2) INTO geomin from b;
WITH
convexhull_1 as (
SELECT
ST_ConvexHull(ST_Collect(geomin)) as g,
buffer * |/ (st_area(ST_ConvexHull(ST_Collect(geomin)))/PI()) as r
),
clipper as(
SELECT
st_buffer(ST_MinimumBoundingCircle(a.g), buffer*a.r) as g
FROM convexhull_1 a
),
env0 as (
SELECT
(st_dumppoints(st_expand(a.g, buffer*a.r))).geom as e
FROM convexhull_1 a
),
env as (
SELECT
array_agg(env0.e) as e
FROM env0
),
sample AS (
SELECT
ST_Collect(geomin || env.e) as geom
FROM env
),
convexhull as (
SELECT
ST_ConvexHull(ST_Collect(geomin)) as cg
),
tin as (
SELECT
ST_Dump(ST_DelaunayTriangles(geom, tolerance, 0)) as gd
FROM
sample
),
tin_polygons as (
SELECT
(gd).Path as id,
(gd).Geom as pg,
ST_Centroid(ST_MinimumBoundingCircle((gd).Geom, 180)) as ct
FROM tin
),
tin_lines as (
SELECT
id,
ST_ExteriorRing(pg) as lg
FROM tin_polygons
),
tin_nodes as (
SELECT
id,
ST_PointN(lg,1) p1,
ST_PointN(lg,2) p2,
ST_PointN(lg,3) p3
FROM tin_lines
),
tin_edges AS (
SELECT
p.id,
UNNEST(ARRAY[
ST_MakeLine(n.p1,n.p2) ,
ST_MakeLine(n.p2,n.p3) ,
ST_MakeLine(n.p3,n.p1)]) as Edge,
ST_Force2D(cdb_crankshaft._Find_Circle(n.p1,n.p2,n.p3)) as ct,
CASE WHEN st_distance(p.ct, ST_ExteriorRing(p.pg)) < tolerance THEN
TRUE
ELSE FALSE END AS ctx,
p.pg,
ST_within(p.ct, convexhull.cg) as ctin
FROM
tin_polygons p,
tin_nodes n,
convexhull
WHERE p.id = n.id
),
voro_nodes as (
SELECT
CASE WHEN x.ctx = TRUE THEN
ST_Centroid(x.edge)
ELSE
x.ct
END as xct,
CASE WHEN y.id is null THEN
CASE WHEN x.ctin = TRUE THEN
ST_SetSRID(ST_MakePoint(
ST_X(x.ct) + ((ST_X(ST_Centroid(x.edge)) - ST_X(x.ct)) * (1+buffer)),
ST_Y(x.ct) + ((ST_Y(ST_Centroid(x.edge)) - ST_Y(x.ct)) * (1+buffer))
), ST_SRID(x.ct))
END
ELSE
y.ct
END as yct
FROM
tin_edges x
LEFT OUTER JOIN
tin_edges y
ON x.id <> y.id AND ST_Equals(x.edge, y.edge)
),
voro_edges as(
SELECT
ST_LineMerge(ST_Collect(ST_MakeLine(xct, yct))) as v
FROM
voro_nodes
),
voro_cells as(
SELECT
ST_Polygonize(
ST_Node(
ST_LineMerge(
ST_Union(v, ST_ExteriorRing(
ST_Convexhull(v)
)
)
)
)
) as g
FROM
voro_edges
),
voro_set as(
SELECT
(st_dump(v.g)).geom as g
FROM voro_cells v
),
clipped_voro as(
SELECT
ST_intersection(c.g, v.g) as g
FROM
voro_set v,
clipper c
WHERE
ST_GeometryType(v.g) = 'ST_Polygon'
)
SELECT
st_collect(
ST_Transform(
ST_ConvexHull(g),
4326
)
)
INTO geomout
FROM
clipped_voro;
RETURN geomout;
END;
$$ language plpgsql IMMUTABLE;
/** ----------------------------------------------------------------------------------------
* @function : FindCircle
* @precis : Function that determines if three points form a circle. If so a table containing
* centre and radius is returned. If not, a null table is returned.
* @version : 1.0.1
* @param : p_pt1 : First point in curve
* @param : p_pt2 : Second point in curve
* @param : p_pt3 : Third point in curve
* @return : geometry : In which X,Y ordinates are the centre X, Y and the Z being the radius of found circle
* or NULL if three points do not form a circle.
* @history : Simon Greener - Feb 2012 - Original coding.
* Rafa de la Torre - Aug 2016 - Small fix for type checking
* @copyright : Simon Greener @ 2012
* Licensed under a Creative Commons Attribution-Share Alike 2.5 Australia License. (http://creativecommons.org/licenses/by-sa/2.5/au/)
**/
CREATE OR REPLACE FUNCTION _Find_Circle(
IN p_pt1 geometry,
IN p_pt2 geometry,
IN p_pt3 geometry)
RETURNS geometry AS
$BODY$
DECLARE
v_Centre geometry;
v_radius NUMERIC;
v_CX NUMERIC;
v_CY NUMERIC;
v_dA NUMERIC;
v_dB NUMERIC;
v_dC NUMERIC;
v_dD NUMERIC;
v_dE NUMERIC;
v_dF NUMERIC;
v_dG NUMERIC;
BEGIN
IF ( p_pt1 IS NULL OR p_pt2 IS NULL OR p_pt3 IS NULL ) THEN
RAISE EXCEPTION 'All supplied points must be not null.';
RETURN NULL;
END IF;
IF ( ST_GeometryType(p_pt1) <> 'ST_Point' OR
ST_GeometryType(p_pt2) <> 'ST_Point' OR
ST_GeometryType(p_pt3) <> 'ST_Point' ) THEN
RAISE EXCEPTION 'All supplied geometries must be points.';
RETURN NULL;
END IF;
v_dA := ST_X(p_pt2) - ST_X(p_pt1);
v_dB := ST_Y(p_pt2) - ST_Y(p_pt1);
v_dC := ST_X(p_pt3) - ST_X(p_pt1);
v_dD := ST_Y(p_pt3) - ST_Y(p_pt1);
v_dE := v_dA * (ST_X(p_pt1) + ST_X(p_pt2)) + v_dB * (ST_Y(p_pt1) + ST_Y(p_pt2));
v_dF := v_dC * (ST_X(p_pt1) + ST_X(p_pt3)) + v_dD * (ST_Y(p_pt1) + ST_Y(p_pt3));
v_dG := 2.0 * (v_dA * (ST_Y(p_pt3) - ST_Y(p_pt2)) - v_dB * (ST_X(p_pt3) - ST_X(p_pt2)));
-- If v_dG is zero then the three points are collinear and no finite-radius
-- circle through them exists.
IF ( v_dG = 0 ) THEN
RETURN NULL;
ELSE
v_CX := (v_dD * v_dE - v_dB * v_dF) / v_dG;
v_CY := (v_dA * v_dF - v_dC * v_dE) / v_dG;
v_Radius := SQRT(POWER(ST_X(p_pt1) - v_CX,2) + POWER(ST_Y(p_pt1) - v_CY,2) );
END IF;
RETURN ST_SetSRID(ST_MakePoint(v_CX, v_CY, v_radius),ST_Srid(p_pt1));
END;
$BODY$
LANGUAGE plpgsql VOLATILE STRICT;

View File

@ -0,0 +1,7 @@
SET client_min_messages TO WARNING;
\set ECHO none
avg_area
----------------------
0.000178661700690617
(1 row)

View File

@ -0,0 +1,15 @@
SET client_min_messages TO WARNING;
\set ECHO none
WITH a AS (
SELECT
ARRAY[ST_GeomFromText('POINT(2.1744 41.403)', 4326),ST_GeomFromText('POINT(2.1228 41.380)', 4326),ST_GeomFromText('POINT(2.1511 41.374)', 4326),ST_GeomFromText('POINT(2.1528 41.413)', 4326),ST_GeomFromText('POINT(2.165 41.391)', 4326),ST_GeomFromText('POINT(2.1498 41.371)', 4326),ST_GeomFromText('POINT(2.1533 41.368)', 4326),ST_GeomFromText('POINT(2.131386 41.41399)', 4326)] AS geomin
),
b as(
SELECT
(st_dump(cdb_crankshaft.CDB_voronoi(geomin, 0.2, 1e-9))).geom as result
FROM a
)
SELECT
avg(st_area(result)) as avg_area
FROM b;