add files
This commit is contained in:
abelvm
parent
2fe50e1f71
commit
55f8626901
35
doc/09_voronoi.md
Normal file
35
doc/09_voronoi.md
Normal file
@ -0,0 +1,35 @@
|
||||
## Spacial interpolation
|
||||
|
||||
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 | geometry | The target point to calc the value |
|
||||
|
||||
### 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)'),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 geomin
|
||||
)
|
||||
SELECT
|
||||
CDB_voronoi(geomin, 0.2, 1e-9) as result
|
||||
FROM a;
|
||||
```
|
||||
205
src/pg/sql/09_voronoi.sql
Normal file
205
src/pg/sql/09_voronoi.sql
Normal file
@ -0,0 +1,205 @@
|
||||
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
|
||||
WITH
|
||||
convexhull_1 as (
|
||||
SELECT
|
||||
ST_ConvexHull(ST_Collect(geomin)) as g,
|
||||
buffer * |/ st_area(ST_ConvexHull(ST_Collect(geomin))) 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_Force_2D(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
|
||||
)
|
||||
SELECT
|
||||
st_collect(ST_intersection(c.g, v.g))
|
||||
-- ST_intersection(c.g, v.g)
|
||||
INTO geomout
|
||||
FROM
|
||||
voro_set v,
|
||||
clipper c;
|
||||
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
|
||||
* @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.
|
||||
* @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(p_pt1 geometry, p_pt2 geometry, 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_pt1) <> 'ST_Point' OR
|
||||
ST_GeometryType(p_pt1) <> '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;
|
||||
5
src/pg/test/expected/09_voronoi_test.out
Normal file
5
src/pg/test/expected/09_voronoi_test.out
Normal file
File diff suppressed because one or more lines are too long
7
src/pg/test/sql/09_voronoi_test.sql
Normal file
7
src/pg/test/sql/09_voronoi_test.sql
Normal file
@ -0,0 +1,7 @@
|
||||
WITH a AS (
|
||||
SELECT
|
||||
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 geomin
|
||||
)
|
||||
SELECT
|
||||
CDB_voronoi(geomin, 0.2, 1e-9) as result
|
||||
FROM a;
|
||||
Loading…
Reference in New Issue
Block a user