Dit is de code:
drop table if exists daken_eindhoven;
with dump as
(
select gid, ST_ZMin( geometrie ) floorz,
(st_dump(geometrie)).geom::geometry(polygonz,7415) as geom3d
from bag3d_eindhoven
),
data as
(
select gid,
floorz,
st_3dSlope(geom3d)::integer as slope,
st_3dAzimuth(geom3d)::integer as azimuth,
st_3dHeight(geom3d) as height,
ST_ZMin( geom3d ) minz,
ST_ZMax( geom3d ) maxz,
--geom3d,
st_transform( st_force2d( geom3d ), 28992)::geometry(polygon,28992) geom2d
from dump
)
select *
into daken_eindhoven
from data
where data.minz > data.floorz+2 -- anders krijg je ook het vloeropp
;
En die query maakt gebruikt de volgende functies:
DROP FUNCTION st_3DSlope(geometry);
DROP FUNCTION st_3DAzimuth(geometry);
DROP FUNCTION st_3DHeight(geometry);
DROP FUNCTION ST_3DPointOnSurface(geometry);
DROP FUNCTION st_3Dinfo(geometry) ;
--******************************
CREATE OR REPLACE FUNCTION public.st_3dheight( geom3d geometry ) RETURNS float LANGUAGE plpgsql
AS $function$
begin
return abs( st_zmax(geom3d)-st_zmin(geom3d) );
END; $function$;
--******************************
CREATE OR REPLACE FUNCTION public.st_3dslope( geom3d geometry ) RETURNS float LANGUAGE plpgsql
AS $function$
declare
info float[2];
begin
info := st_3dinfo( geom3d );
return info[1];
END; $function$;
--******************************
CREATE OR REPLACE FUNCTION public.st_3dazimuth( geom3d geometry ) RETURNS float LANGUAGE plpgsql
AS $function$
declare
info float[2];
begin
info := st_3dinfo( geom3d );
return info[2];
END; $function$;
--******************************
CREATE OR REPLACE FUNCTION public.ST_3DPointOnSurface( geom3d geometry )
RETURNS geometry LANGUAGE plpgsql
AS $function$
declare
z_max float;
z_min float;
geom2d geometry;
begin
z_max = ST_ZMax( geom3d );
z_min = ST_ZMin( geom3d );
geom2d := st_force2d( geom3d );
return ST_3DClosestPoint(
geom3d,
st_force3D(
ST_PointOnSurface( geom2d ),
(z_max+z_min)/2)
);
END; $function$;
--******************************
CREATE OR REPLACE FUNCTION public.st_3dinfo( geom3d geometry )
RETURNS float[][]
LANGUAGE plpgsql
AS $function$
declare
midpoint geometry;
normalvector geometry;
geomstart geometry;
geomend geometry;
h float;
d float;
slope float;
azimuth float;
p1 geometry;
p2 geometry;
p3 geometry;
try integer;
isok boolean;
begin
isok := true;
slope := 0;
azimuth := null;
--check for 3d
if not ST_NDims( geom3d ) = 3 then
isok := false;
end if;
--check for minimal slope (make it a default parameter?)
if ST_3DHeight( geom3d) < 0.1 then
isok := false;
end if;
--do the complicated stuff
if isok then
--**************************************
--create normal vector
midpoint := ST_3DPointOnSurface( geom3d );
midpoint := st_force3D(st_force2d(midpoint), st_z(midpoint)+0.1); --raise point with 10cm
normalvector := ST_3DShortestLine( midpoint, geom3d);
--**************************************
--get azimuth = direction of line
geomstart := ST_StartPoint(normalvector);
geomend := ST_EndPoint(normalvector);
azimuth := degrees(ST_Azimuth( geomend,geomstart )) ;
--**************************************
--get slope = angle between two 3D points
h := abs( st_z( geomstart ) - st_z( geomend ) );
d := st_distance( st_force2d(geomstart), st_force2d(geomend) );
p1 := st_makePoint( 0,0 );
p2 := st_makePoint( 0,h );
p3 := st_makePoint( d,0 );
slope := 90-degrees(st_angle(p1,p3,p2));
end if;
return array[ slope, azimuth ];
--EXCEPTION WHEN others THEN
-- raise notice 'ST_slope3d error';
-- return false;
END; $function$
;
--**************************************************
