Ik heb de BAG3D van TU Delft ingelezen in een postgreSQL PostGIS database. Nu probeer ik de helling en de richting (slope en azimuth) van de (dak) vlakken te determineren. Ik kan helaas geen postgis functie vinden die dit snel en gemakkelijk teruggeeft per 3Dvlak.
Ik heb nu een eigen st_slope3D en st_azimuth3D functie gemaakt die de punten van een vlak dumpt en op basis van het laagste en hoogste punt (st_minz en st_maxz) wat uitrekent, Maar dat gaat nog niet helemaal goed, en is ook wat sloom.
Daarom toch hier ook nog even gevraagd: Zijn er snellere en gemakkelijkere manieren om de slope en de azimuth uit te rekenen van een 3D polygoon?
Ik weet zo gauw niet een PostGIS-functie hiervoor, maar volgens mij is je uitgangspunt met minz en maxz (dz) niet juist. Je moet denk ik de steilste lijn vinden, dus die waarvan de verhouding afstand-xy tot dz het grootst is. Maar goed, dat bedenk ik bij mijn eerste kopje koffie op maandagochtend. 
Voor de snelheid kun je misschien 1 functie maken die zowel slope als azimuth berekent, omdat de functies nu 2x bijna hetzelfde doen. Dat scheelt weer een beetje.
Nog een bakkie verder bedenk ik dat het toch nog iets lastiger is. De steilste lijn bij een driehoek (bijvoorbeeld zijvlak piramide) loopt niet naar 1 van de onderste punten maar naar het midden van het de onderste zijde…
Volgens mij staat de uitwerking hier: mathematics - Calculate the slope of a 3D triangle? - Game Development Stack Exchange
En misschien dat CGAL dit kan, die library wordt ook gebruikt in postgresql: https://www.cgal.org/
In mijn beleving bereken je eerst de normaalvector van een vlak: Meetkunde in 3D
Maar eerlijk gezegd denk ik dat je niet zo low-level hoeft te gaan rekenen. Dit is vast eerder gedaan .
Je ziet dat ook in de 3d BAG viewer gebeuren als je op een dakvlak klikt.
Deze link is nog interessant: 3DDistancecalc – PostGIS
Die normaal vector is het denk ik wel. Zit even te zoeken tussen de postgis 3D functies. Ik denk dat de simpelste manier om de normaal vector te maken de functie ST_3DShortestLine() is.
- Maak 3d-punt boven het 3d-vlak = st_makepoint() met st_centroid( 2d vlak) [x,y] en ST_ZMax() [z] (+1?)
- normaalvector = ST_3DShortestLine( 3d-punt, 3d-vlak)
En met die lijn kan ik dan de slope en de azimuth berekenen. Ik weet niet of ik er vandaag aan toekom, maar ik zal melden wat de uitkomst is.
Enorm bedankt in ieder geval voor het meedenken!
Ik denk niet dat dat altijd goed gaat, omdat je soms de afstand naar een vertex of edge zult berekenen en niet naar het vlak. Ook moet je rekening houden met de richting van de normal, omdat een vlak 2 zijden heeft.
Volgens mij staat het scriptje gewoon in de uitwerking die ik eerder stuurde en kun je dat namaken in postgres (of als python functie?)
Ter info: Het resultaat: Moet nog even goed testen of ergens nog gekke dingen gebeuren.
Nogmaals dank voor de input!
4 likes
Zou je in de beste open source gedacht ook willen delen hoe je dit hebt gedaan? Dan heeft een volgende met dezelfde vraag er ook wat aan!
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$
;
--**************************************************
3 likes
@ArendjanvdNeut Bedankt voor het delen! Ik probeer het dakoppervlak en de hellingshoek te bepalen voor een adres. In de BAG 3d viewer is de hoek te zien:
Echter, hoe krijg ik deze uit de resultaten van je functies? Zou dat de Azimuth moeten zijn?
Daarnaast, welke LOD gebruik je? Want ik zie dat je bag3d_eindhoven gebruikt, is dat een custom set op basis van LOD22_3d?
Azimuth is de richting en de Slope is de hellingshoek. Dus deze functies:
st_3dSlope(geom3d) en st_3dAzimuth(geom3d)
Mijn bag3d_eindhoven was een testtukje (zelf gemaakte selectie) uit een grotere set die ik gedownload heb. Wat ik indertijd exact gedownload heb weet ik niet meer, maar je hebt 3D (dak) geometrien nodig.
