Created
October 1, 2012 10:01
-
-
Save ppKrauss/3810651 to your computer and use it in GitHub Desktop.
PostGIS shape descriptor, "as rectangle" metrics vector
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
CREATE FUNCTION shapedescr_sizes( | |
-- Shape-descriptor "as rectangle" for geometry description by sizes. | |
gbase geometry, -- input | |
-- p_seqs integer DEFAULT 8, -- deprecated? for st_buffer(g,w,p_seqs) or point-buffer inference | |
-- p_shape varchar DEFAULT '', -- will be endcap indicator | |
p_decplacesof_zero integer DEFAULT 6, -- precision of zero when rounding delta | |
p_dwmin float DEFAULT 99999999.0, -- change to ex. 0.0001, if to use. | |
p_deltaimpact float DEFAULT 9999.0 -- internal (maximized by probability of negative delta) | |
) RETURNS float[] AS $f$ | |
DECLARE | |
ret float[]; | |
dw float; | |
b float; | |
L_estim float; | |
H_estim float; | |
aorig float; | |
gaux geometry; | |
g1 geometry; | |
A0 float; | |
A1 float; | |
c float; | |
delta float; | |
per float; | |
errcod float; | |
BEGIN | |
errcod=0.0; | |
IF gbase IS NULL OR NOT(ST_IsClosed(gbase)) THEN | |
errcod=1; -- ERROR1 (die) | |
RAISE EXCEPTION 'error %: invalid input geometry',errcod; | |
END IF; | |
A0 := ST_Area(gbase); | |
per := st_perimeter(gbase); | |
dw := sqrt(A0)/p_deltaimpact; | |
IF dw>p_dwmin THEN dw:=p_dwmin; END IF; | |
g1 = ST_Buffer(gbase,dw); | |
A1 = ST_area(g1); | |
IF A0>A1 THEN | |
errcod=10; -- ERROR2 (die) | |
RAISE EXCEPTION 'error %: invalid buffer/geometry with A0=% g.t. A1=%',errcod,A0,A1; | |
END IF; | |
IF (A1-A0)>1.001*dw*per THEN | |
gaux := ST_Buffer(g1,-dw); -- closing operation. | |
A0 = ST_Area(gaux); -- changed area | |
per := ST_Perimeter(gaux); -- changed | |
errcod:=errcod + 0.1; -- Warning3 | |
END IF; | |
C := 2.0*dw; | |
b := -(A1-A0)/C+C; | |
delta := b^2-4.0*A0; | |
IF delta<0.0 AND round(delta,p_decplacesof_zero)<=0.0 THEN | |
delta=0.0; -- for regular shapes like the square | |
errcod:=errcod + 0.01; -- Warning2 | |
END IF; | |
IF delta<0.0 THEN | |
L_estim := NULL; | |
H_estim := NULL; | |
errcod:=errcod+100; -- ERROR3 | |
ELSE | |
L_estim := (-b + sqrt(delta))/2.0; | |
H_estim := (-b - sqrt(delta))/2.0; | |
END IF; | |
IF abs(A0-L_estim*H_estim)>0.001 THEN | |
errcod:=errcod + 0.001; -- Warning1 | |
END IF; | |
ret := array[L_estim,H_estim,a0,per,dw,errcod]; | |
return ret; | |
END | |
$f$ LANGUAGE plpgsql IMMUTABLE; | |
CREATE or replace FUNCTION shapedescr_sizes_tr( | |
-- Translator for shapedescr_sizes(). Uses ROUND(float). | |
float[], -- shapedescr_sizes() returned vector | |
integer DEFAULT 0, -- general round parameter | |
integer DEFAULT 3 -- number of "decimal warnings" | |
) RETURNS varchar[] -- length, height, area, perimeter, dw, radius, err_message | |
AS $f$ | |
SELECT array[ | |
round(L,$2+1)::varchar, round(H,$2+1)::varchar, round(area,$2)::varchar, | |
round(perim,$2+1)::varchar, round(dw,$2+3)::varchar, | |
round(sqrt(L*H/pi()),$2+1)::varchar, -- radius for "shape as circle" | |
CASE WHEN errcod>3.0 THEN 'ERROR ' ||round(errcod-$3) | |
WHEN errcod>0.0 THEN 'WARNING '||round(errcod)||CASE | |
WHEN round(10^(-errcod),$3-1)!=($1)[6] THEN '.'|| ($1)[6]*10^$3 | |
ELSE '' | |
END | |
ELSE '' | |
END::varchar | |
] | |
FROM ( | |
SELECT ($1)[1] as L, ($1)[2] as H, ($1)[3] as area, ($1)[4] as perim, | |
($1)[5] as dw, log(($1)[6]*10.0^($3+1)+1.0) as errcod | |
) as t; | |
$f$ LANGUAGE SQL IMMUTABLE; | |
CREATE FUNCTION shapedescr_sizes_tr(geometry, integer DEFAULT 0, integer DEFAULT 3) | |
RETURNS varchar[] AS $f$ | |
SELECT shapedescr_sizes_tr(shapedescr_sizes($1),$2,$3); | |
$f$ LANGUAGE SQL IMMUTABLE; | |
-- visual check: | |
SELECT shapedescr_sizes_tr(ST_MakePolygon(ST_GeomFromText('LINESTRING(75 30,77 29,78 29.5, 75 30)'))); -- 0.6 | |
-- simulate areas and its central lines, and calculates estimatives. | |
SELECT tipo, round("L",1) as "L", 2*w as "H", | |
round(descr[1],1) as "L_estim", round(descr[2],1) as "H_estim", | |
round(ST_Area(g0)/"L",1) as "HbyL_Area", round(0.5*ST_Perimeter(g0)-"L",1) as "HbyL_perim" | |
FROM ( | |
SELECT *, st_length(g) as "L", shapedescr_sizes(g0) AS descr | |
FROM ( | |
SELECT *, st_buffer(g,w,tipo) as g0 | |
FROM ( | |
SELECT ST_GeomFromText('LINESTRING(50 50,150 150,160 180,200 200)') AS g, | |
unnest(array[2.0,10.0,20.0,50.0]) AS w | |
) AS t, ( | |
SELECT unnest(array['','endcap=square','endcap=flat join=bevel']) AS tipo | |
) AS t2 | |
) as t3 | |
) as t4; | |
-- tests with database: | |
--SELECT shapedescr_sizes_tr(the_geom) as descr FROM quadras LIMIT 100; |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
-- -- -- -- -- -- -- -- -- -- -- | |
-- -- -- -- -- -- -- -- -- -- -- | |
-- ADDING TESTS OF SECOND METHOD | |
CREATE or replace FUNCTION shapedescr_sizes_test( | |
-- TESTING grow of precision: result NO. | |
-- Shape-descriptor "as rectangle" for geometry description by sizes. | |
gbase geometry, -- input | |
-- p_seqs integer DEFAULT 8, -- deprecated? for st_buffer(g,w,p_seqs) or point-buffer inference | |
-- p_shape varchar DEFAULT '', -- will be endcap indicator | |
p_decplacesof_zero integer DEFAULT 4, -- precision of zero when rounding delta | |
p_dwmin float DEFAULT 99999999.0, -- change to ex. 0.0001, if to use. | |
p_deltaimpact float DEFAULT 9999.0 -- internal (maximized by probability of negative delta) | |
) RETURNS float[] AS $f$ | |
DECLARE | |
ret float[]; | |
dw float; | |
b float; | |
L_estim float; | |
H_estim float; | |
aorig float; | |
gaux geometry; | |
g1 geometry; | |
A0 float; | |
A1 float; | |
A2 float; | |
delta float; | |
per float; | |
errcod float; | |
BEGIN | |
errcod=0.0; | |
IF gbase IS NULL OR NOT(ST_IsClosed(gbase)) THEN | |
errcod=1; -- ERROR1 (die) | |
RAISE EXCEPTION 'error %: invalid input geometry',errcod; | |
END IF; | |
A0 := ST_Area(gbase); | |
per := st_perimeter(gbase); -- only for dw*per cheching, and for output | |
dw := sqrt(A0)/p_deltaimpact; | |
IF dw>p_dwmin THEN dw:=p_dwmin; END IF; | |
g1 = ST_Buffer(gbase,dw); | |
A1 = ST_area(g1); | |
IF A0>A1 THEN | |
errcod=10; -- ERROR2 (die) | |
RAISE EXCEPTION 'error %: invalid buffer/geometry with A0=% g.t. A1=%',errcod,A0,A1; | |
END IF; | |
IF (A1-A0)>1.001*dw*per THEN | |
gaux := ST_Buffer(g1,-dw); -- closing operation. | |
A0 = ST_Area(gaux); -- changed area | |
per := ST_Perimeter(gaux); -- changed | |
errcod:=errcod + 0.1; -- Warning3 | |
A2 = ST_Area(ST_Buffer(gaux,-dw)); -- second changed area | |
ELSE | |
A2 = ST_Area(ST_Buffer(gbase,-dw)); | |
END IF; | |
-- C := 2.0*dw; | |
b := -(A1-A2)/(4.0*dw); -- -(A1-A0)/C+C; | |
delta := b^2-4.0*A0; | |
IF delta<0.0 AND round(delta,p_decplacesof_zero)<=0.0 THEN | |
delta=0.0; -- for regular shapes like the square | |
errcod:=errcod + 0.01; -- Warning2 | |
END IF; | |
IF delta<0.0 THEN | |
L_estim := NULL; | |
H_estim := NULL; | |
errcod:=errcod+100; -- ERROR3 | |
ELSE | |
L_estim := (-b + sqrt(delta))/2.0; | |
H_estim := (-b - sqrt(delta))/2.0; | |
END IF; | |
IF abs(A0-L_estim*H_estim)>0.001 THEN | |
errcod:=errcod + 0.001; -- Warning1 | |
END IF; | |
ret := array[L_estim,H_estim,a0,per,dw,errcod]; | |
return ret; | |
END | |
$f$ LANGUAGE plpgsql IMMUTABLE; | |
-- COMPARING METHODS | |
SELECT tipo, round("L",1) as "L", 2*w as "H", | |
round(descr[1],3) as "L_estim", round(descr[2],3) as "H_estim", | |
round(descr2[1],3) as "L_estim_v2", round(descr2[2],3) as "H_estim_v2" | |
FROM ( | |
SELECT *, st_length(g) as "L", shapedescr_sizes(g0) AS descr, | |
shapedescr_sizes_test(g0) AS descr2 | |
FROM ( | |
SELECT *, st_buffer(g,w,tipo) as g0 | |
FROM ( | |
SELECT ST_GeomFromText('LINESTRING(50 50,150 150,160 180,200 200)') AS g, | |
unnest(array[2.0,10.0,20.0,50.0]) AS w, 0.005 as dw | |
) AS t, ( | |
SELECT unnest(array['','endcap=square','endcap=flat join=bevel']) AS tipo | |
) AS t2 | |
) as t3 | |
) as t4; -- result = very little gain, NOT justify heavy shapedescr_sizes_test(). |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment