Skylight

This script calculates the skylight acces and the skylight blocking on the context of an envelope.

It creates a sphere around the envelope and then castst a ray from each voxel to each of the points on the sphere. This ray could be blocked by the context. If the ray is not blocked by the context, it then also casts an inverse ray from the voxel towards. If this intersects the context, then the particular voxel blocks skylight from the context.

As an iput this script takes a envelope_highres and envelope_lowres. As an output it gives a sky_shadowcasting_highres and sky_access highres lattice.

0. Initialization

Importing all necessary libraries and specifying the inputs

import os
import topogenesis as tg
import pyvista as pv
import trimesh as tm
import numpy as np
from ladybug.sunpath import Sunpath
from scipy.interpolate import RegularGridInterpolator

1. Import Meshes (context + envelope)

1.1. Load Meshes

envelope_path = os.path.relpath("../data/compulsory_envelope.obj")
context_path = os.path.relpath("../data/immediate_context.obj")

# load the mesh from file
envelope_mesh = tm.load(envelope_path)
context_mesh = tm.load(context_path)

# Check if the mesh is watertight
print(envelope_mesh.is_watertight)

1.2. Visualize Meshes (with pyvista)

# convert trimesh to pv_mesh
def tri_to_pv(tri_mesh):
    faces = np.pad(tri_mesh.faces, ((0, 0),(1,0)), 'constant', constant_values=3)
    pv_mesh = pv.PolyData(tri_mesh.vertices, faces)
    return pv_mesh

# initiating the plotter
p = pv.Plotter(notebook=True)

# adding the meshes
p.add_mesh(tri_to_pv(envelope_mesh), color='#abd8ff')
p.add_mesh(tri_to_pv(context_mesh), color='#aaaaaa')

# plotting
p.show(use_ipyvtk=True)

2. Import Lattice (envelope)

2.1. Load the Envelope Lattice

# loading the lattice from csv
lattice_path = os.path.relpath('../data/voxelized_envelope_lowres.csv')
envelope_lattice = tg.lattice_from_csv(lattice_path)

full_lattice = envelope_lattice * 0 + 1

# loading the lattice from csv
lattice_path = os.path.relpath('../data/voxelized_envelope_highres.csv')
avail_lattice_highres = tg.lattice_from_csv(lattice_path)

2.2. Visualize the Context Mesh + Envelope Lattice

# convert mesh to pv_mesh
def tri_to_pv(tri_mesh):
    faces = np.pad(tri_mesh.faces, ((0, 0),(1,0)), 'constant', constant_values=3)
    pv_mesh = pv.PolyData(tri_mesh.vertices, faces)
    return pv_mesh

# initiating the plotter
p = pv.Plotter(notebook=True)

# fast visualization of the lattice
envelope_lattice.fast_vis(p)

# adding the meshes
p.add_mesh(tri_to_pv(context_mesh), color='#aaaaaa')

# plotting
p.show(use_ipyvtk=True)

3. Skylight factors

3.1. Compute Skylight Vectors

#Create a sphere to put points on that represent the sky 
sphere_mesh = tm.creation.icosphere(subdivisions=3, radius= 300.0, color=None)
sphere_vectors = np.copy(sphere_mesh.vertices)
sky_vectors = []
for v in sphere_vectors:
    if v[2] > 0.0:
        sky_vectors.append(v)


sky_vectors = np.array(sky_vectors)

# convert mesh to pv_mesh
def tri_to_pv(tri_mesh):
    faces = np.pad(tri_mesh.faces, ((0, 0),(1,0)), 'constant', constant_values=3)
    pv_mesh = pv.PolyData(tri_mesh.vertices, faces)
    return pv_mesh


# initiating the plotter
p = pv.Plotter(notebook=True)

# fast visualization of the lattice
envelope_lattice.fast_vis(p)

# adding the meshes
p.add_mesh(tri_to_pv(sphere_mesh), color='#aaaaaa')

# add the sun locations, color orange
p.add_points( sky_vectors * 300, color='#0033ff')

# plotting
p.show(use_ipyvtk=True)

4. Compute Intersection of Sun Rays with Context Mesh

4.1. Preparing the List of Ray Directions and Origins

# constructing the sky direction from the sky vectors in a numpy array
sky_dirs = np.array(sky_vectors)

# exract the centroids of the envelope voxels
full_lattice = envelope_lattice * 0 + 1
vox_cens = full_lattice.centroids

# next step we need to shoot in all of the sky directions from all of the voxels 
ray_dir = []
ray_src = []
for v_cen in vox_cens:
    for s_dir in sky_dirs:
        #if s_dir = envelope_mesh
        ray_src.append(v_cen)
        ray_dir.append(s_dir)

# converting the list of directions and sources to numpy array
ray_dir = np.array(ray_dir)
ray_src = np.array(ray_src)

print("number of voxels to shoot rays from :",vox_cens.shape)
print("number of rays per each voxel :",sky_dirs.shape)
print("number of rays to be shot :",ray_src.shape)

4.2. Computing the Intersection

# computing the intersections of rays with the context from voxel to sky 
f_tri_id, f_ray_id = context_mesh.ray.intersects_id(ray_origins=ray_src, ray_directions=ray_dir, multiple_hits=False)

# computing the intersections of rays with the context backwards from voxel to context 
b_tri_id, b_ray_id = context_mesh.ray.intersects_id(ray_origins=ray_src, ray_directions=-ray_dir, multiple_hits=False)

5. Aggregate Simulation Result in the Sun Access Lattice

5.1. Compute the percentage of time that each voxel sees the sun

# initializing the hits list full of zeros
f_hits = [0]*len(ray_dir)
b_hits = [0]*len(ray_dir)

for id in f_ray_id:
    f_hits[id] = 1

for id in b_ray_id:
    if not f_hits[id]:
        b_hits[id] = 1

sky_count = len(sky_dirs)
vox_count = len(vox_cens)
# initiating the list of ratio
vox_sky_acc = []
# iterate over the voxels
for v_id in range(vox_count):
    # counter for the intersection
    int_count = 0
    # iterate over the sun rays
    for s_id in range(sky_count):
        # computing the ray id from voxel id and sun id
        r_id = s_id + v_id * sky_count 

        # summing the intersections
        int_count += f_hits[r_id]

    # computing the percentage of the rays that DID NOT have 
    # an intersection (aka could see the sun)
    sky_access =   1 - int_count/sky_count

    # add the ratio to list
    vox_sky_acc.append(sky_access)

# hits = np.array(hits)
vox_sky_acc = np.array(vox_sky_acc)

5.1. Compute the shadow casting

# initiating the list of ratio
vox_sky_shadowing = []
# iterate over the voxels
for v_id in range(vox_count):
    # counter for the intersection
    int_count = 0
    # iterate over the sun rays
    for s_id in range(sky_count):
        # computing the ray id from voxel id and sun id
        r_id = s_id + v_id * sky_count 

        # summing the intersections
        int_count += b_hits[r_id]

    # computing the percentage of the rays that DID NOT have an intersection
    shadowing = int_count  /  sky_count
    # add the ratio to list
    vox_sky_shadowing.append(shadowing)

vox_sky_shadowing = np.array(vox_sky_shadowing)

5.2. Store sun access information in a Lattice

# getting the condition of all voxels: are they inside the envelop or not
env_all_vox = full_lattice.flatten()

# all voxels sun access
all_vox_sky_acc = []

# v_id: voxel id in the list of only interior voxels
v_id = 0

# for all the voxels, place the interiority condition of each voxel in "vox_in"
for vox_in in env_all_vox:
    if vox_in == True:
        # read its value of sun access and append it to the list of all voxel sun access
        all_vox_sky_acc.append(vox_sky_acc[v_id])
        # add one to the voxel id so the next time we read the next voxel
        v_id += 1 
    else:
        # add 0.0 for its sun access
        all_vox_sky_acc.append(0.0)

# convert to array
sky_acc_array = np.array(all_vox_sky_acc)

# reshape to lattice shape
sky_acc_array = sky_acc_array.reshape(envelope_lattice.shape)

# convert to lattice
sky_acc_lattice = tg.to_lattice(sky_acc_array, envelope_lattice)

5.2. Store shadow casting information in a Lattice

# getting the condition of all voxels: are they inside the envelop or not
env_all_vox = full_lattice.flatten()

# all voxels sun access
all_vox_sky_shad_cast = []

# v_id: voxel id in the list of only interior voxels
v_id = 0

# for all the voxels, place the interiority condition of each voxel in "vox_in"
for vox_in in env_all_vox:
    if vox_in == True:
        # read its value of sun access and append it to the list of all voxel sun access
        all_vox_sky_shad_cast.append(vox_sky_shadowing[v_id])
        # add one to the voxel id so the next time we read the next voxel
        v_id += 1 
    else:
        # add 0.0 for its sun access
        all_vox_sky_shad_cast.append(0.0)


# convert to array
sky_shadcast_array = np.array(all_vox_sky_shad_cast)

# reshape to lattice shape
sky_shadcast_array = sky_shadcast_array.reshape(envelope_lattice.shape)

# convert to lattice
sky_shadcast_lattice = tg.to_lattice(sky_shadcast_array, envelope_lattice)

5.3. Visualize the sun access lattice

# initiating the plotter
p = pv.Plotter(notebook=True)

vis_lattice = sky_acc_lattice
# Create the spatial reference
grid = pv.UniformGrid()

# Set the grid dimensions: shape because we want to inject our values
grid.dimensions = vis_lattice.shape
# The bottom left corner of the data set
grid.origin = vis_lattice.minbound
# These are the cell sizes along each axis
grid.spacing = vis_lattice.unit

# Add the data values to the cell data
grid.point_arrays["Sky Access"] = vis_lattice.flatten(order="F")  # Flatten the Lattice

# adding the meshes
p.add_mesh(tri_to_pv(context_mesh), opacity=0.1, style='wireframe')

# adding the volume
opacity = np.array([0,0.6,0.6,0.6,0.6,0.6,0.6])*1.5
p.add_volume(grid, cmap="coolwarm", clim=[0.0, 1.0],opacity=opacity, shade=True)

# plotting
p.show(use_ipyvtk=True)

Visualize the shadow casting lattice

# initiating the plotter
p = pv.Plotter(notebook=True)

vis_lattice = sky_shadcast_lattice
# Create the spatial reference
grid = pv.UniformGrid()

# Set the grid dimensions: shape because we want to inject our values
grid.dimensions = vis_lattice.shape
# The bottom left corner of the data set
grid.origin = vis_lattice.minbound
# These are the cell sizes along each axis
grid.spacing = vis_lattice.unit

# Add the data values to the cell data
grid.point_arrays["Skylight Blocking"] = vis_lattice.flatten(order="F")  # Flatten the Lattice

# adding the meshes
p.add_mesh(tri_to_pv(context_mesh), opacity=0.1, style='wireframe')

# adding the volume
opacity = np.array([0,0.6,0.6,0.6,0.6,0.6,0.6])*1.5
p.add_volume(grid, cmap="coolwarm", clim=[0., 1.0],opacity=opacity, shade=True)

# plotting
p.show(use_ipyvtk=True)

5.4 interpolation

interpolation sky factor

def interpolate(info_lowres, base_highres):
    # line spaces
    x_space = np.linspace(info_lowres.minbound[0], info_lowres.maxbound[0],info_lowres.shape[0])
    y_space = np.linspace(info_lowres.minbound[1], info_lowres.maxbound[1],info_lowres.shape[1])
    z_space = np.linspace(info_lowres.minbound[2], info_lowres.maxbound[2],info_lowres.shape[2])

    # interpolation function
    interpolating_function = RegularGridInterpolator((x_space, y_space, z_space), info_lowres, bounds_error=False, fill_value=None)

    # high_res lattice
    envelope_lattice = base_highres + 1

    # sample points
    sample_points = envelope_lattice.centroids

    # interpolation
    interpolated_values = interpolating_function(sample_points)

    # lattice construction
    info_highres = tg.to_lattice(interpolated_values.reshape(base_highres.shape), base_highres)

    # nulling the unavailable cells
    info_highres *= base_highres

    return info_highres

#interpolate the lattice over highres
sky_acc_highres = interpolate(sky_acc_lattice, avail_lattice_highres)

interpolation sky blocking factor

def interpolate(info_lowres, base_highres):
    # line spaces
    x_space = np.linspace(info_lowres.minbound[0], info_lowres.maxbound[0],info_lowres.shape[0])
    y_space = np.linspace(info_lowres.minbound[1], info_lowres.maxbound[1],info_lowres.shape[1])
    z_space = np.linspace(info_lowres.minbound[2], info_lowres.maxbound[2],info_lowres.shape[2])

    # interpolation function
    interpolating_function = RegularGridInterpolator((x_space, y_space, z_space), info_lowres, bounds_error=False, fill_value=None)

    # high_res lattice
    envelope_lattice = base_highres + 1

    # sample points
    sample_points = envelope_lattice.centroids

    # interpolation
    interpolated_values = interpolating_function(sample_points)

    # lattice construction
    info_highres = tg.to_lattice(interpolated_values.reshape(base_highres.shape), base_highres)

    # nulling the unavailable cells
    info_highres *= base_highres

    return info_highres

#interpolate the lattice over highres
sky_shadcast_highres = interpolate(sky_shadcast_lattice, avail_lattice_highres)

5.5 visualize interpolation

Sky access interpolation

# initiating the plotter
p = pv.Plotter(notebook=True)

vis_lattice = sky_acc_highres
# Create the spatial reference
grid = pv.UniformGrid()

# Set the grid dimensions: shape because we want to inject our values
grid.dimensions = vis_lattice.shape
# The bottom left corner of the data set
grid.origin = vis_lattice.minbound
# These are the cell sizes along each axis
grid.spacing = vis_lattice.unit

# Add the data values to the cell data
grid.point_arrays["Sky access Highres"] = vis_lattice.flatten(order="F")  # Flatten the Lattice

# adding the meshes
p.add_mesh(tri_to_pv(context_mesh), opacity=0.1, style='wireframe')

# adding the volume
opacity = np.array([0,0.6,0.6,0.6,0.6,0.6,0.6])*1.5
p.add_volume(grid, cmap="coolwarm", clim=[0., 1.0],opacity=opacity, shade=True)

# plotting
p.show(use_ipyvtk=True)

sky view blocking interpolation

# initiating the plotter
p = pv.Plotter(notebook=True)

vis_lattice = sky_shadcast_highres
# Create the spatial reference
grid = pv.UniformGrid()

# Set the grid dimensions: shape because we want to inject our values
grid.dimensions = vis_lattice.shape
# The bottom left corner of the data set
grid.origin = vis_lattice.minbound
# These are the cell sizes along each axis
grid.spacing = vis_lattice.unit

# Add the data values to the cell data
grid.point_arrays["SkyBlocking highres"] = vis_lattice.flatten(order="F")  # Flatten the Lattice

# adding the meshes
p.add_mesh(tri_to_pv(context_mesh), opacity=0.1, style='wireframe')

# adding the volume
opacity = np.array([0,0.6,0.6,0.6,0.6,0.6,0.6])*1.5
p.add_volume(grid, cmap="coolwarm", clim=[0., 1.0],opacity=opacity, shade=True)

# plotting
p.show(use_ipyvtk=True)

6. Save Sun Access Lattice into a CSV

#save the sky access latice to csv
csv_path = os.path.relpath('../data/sky_access_highres.csv')
sky_acc_highres.to_csv(csv_path)

#save the skyblocking latice to csv
csv_path = os.path.relpath('../data/sky_shadowcasting_highres.csv')
sky_shadcast_highres.to_csv(csv_path)

Credits

__author__ = "Shervin Azadi and Pirouz Nourian"
__editor__ = "Siebren Meines"
__version__ = "1.0"
__url__ = "https://github.com/shervinazadi/spatial_computing_workshops"
__summary__ = "Skylight and skylight blocking calculation"