# OSM skillion roof with bmesh… need an alternative to ray_cast

OSM buildings with the skillion roof shape, are often tagged with only direction and levels. {'direction': '357', 'shape': 'skillion', 'levels': '1'}
{'direction': '267', 'shape': 'skillion', 'levels': '1'}
{'direction': '87', 'shape': 'skillion', 'levels': '1'}


The direction Vector can be obtained using

tag_direction = 270
v = Vector((sin(a), cos(a), 0)) # west (-1, 0, 0)


For example the default Cube could be considered a building, and given a 1 bu height, direction 0, skillion roof would look like Here is some code to produce skillion roof from default cube, alter height and direction d to test.

The code adds a new cube, uses its top face as the roof, from the centre of the face, ray casts in the slope direction for low point and in negative slope direction for the high point.

These intersection points are then used to project the roof verts to the pitched plane which is defined by low point and the resulting normal of the pitched plane.

import bpy
import bmesh
from math import sin, cos, radians, degrees
from mathutils import Vector

from bpy import context
cube = context.object
mesh = cube.data

d = 90
b = radians(d + 90) # ortho to direction
height = 1

up = Vector((0, 0, 1))
# direction vector
v = Vector((sin(a), cos(a), 0))
# ortho to direction vector
vo = Vector((sin(b), cos(b), 0))
bm = bmesh.new()
bm.from_mesh(mesh)
# roof face is index 5
bm.faces.ensure_lookup_table()
roof = bm.faces
centre = roof.calc_center_bounds()
# use ray cast to find low and high point.
low_pt = cube.ray_cast(centre, v)
high_pt = cube.ray_cast(centre, -v)
print(low_pt, high_pt)
# add height to high point
high_pt.z += height
#Vector Across roof
V = low_pt - high_pt
print("Roof Pitch", degrees(v.angle(V)))

# define a plane at low point with normal n
# n is normal vector of pitched roof
n = vo.cross(V).normalized()
p = low_pt # low point
# project verts in up direction to plane.
for v in roof.verts:
v.co = v.co - (v.co - p).dot(n) * up
# check out result
bm.to_mesh(mesh)
mesh.update()


The rub is that I don't believe I can use ray_cast when building the objects as a bmesh, without having to to_mesh, which for potentially thousands of buildings could be expensive.

What alternatives would I have to the ray_cast call that essentially find the intersection of a line with a polygon.

Finding the intersection between a plane and a line could be done using the mathutils module. There are a few nice functions in there, the most interesting one for you is mathutils.geometry.intersect_line_plane(). It basically takes four arguments, all of them Vectors, describing the setup as follows:

mathutils.geometry.intersect_line_plane(VecLineStart, VecLineEnd, VecPlaneOrigin, VecPlaneNormal)


it returns either the intersection point or None (if line and Plane are parallel).

Example:

x = mathutils.geometry.intersect_line_plane((0,0,0), (0,1,0), (0.33, 1.12, 2.44), (0,1,0))


gives you

Vector((0.0, 1.1200000047683716, 0.0))


which is the Intersection Point between the line from (0,0,0) to (0,1,0) - extended to infinity - and the plane defined by the origin at (0.33, 1.12, 2.44) with its normal pointing at (0,1,0).

• Cheers, Have done similar using mathutils.geometry.intersect_line_line_2d(...) which can use the (lat, lon) coords from the OSM data. Then I can step around these edges and find both the d and -d direction vectors intersection point. The hassle with intersect_line_plane (as I see it is..) I need the intersection points of the poly, to make the plane (the roof truss triangle), which is still a guestimation. I'm thinking a bounding box method might be better to make the slope plane, as the outer ring poly can be concave in parts. – batFINGER Oct 7 '16 at 7:03