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

Merge branch 'develop' of github.com:CartoDB/crankshaft into add-interpolation

Browse Source
This commit is contained in:
abelvm
2016-09-20 12:57:51 +02:00
parent 2ff20e596e 2261d11de0
commit 5f1cf951ea
107 changed files with 14993 additions and 23 deletions
Download Patch File Download Diff File Expand all files Collapse all files
src/pg/crankshaft.control
+1 -1
View File
@@ -1,5 +1,5 @@
comment = 'CartoDB Spatial Analysis extension'
default_version = '0.2.0'
default_version = '0.4.0'
requires = 'plpythonu, postgis'
superuser = true
schema = cdb_crankshaft
src/pg/sql/09_voronoi.sql
+236
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;
src/pg/sql/13_PIA.sql
+123
View 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;
src/pg/sql/14_densify.sql
+67
View 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;
src/pg/sql/15_tinmap.sql
+43
View 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;
src/pg/sql/19_contour.sql
+205
View File
@@ -0,0 +1,205 @@
CREATE OR REPLACE FUNCTION CDB_Contour(
IN geomin geometry[],
IN colin numeric[],
IN buffer numeric,
IN intmethod integer,
IN classmethod integer,
IN steps integer,
IN max_time integer DEFAULT 60000
)
RETURNS TABLE(
the_geom geometry,
bin integer,
min_value numeric,
max_value numeric,
avg_value numeric
) AS $$
DECLARE
cell_count integer;
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;
-- calc the optimal number of cells for the current dataset
SELECT
CASE intmethod
WHEN 0 THEN round(3.7745903782 * max_time - 9.4399210051 * array_length(geomin,1) - 1350.8778213073)
WHEN 1 THEN round(2.2855592156 * max_time - 87.285217133 * array_length(geomin,1) + 17255.7085601797)
WHEN 2 THEN round(0.9799471999 * max_time - 127.0334085369 * array_length(geomin,1) + 22707.9579721218)
ELSE 10000
END INTO cell_count;
-- 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
),
resolution as(
SELECT
round(|/ (
ST_area(geom) / cell_count
)) as cell
FROM envelope3857
),
grid as(
SELECT
ST_Transform(cdb_crankshaft.CDB_RectangleGrid(e.geom, r.cell, r.cell), 4326) as geom
FROM envelope3857 e, resolution r
),
interp as(
SELECT
geom,
CASE
WHEN intmethod=1 THEN cdb_crankshaft._interp_in_tin(geomin, colin, tin, ST_Centroid(geom))
ELSE cdb_crankshaft.CDB_SpatialInterpolation(geomin, colin, ST_Centroid(geom), intmethod)
END as val
FROM grid
),
classes as(
SELECT CASE
WHEN classmethod = 0 THEN
cdb_crankshaft.CDB_EqualIntervalBins(array_agg(val), steps)
WHEN classmethod = 1 THEN
cdb_crankshaft.CDB_HeadsTailsBins(array_agg(val), steps)
WHEN classmethod = 2 THEN
cdb_crankshaft.CDB_JenksBins(array_agg(val), steps)
ELSE
cdb_crankshaft.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 _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;
src/pg/sql/cdb_utils.sql
+447
View File
@@ -0,0 +1,447 @@
--
-- Fill given extent with a rectangular coverage
--
-- @param ext Extent to fill. Only rectangles with center point falling
-- inside the extent (or at the lower or leftmost edge) will
-- be emitted. The returned hexagons will have the same SRID
-- as this extent.
--
-- @param width With of each rectangle
--
-- @param height Height of each rectangle
--
-- @param origin Optional origin to allow for exact tiling.
-- If omitted the origin will be 0,0.
-- The parameter is checked for having the same SRID
-- as the extent.
--
--
CREATE OR REPLACE FUNCTION CDB_RectangleGrid(ext GEOMETRY, width FLOAT8, height FLOAT8, origin GEOMETRY DEFAULT NULL)
RETURNS SETOF GEOMETRY
AS $$
DECLARE
h GEOMETRY; -- rectangle cell
hstep FLOAT8; -- horizontal step
vstep FLOAT8; -- vertical step
hw FLOAT8; -- half width
hh FLOAT8; -- half height
vstart FLOAT8;
hstart FLOAT8;
hend FLOAT8;
vend FLOAT8;
xoff FLOAT8;
yoff FLOAT8;
xgrd FLOAT8;
ygrd FLOAT8;
x FLOAT8;
y FLOAT8;
srid INTEGER;
BEGIN
srid := ST_SRID(ext);
xoff := 0;
yoff := 0;
IF origin IS NOT NULL THEN
IF ST_SRID(origin) != srid THEN
RAISE EXCEPTION 'SRID mismatch between extent (%) and origin (%)', srid, ST_SRID(origin);
END IF;
xoff := ST_X(origin);
yoff := ST_Y(origin);
END IF;
--RAISE DEBUG 'X offset: %', xoff;
--RAISE DEBUG 'Y offset: %', yoff;
hw := width/2.0;
hh := height/2.0;
xgrd := hw;
ygrd := hh;
--RAISE DEBUG 'X grid size: %', xgrd;
--RAISE DEBUG 'Y grid size: %', ygrd;
hstep := width;
vstep := height;
-- Tweak horizontal start on hstep grid from origin
hstart := xoff + ceil((ST_XMin(ext)-xoff)/hstep)*hstep;
--RAISE DEBUG 'hstart: %', hstart;
-- Tweak vertical start on vstep grid from origin
vstart := yoff + ceil((ST_Ymin(ext)-yoff)/vstep)*vstep;
--RAISE DEBUG 'vstart: %', vstart;
hend := ST_XMax(ext);
vend := ST_YMax(ext);
--RAISE DEBUG 'hend: %', hend;
--RAISE DEBUG 'vend: %', vend;
x := hstart;
WHILE x < hend LOOP -- over X
y := vstart;
h := ST_MakeEnvelope(x-hw, y-hh, x+hw, y+hh, srid);
WHILE y < vend LOOP -- over Y
RETURN NEXT h;
h := ST_Translate(h, 0, vstep);
y := yoff + round(((y + vstep)-yoff)/ygrd)*ygrd; -- round to grid
END LOOP;
x := xoff + round(((x + hstep)-xoff)/xgrd)*xgrd; -- round to grid
END LOOP;
RETURN;
END
$$ LANGUAGE 'plpgsql' IMMUTABLE;
--
-- Calculate the equal interval bins for a given column
--
-- @param in_array A numeric array of numbers to determine the best
-- to determine the bin boundary
--
-- @param breaks The number of bins you want to find.
--
--
-- Returns: upper edges of bins
--
--
CREATE OR REPLACE FUNCTION CDB_EqualIntervalBins ( in_array NUMERIC[], breaks INT ) RETURNS NUMERIC[] as $$
DECLARE
diff numeric;
min_val numeric;
max_val numeric;
tmp_val numeric;
i INT := 1;
reply numeric[];
BEGIN
SELECT min(e), max(e) INTO min_val, max_val FROM ( SELECT unnest(in_array) e ) x WHERE e IS NOT NULL;
diff = (max_val - min_val) / breaks::numeric;
LOOP
IF i < breaks THEN
tmp_val = min_val + i::numeric * diff;
reply = array_append(reply, tmp_val);
i := i+1;
ELSE
reply = array_append(reply, max_val);
EXIT;
END IF;
END LOOP;
RETURN reply;
END;
$$ language plpgsql IMMUTABLE;
--
-- Determine the Heads/Tails classifications from a numeric array
--
-- @param in_array A numeric array of numbers to determine the best
-- bins based on the Heads/Tails method.
--
-- @param breaks The number of bins you want to find.
--
--
CREATE OR REPLACE FUNCTION CDB_HeadsTailsBins ( in_array NUMERIC[], breaks INT) RETURNS NUMERIC[] as $$
DECLARE
element_count INT4;
arr_mean numeric;
i INT := 2;
reply numeric[];
BEGIN
-- get the total size of our row
element_count := array_upper(in_array, 1) - array_lower(in_array, 1);
-- ensure the ordering of in_array
SELECT array_agg(e) INTO in_array FROM (SELECT unnest(in_array) e ORDER BY e) x;
-- stop if no rows
IF element_count IS NULL THEN
RETURN NULL;
END IF;
-- stop if our breaks are more than our input array size
IF element_count < breaks THEN
RETURN in_array;
END IF;
-- get our mean value
SELECT avg(v) INTO arr_mean FROM ( SELECT unnest(in_array) as v ) x;
reply = Array[arr_mean];
-- slice our bread
LOOP
IF i > breaks THEN EXIT; END IF;
SELECT avg(e) INTO arr_mean FROM ( SELECT unnest(in_array) e) x WHERE e > reply[i-1];
IF arr_mean IS NOT NULL THEN
reply = array_append(reply, arr_mean);
END IF;
i := i+1;
END LOOP;
RETURN reply;
END;
$$ language plpgsql IMMUTABLE;
--
-- Determine the Jenks classifications from a numeric array
--
-- @param in_array A numeric array of numbers to determine the best
-- bins based on the Jenks method.
--
-- @param breaks The number of bins you want to find.
--
-- @param iterations The number of different starting positions to test.
--
-- @param invert Optional wheter to return the top of each bin (default)
-- or the bottom. BOOLEAN, default=FALSE.
--
--
CREATE OR REPLACE FUNCTION CDB_JenksBins ( in_array NUMERIC[], breaks INT, iterations INT DEFAULT 5, invert BOOLEAN DEFAULT FALSE) RETURNS NUMERIC[] as $$
DECLARE
element_count INT4;
arr_mean NUMERIC;
bot INT;
top INT;
tops INT[];
classes INT[][];
i INT := 1; j INT := 1;
curr_result NUMERIC[];
best_result NUMERIC[];
seedtarget TEXT;
quant NUMERIC[];
shuffles INT;
BEGIN
-- get the total size of our row
element_count := array_length(in_array, 1); --array_upper(in_array, 1) - array_lower(in_array, 1);
-- ensure the ordering of in_array
SELECT array_agg(e) INTO in_array FROM (SELECT unnest(in_array) e ORDER BY e) x;
-- stop if no rows
IF element_count IS NULL THEN
RETURN NULL;
END IF;
-- stop if our breaks are more than our input array size
IF element_count < breaks THEN
RETURN in_array;
END IF;
shuffles := LEAST(GREATEST(floor(2500000.0/(element_count::float*iterations::float)), 1), 750)::int;
-- get our mean value
SELECT avg(v) INTO arr_mean FROM ( SELECT unnest(in_array) as v ) x;
-- assume best is actually Quantile
SELECT cdb_crankshaft.CDB_QuantileBins(in_array, breaks) INTO quant;
-- if data is very very large, just return quant and be done
IF element_count > 5000000 THEN
RETURN quant;
END IF;
-- change quant into bottom, top markers
LOOP
IF i = 1 THEN
bot = 1;
ELSE
-- use last top to find this bot
bot = top+1;
END IF;
IF i = breaks THEN
top = element_count;
ELSE
SELECT count(*) INTO top FROM ( SELECT unnest(in_array) as v) x WHERE v <= quant[i];
END IF;
IF i = 1 THEN
classes = ARRAY[ARRAY[bot,top]];
ELSE
classes = ARRAY_CAT(classes,ARRAY[bot,top]);
END IF;
IF i > breaks THEN EXIT; END IF;
i = i+1;
END LOOP;
best_result = cdb_crankshaft.CDB_JenksBinsIteration( in_array, breaks, classes, invert, element_count, arr_mean, shuffles);
--set the seed so we can ensure the same results
SELECT setseed(0.4567) INTO seedtarget;
--loop through random starting positions
LOOP
IF j > iterations-1 THEN EXIT; END IF;
i = 1;
tops = ARRAY[element_count];
LOOP
IF i = breaks THEN EXIT; END IF;
SELECT array_agg(distinct e) INTO tops FROM (SELECT unnest(array_cat(tops, ARRAY[floor(random()*element_count::float)::int])) as e ORDER BY e) x WHERE e != 1;
i = array_length(tops, 1);
END LOOP;
i = 1;
LOOP
IF i > breaks THEN EXIT; END IF;
IF i = 1 THEN
bot = 1;
ELSE
bot = top+1;
END IF;
top = tops[i];
IF i = 1 THEN
classes = ARRAY[ARRAY[bot,top]];
ELSE
classes = ARRAY_CAT(classes,ARRAY[bot,top]);
END IF;
i := i+1;
END LOOP;
curr_result = cdb_crankshaft.CDB_JenksBinsIteration( in_array, breaks, classes, invert, element_count, arr_mean, shuffles);
IF curr_result[1] > best_result[1] THEN
best_result = curr_result;
j = j-1; -- if we found a better result, add one more search
END IF;
j = j+1;
END LOOP;
RETURN (best_result)[2:array_upper(best_result, 1)];
END;
$$ language plpgsql IMMUTABLE;
--
-- Perform a single iteration of the Jenks classification
--
CREATE OR REPLACE FUNCTION CDB_JenksBinsIteration ( in_array NUMERIC[], breaks INT, classes INT[][], invert BOOLEAN, element_count INT4, arr_mean NUMERIC, max_search INT DEFAULT 50) RETURNS NUMERIC[] as $$
DECLARE
tmp_val numeric;
new_classes int[][];
tmp_class int[];
i INT := 1;
j INT := 1;
side INT := 2;
sdam numeric;
gvf numeric := 0.0;
new_gvf numeric;
arr_gvf numeric[];
class_avg numeric;
class_max_i INT;
class_min_i INT;
class_max numeric;
class_min numeric;
reply numeric[];
BEGIN
-- Calculate the sum of squared deviations from the array mean (SDAM).
SELECT sum((arr_mean - e)^2) INTO sdam FROM ( SELECT unnest(in_array) as e ) x;
--Identify the breaks for the lowest GVF
LOOP
i = 1;
LOOP
-- get our mean
SELECT avg(e) INTO class_avg FROM ( SELECT unnest(in_array[classes[i][1]:classes[i][2]]) as e) x;
-- find the deviation
SELECT sum((class_avg-e)^2) INTO tmp_val FROM ( SELECT unnest(in_array[classes[i][1]:classes[i][2]]) as e ) x;
IF i = 1 THEN
arr_gvf = ARRAY[tmp_val];
-- init our min/max map for later
class_max = arr_gvf[i];
class_min = arr_gvf[i];
class_min_i = 1;
class_max_i = 1;
ELSE
arr_gvf = array_append(arr_gvf, tmp_val);
END IF;
i := i+1;
IF i > breaks THEN EXIT; END IF;
END LOOP;
-- calculate our new GVF
SELECT sdam-sum(e) INTO new_gvf FROM ( SELECT unnest(arr_gvf) as e ) x;
-- if no improvement was made, exit
IF new_gvf < gvf THEN EXIT; END IF;
gvf = new_gvf;
IF j > max_search THEN EXIT; END IF;
j = j+1;
i = 1;
LOOP
--establish directionality (uppward through classes or downward)
IF arr_gvf[i] < class_min THEN
class_min = arr_gvf[i];
class_min_i = i;
END IF;
IF arr_gvf[i] > class_max THEN
class_max = arr_gvf[i];
class_max_i = i;
END IF;
i := i+1;
IF i > breaks THEN EXIT; END IF;
END LOOP;
IF class_max_i > class_min_i THEN
class_min_i = class_max_i - 1;
ELSE
class_min_i = class_max_i + 1;
END IF;
--Move from higher class to a lower gid order
IF class_max_i > class_min_i THEN
classes[class_max_i][1] = classes[class_max_i][1] + 1;
classes[class_min_i][2] = classes[class_min_i][2] + 1;
ELSE -- Move from lower class UP into a higher class by gid
classes[class_max_i][2] = classes[class_max_i][2] - 1;
classes[class_min_i][1] = classes[class_min_i][1] - 1;
END IF;
END LOOP;
i = 1;
LOOP
IF invert = TRUE THEN
side = 1; --default returns bottom side of breaks, invert returns top side
END IF;
reply = array_append(reply, in_array[classes[i][side]]);
i = i+1;
IF i > breaks THEN EXIT; END IF;
END LOOP;
RETURN array_prepend(gvf, reply);
END;
$$ language plpgsql IMMUTABLE;
--
-- Determine the Quantile classifications from a numeric array
--
-- @param in_array A numeric array of numbers to determine the best
-- bins based on the Quantile method.
--
-- @param breaks The number of bins you want to find.
--
--
CREATE OR REPLACE FUNCTION CDB_QuantileBins ( in_array NUMERIC[], breaks INT) RETURNS NUMERIC[] as $$
DECLARE
element_count INT4;
break_size numeric;
tmp_val numeric;
i INT := 1;
reply numeric[];
BEGIN
-- sort our values
SELECT array_agg(e) INTO in_array FROM (SELECT unnest(in_array) e ORDER BY e ASC) x;
-- get the total size of our data
element_count := array_length(in_array, 1);
break_size := element_count::numeric / breaks;
-- slice our bread
LOOP
IF i < breaks THEN
IF break_size * i % 1 > 0 THEN
SELECT e INTO tmp_val FROM ( SELECT unnest(in_array) e LIMIT 1 OFFSET ceil(break_size * i) - 1) x;
ELSE
SELECT avg(e) INTO tmp_val FROM ( SELECT unnest(in_array) e LIMIT 2 OFFSET ceil(break_size * i) - 1 ) x;
END IF;
ELSIF i = breaks THEN
-- select the last value
SELECT max(e) INTO tmp_val FROM ( SELECT unnest(in_array) e ) x;
ELSE
EXIT;
END IF;
reply = array_append(reply, tmp_val);
i := i+1;
END LOOP;
RETURN reply;
END;
$$ language plpgsql IMMUTABLE;
src/pg/test/expected/09_voronoi_test.out
+7
View File
@@ -0,0 +1,7 @@
SET client_min_messages TO WARNING;
\set ECHO none
avg_area
----------------------
0.000178661700690617
(1 row)
src/pg/test/expected/13_pia_test.out
+7
View File
@@ -0,0 +1,7 @@
SET client_min_messages TO WARNING;
\set ECHO none
st_astext
-------------------------------------------
POINT(-3.67484492582767 40.4395084885993)
(1 row)
src/pg/test/expected/14_densify_test.out
+50
View File
@@ -0,0 +1,50 @@
SET client_min_messages TO WARNING;
\set ECHO none
cdb_densify
-----------------------------------------------------------------
(01010000001361C3D32B6501403255302AA9B34440,7.0)
(01010000002497FF907EFB0040F085C954C1B04440,8.0)
(0101000000A167B3EA73350140E4141DC9E5AF4440,1.0)
(010100000062A1D634EF38014014D044D8F0B44440,2.0)
(010100000052B81E85EB510140EEEBC03923B24440,3.0)
(0101000000C286A757CA3201409D8026C286AF4440,5.0)
(01010000007DD0B359F5390140F38E537424AF4440,6.0)
(0101000000D237691A140D0140014EEFE2FDB44440,4.0)
(01010000003312B4DCAC14014047F8F1AAE1B14440,4.3333333333333333)
(010100000048C0FBBD27290140A9BBC5D646B34440,2.3333333333333333)
(01010000001DEBE2361A400140F79A0B4953B24440,2.0000000000000000)
(010100000098933D2F02500140115B676994B34440,4.0000000000000000)
(01010000004BA3DC90595001405C456C9DA5B14440,5.3333333333333333)
(0101000000D1FA8198714001404285107D64B04440,3.3333333333333333)
(01010000004AEA043411360140D2B687AA85AF4440,4.0000000000000000)
(0101000000CCA4736BBF2201402B8716D9CEAF4440,6.3333333333333333)
(0101000000832C1EF13E2101402609AF4A0FB04440,4.6666666666666667)
(010100000063A009D8BF090140BEEE38F68AB24440,5.4444444444444444)
(010100000019AE5D3CF8180140A5008D2162B34440,3.5555555555555555)
(01010000007E88BE590E250140EB9DA83067B44440,2.7777777777777778)
(0101000000C05105B65D3B014044A2D0B2EEB34440,2.7777777777777778)
(010100000005329D125F4F01401C8049790FB44440,4.3333333333333333)
(010100000054E45F2DB3570140BBDE724420B34440,4.6666666666666667)
(0101000000E53EEA4DD05701407FD7C9557BB24440,5.1111111111111111)
(01010000004A9CC694D34F01402B63A5BE7BB14440,6.1111111111111111)
(0101000000DD24068195430140317345DA64B04440,4.8888888888888889)
(010100000033E768B7D23A014002994E89AFAF4440,4.4444444444444444)
(010100000083C0CAA145360140CAEC55A065AF4440,5.0000000000000000)
(010100000004549A09D52F01403E3230057EAF4440,5.7777777777777778)
(010100000024040D72661D0140B0DEBBE0E6AF4440,6.7777777777777778)
(01010000007CCD85A42915014017B22F2835B04440,6.3333333333333333)
(0101000000F5F145CA781001401F2DCE18E6B04440,5.6666666666666667)
(0101000000F10DE756572701402134E4146CB14440,3.6666666666666667)
(01010000008894DB45FA290140F8C4EB987EB24440,2.8888888888888889)
(0101000000AABF5E61C13D0140E6E512830FB34440,2.7777777777777778)
(0101000000581215F9574B0140FDF5BB4EAEB24440,3.0000000000000000)
(0101000000EA6C9F19754B0140C0EE126009B24440,3.4444444444444444)
(01010000001383C0CAA145014088CC827674B14440,3.5555555555555555)
(0101000000D0B02B40EE350140C90D9905EDB04440,3.3333333333333333)
(010100000051DA1B7C613201406F36F4DA1DB04440,3.0000000000000000)
(0101000000EA6E133D52390140FD1AE7FAEFAF4440,2.7777777777777778)
(01010000003A487527C5340140C76EEE11A6AF4440,3.3333333333333333)
(0101000000BADB448F542E01403AB4C876BEAF4440,4.1111111111111111)
(0101000000FB12176D7B280140BDE1A0F9EBAF4440,4.0000000000000000)
(44 rows)
src/pg/test/expected/15_tinmap_test.out
+83
View File
@@ -0,0 +1,83 @@
SET client_min_messages TO WARNING;
\set ECHO none
cdb_tinmap
---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
(01030000000100000004000000D237691A140D0140014EEFE2FDB444402497FF907EFB0040F085C954C1B0444063A009D8BF090140BEEE38F68AB24440D237691A140D0140014EEFE2FDB44440,5.8148148148148148)
(01030000000100000004000000D237691A140D0140014EEFE2FDB4444063A009D8BF090140BEEE38F68AB2444019AE5D3CF8180140A5008D2162B34440D237691A140D0140014EEFE2FDB44440,4.3333333333333333)
(01030000000100000004000000D237691A140D0140014EEFE2FDB4444019AE5D3CF8180140A5008D2162B344407E88BE590E250140EB9DA83067B44440D237691A140D0140014EEFE2FDB44440,3.4444444444444444)
(01030000000100000004000000D237691A140D0140014EEFE2FDB444407E88BE590E250140EB9DA83067B4444062A1D634EF38014014D044D8F0B44440D237691A140D0140014EEFE2FDB44440,2.9259259259259259)
(0103000000010000000400000062A1D634EF38014014D044D8F0B444407E88BE590E250140EB9DA83067B44440C05105B65D3B014044A2D0B2EEB3444062A1D634EF38014014D044D8F0B44440,2.5185185185185185)
(0103000000010000000400000062A1D634EF38014014D044D8F0B44440C05105B65D3B014044A2D0B2EEB3444005329D125F4F01401C8049790FB4444062A1D634EF38014014D044D8F0B44440,3.0370370370370370)
(0103000000010000000400000062A1D634EF38014014D044D8F0B4444005329D125F4F01401C8049790FB444401361C3D32B6501403255302AA9B3444062A1D634EF38014014D044D8F0B44440,4.4444444444444444)
(010300000001000000040000001361C3D32B6501403255302AA9B3444005329D125F4F01401C8049790FB4444098933D2F02500140115B676994B344401361C3D32B6501403255302AA9B34440,5.1111111111111111)
(010300000001000000040000001361C3D32B6501403255302AA9B3444098933D2F02500140115B676994B3444054E45F2DB3570140BBDE724420B344401361C3D32B6501403255302AA9B34440,5.2222222222222222)
(010300000001000000040000001361C3D32B6501403255302AA9B3444054E45F2DB3570140BBDE724420B34440E53EEA4DD05701407FD7C9557BB244401361C3D32B6501403255302AA9B34440,5.5925925925925926)
(010300000001000000040000001361C3D32B6501403255302AA9B34440E53EEA4DD05701407FD7C9557BB244404A9CC694D34F01402B63A5BE7BB144401361C3D32B6501403255302AA9B34440,6.0740740740740741)
(010300000001000000040000007DD0B359F5390140F38E537424AF44404A9CC694D34F01402B63A5BE7BB14440DD24068195430140317345DA64B044407DD0B359F5390140F38E537424AF4440,5.6666666666666667)
(010300000001000000040000007DD0B359F5390140F38E537424AF4440DD24068195430140317345DA64B0444033E768B7D23A014002994E89AFAF44407DD0B359F5390140F38E537424AF4440,5.1111111111111111)
(010300000001000000040000007DD0B359F5390140F38E537424AF444033E768B7D23A014002994E89AFAF44404AEA043411360140D2B687AA85AF44407DD0B359F5390140F38E537424AF4440,4.8148148148148148)
(010300000001000000040000007DD0B359F5390140F38E537424AF44404AEA043411360140D2B687AA85AF4440C286A757CA3201409D8026C286AF44407DD0B359F5390140F38E537424AF4440,5.0000000000000000)
(010300000001000000040000007DD0B359F5390140F38E537424AF4440C286A757CA3201409D8026C286AF444004549A09D52F01403E3230057EAF44407DD0B359F5390140F38E537424AF4440,5.5925925925925926)
(010300000001000000040000007DD0B359F5390140F38E537424AF444004549A09D52F01403E3230057EAF4440CCA4736BBF2201402B8716D9CEAF44407DD0B359F5390140F38E537424AF4440,6.0370370370370370)
(010300000001000000040000007DD0B359F5390140F38E537424AF4440CCA4736BBF2201402B8716D9CEAF444024040D72661D0140B0DEBBE0E6AF44407DD0B359F5390140F38E537424AF4440,6.3703703703703704)
(0103000000010000000400000024040D72661D0140B0DEBBE0E6AF4440CCA4736BBF2201402B8716D9CEAF4440832C1EF13E2101402609AF4A0FB0444024040D72661D0140B0DEBBE0E6AF4440,5.9259259259259259)
(0103000000010000000400000024040D72661D0140B0DEBBE0E6AF4440832C1EF13E2101402609AF4A0FB044407CCD85A42915014017B22F2835B0444024040D72661D0140B0DEBBE0E6AF4440,5.9259259259259259)
(0103000000010000000400000024040D72661D0140B0DEBBE0E6AF44407CCD85A42915014017B22F2835B044402497FF907EFB0040F085C954C1B0444024040D72661D0140B0DEBBE0E6AF4440,7.0370370370370370)
(010300000001000000040000002497FF907EFB0040F085C954C1B044407CCD85A42915014017B22F2835B04440F5F145CA781001401F2DCE18E6B044402497FF907EFB0040F085C954C1B04440,6.6666666666666667)
(010300000001000000040000002497FF907EFB0040F085C954C1B04440F5F145CA781001401F2DCE18E6B0444063A009D8BF090140BEEE38F68AB244402497FF907EFB0040F085C954C1B04440,6.3703703703703704)
(0103000000010000000400000063A009D8BF090140BEEE38F68AB24440F5F145CA781001401F2DCE18E6B044403312B4DCAC14014047F8F1AAE1B1444063A009D8BF090140BEEE38F68AB24440,5.1481481481481481)
(0103000000010000000400000063A009D8BF090140BEEE38F68AB244403312B4DCAC14014047F8F1AAE1B1444019AE5D3CF8180140A5008D2162B3444063A009D8BF090140BEEE38F68AB24440,4.4444444444444444)
(0103000000010000000400000019AE5D3CF8180140A5008D2162B344403312B4DCAC14014047F8F1AAE1B144408894DB45FA290140F8C4EB987EB2444019AE5D3CF8180140A5008D2162B34440,3.5925925925925926)
(0103000000010000000400000019AE5D3CF8180140A5008D2162B344408894DB45FA290140F8C4EB987EB2444048C0FBBD27290140A9BBC5D646B3444019AE5D3CF8180140A5008D2162B34440,2.9259259259259259)
(0103000000010000000400000019AE5D3CF8180140A5008D2162B3444048C0FBBD27290140A9BBC5D646B344407E88BE590E250140EB9DA83067B4444019AE5D3CF8180140A5008D2162B34440,2.8888888888888889)
(010300000001000000040000007E88BE590E250140EB9DA83067B4444048C0FBBD27290140A9BBC5D646B34440C05105B65D3B014044A2D0B2EEB344407E88BE590E250140EB9DA83067B44440,2.6296296296296296)
(01030000000100000004000000C05105B65D3B014044A2D0B2EEB3444048C0FBBD27290140A9BBC5D646B34440AABF5E61C13D0140E6E512830FB34440C05105B65D3B014044A2D0B2EEB34440,2.6296296296296296)
(01030000000100000004000000C05105B65D3B014044A2D0B2EEB34440AABF5E61C13D0140E6E512830FB3444098933D2F02500140115B676994B34440C05105B65D3B014044A2D0B2EEB34440,3.1851851851851852)
(01030000000100000004000000C05105B65D3B014044A2D0B2EEB3444098933D2F02500140115B676994B3444005329D125F4F01401C8049790FB44440C05105B65D3B014044A2D0B2EEB34440,3.7037037037037037)
(0103000000010000000400000098933D2F02500140115B676994B34440AABF5E61C13D0140E6E512830FB34440581215F9574B0140FDF5BB4EAEB2444098933D2F02500140115B676994B34440,3.2592592592592593)
(0103000000010000000400000098933D2F02500140115B676994B34440581215F9574B0140FDF5BB4EAEB2444054E45F2DB3570140BBDE724420B3444098933D2F02500140115B676994B34440,3.8888888888888889)
(0103000000010000000400000054E45F2DB3570140BBDE724420B34440581215F9574B0140FDF5BB4EAEB24440E53EEA4DD05701407FD7C9557BB2444054E45F2DB3570140BBDE724420B34440,4.2592592592592593)
(01030000000100000004000000E53EEA4DD05701407FD7C9557BB24440581215F9574B0140FDF5BB4EAEB2444052B81E85EB510140EEEBC03923B24440E53EEA4DD05701407FD7C9557BB24440,3.7037037037037037)
(01030000000100000004000000E53EEA4DD05701407FD7C9557BB2444052B81E85EB510140EEEBC03923B244404BA3DC90595001405C456C9DA5B14440E53EEA4DD05701407FD7C9557BB24440,4.4814814814814815)
(01030000000100000004000000E53EEA4DD05701407FD7C9557BB244404BA3DC90595001405C456C9DA5B144404A9CC694D34F01402B63A5BE7BB14440E53EEA4DD05701407FD7C9557BB24440,5.5185185185185185)
(010300000001000000040000004A9CC694D34F01402B63A5BE7BB144404BA3DC90595001405C456C9DA5B144401383C0CAA145014088CC827674B144404A9CC694D34F01402B63A5BE7BB14440,5.0000000000000000)
(010300000001000000040000004A9CC694D34F01402B63A5BE7BB144401383C0CAA145014088CC827674B14440DD24068195430140317345DA64B044404A9CC694D34F01402B63A5BE7BB14440,4.8518518518518518)
(01030000000100000004000000DD24068195430140317345DA64B044401383C0CAA145014088CC827674B14440D1FA8198714001404285107D64B04440DD24068195430140317345DA64B04440,3.9259259259259259)
(01030000000100000004000000DD24068195430140317345DA64B04440D1FA8198714001404285107D64B0444033E768B7D23A014002994E89AFAF4440DD24068195430140317345DA64B04440,4.2222222222222222)
(0103000000010000000400000033E768B7D23A014002994E89AFAF4440D1FA8198714001404285107D64B04440EA6E133D52390140FD1AE7FAEFAF444033E768B7D23A014002994E89AFAF4440,3.5185185185185185)
(0103000000010000000400000033E768B7D23A014002994E89AFAF4440EA6E133D52390140FD1AE7FAEFAF4440A167B3EA73350140E4141DC9E5AF444033E768B7D23A014002994E89AFAF4440,2.7407407407407407)
(0103000000010000000400000033E768B7D23A014002994E89AFAF4440A167B3EA73350140E4141DC9E5AF44403A487527C5340140C76EEE11A6AF444033E768B7D23A014002994E89AFAF4440,2.9259259259259259)
(0103000000010000000400000033E768B7D23A014002994E89AFAF44403A487527C5340140C76EEE11A6AF44404AEA043411360140D2B687AA85AF444033E768B7D23A014002994E89AFAF4440,3.9259259259259259)
(010300000001000000040000004AEA043411360140D2B687AA85AF44403A487527C5340140C76EEE11A6AF4440C286A757CA3201409D8026C286AF44404AEA043411360140D2B687AA85AF4440,4.1111111111111111)
(01030000000100000004000000C286A757CA3201409D8026C286AF44403A487527C5340140C76EEE11A6AF4440BADB448F542E01403AB4C876BEAF4440C286A757CA3201409D8026C286AF4440,4.1481481481481481)
(01030000000100000004000000C286A757CA3201409D8026C286AF4440BADB448F542E01403AB4C876BEAF444004549A09D52F01403E3230057EAF4440C286A757CA3201409D8026C286AF4440,4.9629629629629630)
(0103000000010000000400000004549A09D52F01403E3230057EAF4440BADB448F542E01403AB4C876BEAF4440FB12176D7B280140BDE1A0F9EBAF444004549A09D52F01403E3230057EAF4440,4.6296296296296296)
(0103000000010000000400000004549A09D52F01403E3230057EAF4440FB12176D7B280140BDE1A0F9EBAF4440CCA4736BBF2201402B8716D9CEAF444004549A09D52F01403E3230057EAF4440,5.3703703703703704)
(01030000000100000004000000CCA4736BBF2201402B8716D9CEAF4440FB12176D7B280140BDE1A0F9EBAF4440832C1EF13E2101402609AF4A0FB04440CCA4736BBF2201402B8716D9CEAF4440,5.0000000000000000)
(01030000000100000004000000832C1EF13E2101402609AF4A0FB04440FB12176D7B280140BDE1A0F9EBAF4440F10DE756572701402134E4146CB14440832C1EF13E2101402609AF4A0FB04440,4.1111111111111111)
(01030000000100000004000000832C1EF13E2101402609AF4A0FB04440F10DE756572701402134E4146CB14440F5F145CA781001401F2DCE18E6B04440832C1EF13E2101402609AF4A0FB04440,4.6666666666666667)
(01030000000100000004000000832C1EF13E2101402609AF4A0FB04440F5F145CA781001401F2DCE18E6B044407CCD85A42915014017B22F2835B04440832C1EF13E2101402609AF4A0FB04440,5.5555555555555556)
(01030000000100000004000000F5F145CA781001401F2DCE18E6B04440F10DE756572701402134E4146CB144403312B4DCAC14014047F8F1AAE1B14440F5F145CA781001401F2DCE18E6B04440,4.5555555555555556)
(010300000001000000040000003312B4DCAC14014047F8F1AAE1B14440F10DE756572701402134E4146CB144408894DB45FA290140F8C4EB987EB244403312B4DCAC14014047F8F1AAE1B14440,3.6296296296296296)
(010300000001000000040000008894DB45FA290140F8C4EB987EB24440F10DE756572701402134E4146CB144401DEBE2361A400140F79A0B4953B244408894DB45FA290140F8C4EB987EB24440,2.8518518518518519)
(010300000001000000040000008894DB45FA290140F8C4EB987EB244401DEBE2361A400140F79A0B4953B24440AABF5E61C13D0140E6E512830FB344408894DB45FA290140F8C4EB987EB24440,2.5555555555555556)
(010300000001000000040000008894DB45FA290140F8C4EB987EB24440AABF5E61C13D0140E6E512830FB3444048C0FBBD27290140A9BBC5D646B344408894DB45FA290140F8C4EB987EB24440,2.6666666666666667)
(01030000000100000004000000AABF5E61C13D0140E6E512830FB344401DEBE2361A400140F79A0B4953B24440581215F9574B0140FDF5BB4EAEB24440AABF5E61C13D0140E6E512830FB34440,2.5925925925925926)
(01030000000100000004000000581215F9574B0140FDF5BB4EAEB244401DEBE2361A400140F79A0B4953B24440EA6C9F19754B0140C0EE126009B24440581215F9574B0140FDF5BB4EAEB24440,2.8148148148148148)
(01030000000100000004000000581215F9574B0140FDF5BB4EAEB24440EA6C9F19754B0140C0EE126009B2444052B81E85EB510140EEEBC03923B24440581215F9574B0140FDF5BB4EAEB24440,3.1481481481481481)
(0103000000010000000400000052B81E85EB510140EEEBC03923B24440EA6C9F19754B0140C0EE126009B244404BA3DC90595001405C456C9DA5B1444052B81E85EB510140EEEBC03923B24440,3.9259259259259259)
(010300000001000000040000004BA3DC90595001405C456C9DA5B14440EA6C9F19754B0140C0EE126009B244401383C0CAA145014088CC827674B144404BA3DC90595001405C456C9DA5B14440,4.1111111111111111)
(010300000001000000040000001383C0CAA145014088CC827674B14440EA6C9F19754B0140C0EE126009B244401DEBE2361A400140F79A0B4953B244401383C0CAA145014088CC827674B14440,3.0000000000000000)
(010300000001000000040000001383C0CAA145014088CC827674B144401DEBE2361A400140F79A0B4953B24440D0B02B40EE350140C90D9905EDB044401383C0CAA145014088CC827674B14440,2.9629629629629629)
(010300000001000000040000001383C0CAA145014088CC827674B14440D0B02B40EE350140C90D9905EDB04440D1FA8198714001404285107D64B044401383C0CAA145014088CC827674B14440,3.4074074074074074)
(01030000000100000004000000D1FA8198714001404285107D64B04440D0B02B40EE350140C90D9905EDB0444051DA1B7C613201406F36F4DA1DB04440D1FA8198714001404285107D64B04440,3.2222222222222222)
(01030000000100000004000000D1FA8198714001404285107D64B0444051DA1B7C613201406F36F4DA1DB04440EA6E133D52390140FD1AE7FAEFAF4440D1FA8198714001404285107D64B04440,3.0370370370370370)
(01030000000100000004000000EA6E133D52390140FD1AE7FAEFAF444051DA1B7C613201406F36F4DA1DB04440A167B3EA73350140E4141DC9E5AF4440EA6E133D52390140FD1AE7FAEFAF4440,2.2592592592592593)
(01030000000100000004000000A167B3EA73350140E4141DC9E5AF444051DA1B7C613201406F36F4DA1DB04440BADB448F542E01403AB4C876BEAF4440A167B3EA73350140E4141DC9E5AF4440,2.7037037037037037)
(01030000000100000004000000A167B3EA73350140E4141DC9E5AF4440BADB448F542E01403AB4C876BEAF44403A487527C5340140C76EEE11A6AF4440A167B3EA73350140E4141DC9E5AF4440,2.8148148148148148)
(01030000000100000004000000BADB448F542E01403AB4C876BEAF444051DA1B7C613201406F36F4DA1DB04440FB12176D7B280140BDE1A0F9EBAF4440BADB448F542E01403AB4C876BEAF4440,3.7037037037037037)
(01030000000100000004000000FB12176D7B280140BDE1A0F9EBAF444051DA1B7C613201406F36F4DA1DB04440F10DE756572701402134E4146CB14440FB12176D7B280140BDE1A0F9EBAF4440,3.5555555555555556)
(01030000000100000004000000F10DE756572701402134E4146CB1444051DA1B7C613201406F36F4DA1DB04440D0B02B40EE350140C90D9905EDB04440F10DE756572701402134E4146CB14440,3.3333333333333333)
(01030000000100000004000000F10DE756572701402134E4146CB14440D0B02B40EE350140C90D9905EDB044401DEBE2361A400140F79A0B4953B24440F10DE756572701402134E4146CB14440,3.0000000000000000)
(77 rows)
src/pg/test/expected/19_contour_test.out
+9
View File
@@ -0,0 +1,9 @@
SET client_min_messages TO WARNING;
\set ECHO none
bin|avg_value
0|280.23070673030491816424178
1|413.81702914846213479025305
2|479.6334491374486884098328
3|529.1545236882183479447113
4|614.1132081424930103122037
(5 rows)
src/pg/test/sql/09_voronoi_test.sql
+15
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;
src/pg/test/sql/13_pia_test.sql
+7
View File
@@ -0,0 +1,7 @@
SET client_min_messages TO WARNING;
\set ECHO none
with a as(
select st_geomfromtext('POLYGON((-432540.453078056 4949775.20452642,-432329.947920966 4951361.232584,-431245.028163694 4952223.31516671,-429131.071033529 4951768.00415574,-424622.07505895 4952843.13503987,-423688.327170174 4953499.20752423,-424086.294349759 4954968.38274191,-423068.388925945 4954378.63345336,-423387.653225542 4953355.67417084,-420594.869840519 4953781.00230592,-416026.095299382 4951484.06849063,-412483.018546414 4951024.5410983,-410490.399661215 4954502.24032205,-408186.197521284 4956398.91417441,-407627.262358013 4959300.94633864,-406948.770061627 4959874.85407739,-404949.583326472 4959047.74518163,-402570.908447199 4953743.46829807,-400971.358683991 4952193.11680804,-403533.488084088 4949649.89857885,-406335.177028373 4950193.19571096,-407790.456731515 4952391.46015616,-412060.672398345 4950381.2389307,-410716.93482498 4949156.7509561,-408464.162289794 4943912.8940387,-409350.599394983 4942819.84896006,-408087.791091424 4942451.6711778,-407274.045613725 4940572.4807777,-404446.196589102 4939976.71501489,-402422.964843936 4940450.3670813,-401010.654464241 4939054.8061663,-397647.247369412 4940679.80737878,-395658.413346901 4940528.84765185,-395536.852462953 4938829.79565997,-394268.923462818 4938003.7277717,-393388.720249116 4934757.80596815,-392393.301362444 4934326.71675815,-392573.527618037 4932323.40974412,-393464.640141837 4931903.10653605,-393085.597275686 4931094.7353605,-398426.261165985 4929156.87541607,-398261.174361137 4926238.00816416,-394045.059966834 4925765.18668498,-392982.960705174 4926391.81893628,-393090.272694301 4927176.84692181,-391648.240010564 4924626.06386961,-391889.914625075 4923086.14787613,-394345.177314013 4923235.086036,-395550.878718795 4917812.79243978,-399009.463978251 4912927.7157945,-398948.794855767 4911941.91010796,-398092.636652078 4911806.57392519,-401991.601817112 4911722.9204501,-406225.972607907 4914505.47286319,-411104.994569885 4912569.26941163,-412925.513522316 4913030.3608866,-414630.148884835 4914436.69169949,-414207.691417276 4919205.78028405,-418306.141109809 4917994.9580478,-424184.700779621 4918938.12432889,-426816.961458921 4923664.37379373,-420956.324227126 4923381.98014807,-420186.661267781 4924286.48693378,-420943.411166194 4926812.76394433,-419779.45457046 4928527.43466337,-419768.767899344 4930681.94459216,-421911.668097113 4930432.40620397,-423482.386112205 4933451.28047252,-427272.814773717 4934151.56473242,-427144.908678797 4939731.77191996,-428982.125554848 4940522.84445172,-428986.133056516 4942437.17281266,-431237.792396792 4947309.68284815,-432476.889648814 4947791.74800037,-432540.453078056 4949775.20452642))', 3857) as g
)
SELECT st_astext(cdb_crankshaft.CDB_PIA(g)) from a;
src/pg/test/sql/14_densify_test.sql
+9
View File
@@ -0,0 +1,9 @@
SET client_min_messages TO WARNING;
\set ECHO none
with data as (
select
ARRAY[7.0,8.0,1.0,2.0,3.0,5.0,6.0,4.0] as colin,
ARRAY[ST_GeomFromText('POINT(2.1744 41.4036)'),ST_GeomFromText('POINT(2.1228 41.3809)'),ST_GeomFromText('POINT(2.1511 41.3742)'),ST_GeomFromText('POINT(2.1528 41.4136)'),ST_GeomFromText('POINT(2.165 41.3917)'),ST_GeomFromText('POINT(2.1498 41.3713)'),ST_GeomFromText('POINT(2.1533 41.3683)'),ST_GeomFromText('POINT(2.131386 41.413998)')] as geomin
)
select cdb_crankshaft.CDB_Densify(geomin, colin, 2) from data;
src/pg/test/sql/15_tinmap_test.sql
+9
View File
@@ -0,0 +1,9 @@
SET client_min_messages TO WARNING;
\set ECHO none
with data as (
select
ARRAY[7.0,8.0,1.0,2.0,3.0,5.0,6.0,4.0] as colin,
ARRAY[ST_GeomFromText('POINT(2.1744 41.4036)'),ST_GeomFromText('POINT(2.1228 41.3809)'),ST_GeomFromText('POINT(2.1511 41.3742)'),ST_GeomFromText('POINT(2.1528 41.4136)'),ST_GeomFromText('POINT(2.165 41.3917)'),ST_GeomFromText('POINT(2.1498 41.3713)'),ST_GeomFromText('POINT(2.1533 41.3683)'),ST_GeomFromText('POINT(2.131386 41.413998)')] as geomin
)
select cdb_crankshaft.CDB_TINmap(geomin, colin, 2) from data;
src/pg/test/sql/19_contour_test.sql
+17
View File
@@ -0,0 +1,17 @@
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)',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 g
),
b as(
SELECT
foo.*
FROM
a,
cdb_crankshaft.CDB_contour(a.g, a.vals, 0.0, 1, 3, 5, 60) foo
)
SELECT bin, avg_value from b order by bin;