A drone survey of a roof ends with a point cloud: millions of 3D points that together describe the surface. For solar design and roofing estimates, the numbers people actually need are much simpler. How steep is each roof plane? Which direction does it face? How much surface area does it cover?
In this tutorial we will compute all three with about 60 lines of Python and NumPy. We will fit planes to points, use RANSAC to separate roof planes from each other, and convert the results into the units roofers and solar designers use: degrees, rise over 12 inches, compass bearing and square meters.
What we are computing
Every flat roof plane can be described by its normal vector, the unit vector pointing straight out of the surface. Once you have it, the rest is geometry:
• Tilt (pitch angle): the angle between the normal and straight up, arccos(n_z).
• Rise per 12: the pitch roofers use, 12 * tan(tilt). A 6/12 roof rises 6 inches for every 12 inches of horizontal run.
• Azimuth: the compass direction the plane faces, from the horizontal part of the normal, atan2(n_x, n_y). Zero is north, 90 is east, 180 is south.
• Surface area: the plane's footprint seen from above, divided by cos(tilt).
This assumes your coordinates are in a projected system in meters (or feet), with x pointing east, y north and z up. Pitch and azimuth do not depend on the unit, as long as x, y and z share it. Area comes out in the square of whatever unit you use.
Setup
You need NumPy and SciPy:
pip install numpy scipy
If you want to load real survey data later, add laspy too.
Step 1: Fit a plane to a set of points
Given the points on one roof plane, the best-fit plane passes through their centroid, and its normal is the direction in which the points vary the least. Singular value decomposition (SVD) gives us that directly.
import numpy as np
from scipy.spatial import ConvexHull
def fit_plane(points):
"""Least-squares plane through Nx3 points. Returns (unit normal, centroid)."""
centroid = points.mean(axis=0)
# The normal is the direction of least variance (smallest singular vector)
_, _, vt = np.linalg.svd(points - centroid, full_matrices=False)
normal = vt[-1]
if normal[2] < 0: # make the normal point upward
normal = -normal
return normal, centroid
Flipping the normal when z is negative matters. A plane has two valid normals, and without this step the azimuth could point 180 degrees the wrong way.
Step 2: Find the dominant plane with RANSAC
A real roof is not one plane. It has several, plus chimneys, vents, noise and stray points. A plain least-squares fit over everything would produce a meaningless average.
RANSAC solves this. It repeatedly picks three random points, builds a plane through them, and counts how many other points lie within a small distance of it. The plane with the most supporters wins.
def ransac_plane(points, threshold=0.05, iterations=500, rng=None):
"""Find the dominant plane. Returns a boolean inlier mask."""
rng = rng or np.random.default_rng(42)
best_mask = np.zeros(len(points), dtype=bool)
for _ in range(iterations):
sample = points[rng.choice(len(points), 3, replace=False)]
n = np.cross(sample[1] - sample[0], sample[2] - sample[0])
norm = np.linalg.norm(n)
if norm < 1e-9: # the three points are collinear
continue
n /= norm
dist = np.abs((points - sample[0]) @ n)
mask = dist < threshold
if mask.sum() > best_mask.sum():
best_mask = mask
return best_mask
The threshold is the distance, in your coordinate units, within which a point counts as lying on the plane. With meters, 0.05 means 5 centimeters. Choose it based on your data's noise. Too small and real roof points get rejected; too large and neighboring planes merge.
Step 3: Turn a plane into roofing numbers
Now we refit using only the inliers and compute the metrics.
def plane_metrics(points):
"""Pitch, azimuth and surface area for one roof plane."""
normal, _ = fit_plane(points)
tilt = np.degrees(np.arccos(np.clip(normal[2], -1, 1)))
rise_per_12 = 12 * np.tan(np.radians(tilt))
# Compass bearing of the downslope direction (0 = north, 90 = east)
azimuth = np.degrees(np.arctan2(normal[0], normal[1])) % 360
# Footprint area divided by cos(tilt) gives sloped surface area
footprint = ConvexHull(points[:, :2]).volume # in 2D, .volume is the area
area = footprint / np.cos(np.radians(tilt))
return {"tilt_deg": tilt, "rise_per_12": rise_per_12,
"azimuth_deg": azimuth, "area_m2": area}
One detail that trips people up: for a 2D convex hull, SciPy's .volume attribute is the area, and .area is the perimeter.
Step 4: Peel off one plane at a time
To handle a roof with several planes, run RANSAC, remove the points it found, and repeat on what is left.
def segment_roof(points, max_planes=4, threshold=0.05, min_points=200):
"""Sequential RANSAC: peel off one plane at a time."""
remaining = points.copy()
planes = []
for _ in range(max_planes):
if len(remaining) < min_points:
break
mask = ransac_plane(remaining, threshold)
if mask.sum() < min_points:
break
planes.append(remaining[mask])
remaining = remaining[~mask]
return planes
Try it on a synthetic roof
Before touching real data, test on something where you know the answer. Here is a 12 by 10 meter gable roof with a ridge running east to west and a 6/12 pitch, which is a tilt of about 26.57 degrees. We add 1 centimeter of noise.
rng = np.random.default_rng(1)
n = 4000
x = rng.uniform(0, 12, n)
y = rng.uniform(0, 10, n)
ridge_height = 6.0
z = ridge_height - 0.5 * np.abs(y - 5) + rng.normal(0, 0.01, n)
points = np.column_stack([x, y, z])
for i, plane in enumerate(segment_roof(points, max_planes=2), 1):
m = plane_metrics(plane)
print(f"Plane {i}: tilt {m['tilt_deg']:.1f} deg | "
f"{m['rise_per_12']:.1f}/12 | faces {round(m['azimuth_deg']) % 360} deg | "
f"area {m['area_m2']:.1f} m2")
Running it prints:
Plane 1: tilt 26.5 deg | 6.0/12 | faces 0 deg | area 67.2 m2
Plane 2: tilt 26.6 deg | 6.0/12 | faces 180 deg | area 66.4 m2
That matches the geometry. Each half should be 6/12, one facing north and one facing south. The true surface area of each is 12 x 5 / cos(26.57 degrees), about 67.1 square meters, and the small differences are noise.
Using real survey data
Drone processing software usually exports point clouds as LAS or LAZ files. Load one with laspy:
import laspy
las = laspy.read("roof.laz")
points = np.column_stack([las.x, las.y, las.z])
Then crop it to a single building before segmenting. Running RANSAC on a whole site, with ground and trees, will find the ground first.
Gotchas to watch for
Use a projected coordinate system. If your data is in latitude and longitude, x and y are in degrees while z is in meters, and every number will be wrong. Reproject to a local system such as UTM first.
Remove the ground and vegetation. Classify and filter before segmenting, or crop tightly to the roof. A common step is to keep only points above a minimum height.
Convex hulls overestimate non-convex planes. An L-shaped roof section will report more area than it has. For production work, replace the hull with an alpha shape or rasterize the inlier points and count the cells.
Dormers, vents and chimneys create small planes. The min_points setting filters most of them, but you may want to handle obstructions explicitly, since solar designers need to know where they are.
Check the pitch against a known reference. If you can measure one roof by hand, or you have the plans, compare it. A consistent offset usually points to a scale or coordinate problem.
Use a library for large clouds. Plain NumPy is fine for tens of thousands of points. For millions, Open3D has a built-in segment_plane method that is much faster.
Where this fits in a real workflow
Plane fitting is one step in a longer chain: plan the flight, capture overlapping images, process them into a point cloud, segment the roof, measure each plane, and produce a report or design file. Teams that do this often can either assemble the chain themselves, as in this tutorial, or use a tool that does it for them. Full disclosure: I work on Sightflight, drone software for solar and roofing, which is built around that second option. Either way, understanding what happens under the hood makes it easier to judge any tool's output.
Wrapping up
With SVD, RANSAC and a little trigonometry you can turn a pile of points into the numbers a roofer or solar designer actually uses. Try it on one of your own scans, compare the output with a hand measurement, and tune the threshold until the planes match what you see on the roof.
If you extend this, for example with alpha-shape areas or obstruction detection, share what you built in the comments.
