Skip to content

Script 03: Surface Attributes

File: sql/scripts/main/03_calc_child_feat_attr.sql

Reads from: {city2tabula_schema}.{lod_schema}_child_feature_geom_dump

Writes to: {city2tabula_schema}.{lod_schema}_surface_raw

Purpose

For each individual polygon face (from script 02), this script computes these physical attributes:

  • Tilt: the surface's angle from vertical (0° = vertical wall, 90° = flat horizontal roof). This is the opposite of the usual from-horizontal slope convention; see the box below.
  • Azimuth: the compass direction the surface faces (0°/360° = North, 90° = East, 180° = South, 270° = West).
  • Surface area: the true 3D face area in square metres.
  • Height span: the vertical range of the face (Z max − Z min) in metres.
  • Horizontal and slope extent: the face's extent along the horizontal and along the line of steepest slope in its own plane, in metres.

It then marks the faces that lie between two solids of the same building as internal (see Internal faces).

Every attribute except height span is derived from the surface normal: a vector that points perpendicularly outward from the face. The bulk of this script's CTEs are dedicated to computing that normal correctly.

Tilt convention: 0° = wall, 90° = flat roof

tilt is DEGREES(ASIN(ABS(nz))), where nz is the vertical component of the unit surface normal. It is the angle between the surface normal and the horizontal plane, or equivalently the angle between the surface and the vertical. A vertical wall is 0°, a flat roof or floor is 90°, a 45° pitched roof is 45°.

Most building-energy calculations (including what ignis and BuEM expect) use the opposite reference: the surface tilt measured from the horizontal, where a flat roof is 0° and a wall is 90°. The two are complements, so convert with 90 − tilt. They coincide only at 45°. Whoever feeds c2t surface output to ignis or BuEM must apply the conversion. The same note is on the Surface type in internal/onrequest/query.go and in docs/api.md.

Background: what is a surface normal?

The surface normal is a stick set perpendicular to a face and pointing away from the building. It encodes two things:

  • Its vertical component (how much it points up or down) determines tilt.
  • Its horizontal direction (which compass bearing it points toward) determines azimuth.

For the normal to be useful, it must:

  1. Be correctly computed from the polygon's vertices.
  2. Point outward (away from the building), not inward.

Both of these requirements require care when working with real-world CityGML data, which is why script 03 has eight CTEs.

Background: Newell's method

Why not a simple cross product to compute the normal vector?

The simplest way to compute a surface normal is to take three polygon vertices and compute their cross product. For a flat (planar) polygon this works, but CityGML surfaces are sometimes slightly non-planar due to data quality issues (vertices don't lie perfectly on a plane). Picking an arbitrary triplet on a warped surface gives an inconsistent result.

Newell's method

Newell's method solves this by accumulating a contribution from every consecutive edge pair around the polygon ring:

n_x = SUM over all edges: (y_i − y_{i+1}) × (z_i + z_{i+1})
n_y = SUM over all edges: (z_i − z_{i+1}) × (x_i + x_{i+1})
n_z = SUM over all edges: (x_i − x_{i+1}) × (y_i + y_{i+1})

For a perfectly flat polygon this gives the same answer as a simple cross product. For a warped polygon it gives the best-fit normal that minimises the error across all vertices: equivalent to fitting a plane by least squares. It runs in O(n) time (one pass over the edges) and handles both cases uniformly.

Background: vertex winding order

The sign of a cross-product normal depends on the order in which vertices are listed: tracing the polygon boundary clockwise points the normal one way, counter-clockwise the other. CityGML datasets from different providers may use either convention, and sometimes even mix them. Without correction, a north-facing wall might be computed as south-facing (180° error) purely because its vertices happen to be listed in a different order.

This script corrects for this by comparing the computed normal against a known reference point inside the solid the face bounds (the interior point, derived from that solid's GroundSurface). If the normal points toward the interior rather than away from it, it is flipped.

CTE walkthrough

Step 1: new_buildings

SELECT DISTINCT building_feature_id
FROM {city2tabula_schema}.{lod_schema}_child_feature_geom_dump
WHERE building_feature_id IN {building_ids}
  AND building_feature_id NOT IN (
    SELECT DISTINCT building_feature_id
    FROM {city2tabula_schema}.{lod_schema}_surface_raw
  )

Idempotency guard: selects only buildings not yet present in the output table.

Step 2: owner_interior_pts

SELECT
  owner_feature_id,
  ST_PointOnSurface(ST_Collect(ST_Force2D(geom))) AS interior_pt
FROM ...
WHERE classname = 'GroundSurface'
GROUP BY owner_feature_id

For each solid (the Building or one of its BuildingParts), computes a single 2D reference point guaranteed to lie inside its ground footprint. It is taken per solid because, in a building of several parts, a point inside one part would flip the outward test for walls of another.

  • ST_Force2D strips Z coordinates, since only the horizontal position is needed.
  • ST_Collect merges all GroundSurface polygons into one geometry without dissolving them.
  • ST_PointOnSurface returns a point on the surface that is guaranteed to be inside, even for non-convex shapes like L-shaped footprints (where the centroid can fall outside).

This point is used in oriented_normals to check whether a wall's computed normal points outward or inward.

Step 3: raw_surfaces

SELECT
  ...,
  gd.geom AS valid_geom,
  ST_IsPlanar(gd.geom) AS is_planar
FROM {city2tabula_schema}.{lod_schema}_child_feature_geom_dump gd
INNER JOIN new_buildings nb ON gd.building_feature_id = nb.building_feature_id

One row per polygon face from the geometry dump. ST_IsPlanar checks whether all vertices lie on a single plane. This flag is stored in the output for diagnostic purposes but does not change how the normal is computed: Newell's method handles both planar and non-planar polygons identically.

Step 4: surface_points

SELECT
  *,
  (ST_DumpPoints(valid_geom)).geom AS point_geom,
  (ST_DumpPoints(valid_geom)).path[2] AS pt_idx
FROM raw_surfaces

ST_DumpPoints extracts every individual vertex from the polygon ring as a separate row. path[2] gives the sequential index of each vertex within the ring, so they can be walked in order.

Example: a quadrilateral face with vertices v1, v2, v3, v4 (plus v1 again to close the ring) produces five rows with pt_idx values 1, 2, 3, 4, 5.

Step 5: surface_edges

SELECT
  ...,
  point_geom,
  LEAD(point_geom) OVER (PARTITION BY id ORDER BY pt_idx) AS next_pt
FROM surface_points

LEAD is a window function that pairs each vertex row with the next vertex in the sequence. This turns the list of vertices into a list of directed edges (current_vertex → next_vertex). The last vertex in the ring is paired with NULL (there is no next vertex), which is filtered out in the next step.

Result: each row now represents one edge, described by its start and end point, ready for Newell's formula.

Step 6: surface_normals

SELECT
  ...,
  SUM((ST_Y(point_geom) - ST_Y(next_pt)) * (ST_Z(point_geom) + ST_Z(next_pt))) AS n_x,
  SUM((ST_Z(point_geom) - ST_Z(next_pt)) * (ST_X(point_geom) + ST_X(next_pt))) AS n_y,
  SUM((ST_X(point_geom) - ST_X(next_pt)) * (ST_Y(point_geom) + ST_Y(next_pt))) AS n_z
FROM surface_edges
WHERE next_pt IS NOT NULL
GROUP BY id, ...
HAVING COUNT(*) >= 2

Applies Newell's formula across all edges of each polygon face. The three sums give the raw (un-normalised) normal vector (n_x, n_y, n_z).

The outer query then computes the vector's magnitude (sqrt(n_x² + n_y² + n_z²)) and discards surfaces where the magnitude is near zero, which means all vertices are collinear (i.e. the polygon has no area and no well-defined normal).

Step 7: oriented_normals

COALESCE(
  CASE
    WHEN (centroid_x - interior_x) * (n_x / magnitude)
       + (centroid_y - interior_y) * (n_y / magnitude)
    >= 0 THEN 1.0
    ELSE -1.0
  END,
  1.0
) AS surface_flip

For each WallSurface and RoofSurface, computes a dot product between: - The horizontal component of the surface normal (nx, ny), and - The vector from the solid's interior point to the surface centroid.

If the dot product is ≥ 0, the normal already points away from the interior (outward): surface_flip = +1, keep as-is. If the dot product is negative, the normal points inward: surface_flip = −1, flip it.

If no interior point exists (the building has no GroundSurface), COALESCE falls back to +1 (no flip).

Step 8: normalized_normals

CASE
  WHEN classname = 'RoofSurface' AND (n_z / magnitude) < 0
    THEN -(n_x / magnitude)
  WHEN classname = 'WallSurface'
    THEN surface_flip * (n_x / magnitude)
  ELSE n_x / magnitude
END AS nx

Applies per-class rules to produce the final unit normal (nx, ny, nz):

  • RoofSurface: if the Z component is negative (normal points downward), flip all three components. A roof normal must always point upward. The Z sign check is used here rather than the dot-product flip because, for near-flat roofs, the horizontal components (nx, ny) are tiny and sign-unstable.
  • WallSurface: multiply all three components by surface_flip from the previous step. This corrects the CW/CCW winding issue.
  • Other types: divide by magnitude to get a unit vector, no flip.

Final SELECT: computing the attributes

With the unit normal (nx, ny, nz) available, the attributes are:

-- Tilt: angle of the normal above the horizontal plane
CASE
  WHEN ABS(nz) > 0.985 THEN 90.0  -- near-horizontal; snapped to flat
  ELSE DEGREES(ASIN(ABS(nz)))
END AS tilt
tilt diagram tilt diagram
Figure 1: Tilt and normal vector components

nz is the vertical component of the unit normal. For a vertical wall, nz = 0, giving ASIN(0) = 0°. For a flat horizontal surface, |nz| = 1, giving ASIN(1) = 90°. For |nz| > 0.985 (tilt > ~80°), tilt is snapped to exactly 90° rather than the raw computed value. This is the same near-horizontal band that azimuth treats as undefined, and for the same reason: the residual angle there is noise rather than signal.

-- Azimuth: compass bearing of the horizontal normal component
CASE
  WHEN ABS(nz) > 0.985 THEN -1  -- surface is near-horizontal; azimuth is undefined
  ELSE MOD((450.0 - degrees(atan2(ny, nx))) + 360.0, 360.0)
END AS azimuth
azimuth diagram azimuth diagram
Figure 2: Azimuth and normal vector components
azimuth diagram azimuth diagram
Figure 3: Azimuth undefined for a near-horizontal surface

atan2(ny, nx) computes the mathematical angle (counter-clockwise from East). The formula (450 − angle) mod 360 converts this to a compass bearing (clockwise from North). For nearly-flat surfaces (tilt > ~80°), the horizontal components are near zero and atan2 becomes numerically unreliable, so azimuth is set to -1 (undefined).

-- Surface area: projected onto the surface plane before measurement
{city2tabula_schema}.surface_area_corrected_geom(valid_geom, nx, ny, nz)

ST_Area(geom) measures area on the XY plane, which underestimates a tilted surface (the same way a tilted roof looks smaller from directly above than its true size). The surface_area_corrected_geom function projects the polygon onto its own surface plane before measuring, giving the true 3D area.

-- Height: vertical span of the face
(ST_ZMax(valid_geom) - ST_ZMin(valid_geom)) AS height

The simplest attribute: the difference between the highest and lowest Z coordinate of the polygon's vertices.

-- Horizontal and slope extent: projections on two axes in the surface plane
LEFT JOIN LATERAL {city2tabula_schema}.surface_dimensions(valid_geom, nx, ny, nz) d
  ON objectclass_id IN (709, 710, 712)

surface_dimensions fixes two axes in the face's plane from its unit normal n = (nx, ny, nz), with h = √(nx² + ny²):

Axis Direction Column
u (−ny, nx, 0) / h: horizontal, along the face horizontal_extent
v (−nz·nx / h, −nz·ny / h, h): the line of steepest slope, pointing up slope_extent

Each extent is the largest minus the smallest projection of the face's vertices on its axis, so the sign of the normal does not change it. That matters for walls, whose outward flip in oriented_normals points some walls of concave footprints inward. Faces with |nz| > 0.985, where azimuth is −1, have no horizontal axis of their own and get NULL in both columns.

Face horizontal_extent slope_extent height
Wall 10 m wide, 3 m tall 10 3 3
Wall 2 m wide, 8 m tall 2 8 8
Roof, 12 m eave, 30° pitch, 5 m deep in plan 12 5.77 2.89
Gable wall, 8 m base, 3 m to the apex 8 3 3
Flat roof or ground face NULL NULL 0

On a planar face height = slope_extent × cos(tilt). For a wall, horizontal_extent is the run that adds up across coplanar wall segments. For triangles, L-shapes and other non-rectangular faces, horizontal_extent × slope_extent is larger than the face; surface_area / (horizontal_extent × slope_extent) gives the fraction of that rectangle the face fills.

Internal faces

Where two solids of one building touch, the face region between them is interior. A source that models each part as a closed solid (Vienna) has that face on both parts; a source that leaves the parts open on that side (Prague) has it on neither. A second statement after the INSERT finds these faces for the batch.

A face of one solid is internal where a face of another solid of the same building lies on it:

Test Rule
Classes Wall against wall, or roof against ground
Orientation Normals parallel within about 2.6° (\|cos\| ≥ 0.999); the sign of the normals is not used
Plane offset ≤ 0.05 m
Overlap 2D intersection of the two faces in the plane of the first

face_to_plane (in sql/functions/02_face_plane.sql) rotates both faces so the first face's normal points along +Z. The plane offset is then a Z difference and the overlap a 2D intersection. face_from_plane rotates the exposed remainder back into world coordinates.

The sign of the normals is not used because the outward flip in oriented_normals points some walls of concave footprints inward. Two solids cannot both have exterior wall on the same patch, so the overlap is interior whichever way the stored normals point.

The exposed remainder keeps only pieces wider than 0.04 m and larger than 0.01 m², so the tolerance band leaves no slivers. For each face with an overlap of at least 0.01 m²:

  • area_internal = surface_area minus the exposed area.
  • geom_exposed = the exposed pieces, or an empty geometry when the whole face is internal.

Faces with no overlap keep both columns NULL. Faces between different buildings (party walls) are not tested here.

Output columns

Column Description
owner_feature_id Feature id of the solid the face bounds: the Building or one of its BuildingParts
surface_area True 3D area of the face (sqm)
tilt Angle of the surface from vertical (degrees; 0° = wall, 90° = flat roof, snapped to 90° for \|nz\| > 0.985). Complement of the usual from-horizontal slope angle: 90 − tilt.
azimuth Compass bearing of outward-facing normal (degrees; −1 = undefined)
is_valid ST_IsValid result for the polygon
is_planar Whether all vertices lie on a single plane
height Vertical span: ZMax − ZMin (m)
horizontal_extent Extent along the horizontal in the surface plane (m); NULL for faces with \|nz\| > 0.985 and for classes other than wall, roof and ground
slope_extent Extent along the line of steepest slope (m); NULL as for horizontal_extent
normal_x, normal_y, normal_z Unit normal after the orientation rules above
area_internal Part of surface_area that lies against another solid of the same building (sqm); NULL when none
geom_exposed Exterior remainder of an internal face; empty when fully internal; NULL when the face has no internal part
geom Original polygon geometry (carried through)

What comes next

Script 04 reads _surface_raw and aggregates the per-face attributes into one row per solid and one summary row per building.