Skip to content

Shadow casting

More about shadow casting

In this notebook, we will compute which part of our building is casting shadow for the surroundings. Because of our high resolution voxelsize of 3.6m we are using interpolation with a low resolution voxelsize of 10.8m. The purpose of these computations is to remove the voxels that are causing to much shadow.

The inputs of this notebook are the context, the low resolution voxelized envelope and the high resolution voxelized envelope.

0. Initialization

import os
import topogenesis as tg
import pyvista as pv
import trimesh as tm
import numpy as np
from ladybug.sunpath import Sunpath
import pandas as pd

1. Import Meshes (context + envelope)

1.1. Load Meshes

context_path = os.path.relpath('../data/context_with_park2.obj')

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

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

2. Import Lattice

2.1. Load the Envelope Lattice in low resolution

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

3. Sun Vectors

3.1. Compute Sun Vectors

# initiate sunpath
sp = Sunpath(longitude=4.3571, latitude=52.0116)

hoys = []
sun_vectors = []
day_multiples = 30
for d in range(365):
    if d%day_multiples==0:
        for h in range(24):
            i = d*24 + h
            # compute the sun object
            sun = sp.calculate_sun_from_hoy(i)
            # extract the sun vector
            sun_vector = sun.sun_vector.to_array()
            # apparantly, if the Z component of sun vector is positive, 
            # it is under the horizon 
            if sun_vector[2] < 0.0:
                hoys.append(i)
                sun_vectors.append(sun_vector)

sun_vectors = np.array(sun_vectors)
# compute the rotation matrix 
Rz = tm.transformations.rotation_matrix(np.radians(36.324), [0,0,1])
# Rotate the sun vectors to match the site rotation
sun_vectors = tm.transform_points(sun_vectors, Rz)
print(sun_vectors.shape)

full_lattice = envelope_lattice * 0 + 1

# 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
full_lattice.fast_vis(p)

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

# add the sun locations, color orange
p.add_points( - sun_vectors * 300, color='#ffa500')

# 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 sun direction (other direction) from the sun vectors in a numpy array
sun_dirs = np.array(sun_vectors)
# exract the centroids of the envelope voxels

vox_cens = full_lattice.centroids

ray_dir = []
ray_src = []
for v_cen in vox_cens:
    for s_dir in sun_dirs:
        ray_dir.append(s_dir)
        ray_src.append(v_cen)
# 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 :",sun_dirs.shape)
print("number of rays to be shooted :",ray_src.shape)

4.2. Computing the Intersection

# computing the intersections 
tri_id, ray_id = context_mesh.ray.intersects_id(ray_origins=ray_src, ray_directions=ray_dir, multiple_hits=False)

5. Aggregate Simulation Result in the Shadow Casting Lattice

5.1. Compute the percentage of time that each ray intersects

# initializing the hits list full of zeros
hits = [0]*len(ray_dir)
# setting the rays that had an intersection to 1
for id in ray_id:
    hits[id] = 1

sun_count = len(sun_dirs)
vox_count = len(vox_cens)
# initiating the list of ratio
vox_shadow_casting = []
# 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(sun_count):
        # computing the ray id from voxel id and sun id
        r_id = v_id * sun_count + s_id

        # summing the intersections
        int_count += hits[r_id]

    # computing the percentage of the rays that DID NOT have an intersection
    shadow_casting = int_count/sun_count

    # add the ratio to list
    vox_shadow_casting.append(shadow_casting)


hits = np.array(hits)
vox_shadow_casting = np.array(vox_shadow_casting)

5.2. Store shadow casting low resolution 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
all_vox_shadow_casting = []

# 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 the voxel was outside...
    if vox_in == True:
        # read its value of shadow casting and append it to the list of all voxel shadow casting
        all_vox_shadow_casting.append(vox_shadow_casting[v_id])
        # add one to the voxel id so the next time we read the next voxel
        v_id += 1
    # if the voxel was not inside... 
    else:
        # add 0.0 for its shadow casting
        all_vox_shadow_casting.append(0.0)

# convert to array
shadowcasting_array = np.array(all_vox_shadow_casting)

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

# convert to lattice
shadowcasting_lattice = tg.to_lattice(shadowcasting_array, envelope_lattice)

5.3. Visualize the shadow casting lattice in low resolution

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

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

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

# Add the data values to the cell data
grid.point_arrays["Shadow Casting"] = shadowcasting_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.5, 1.0],opacity=opacity, shade=True)

# plotting|
p.show(use_ipyvtk=True)

6. Save the Shadow casting Lattice in low resolution into a CSV

# save the shadow casting latice to csv
csv_path = os.path.relpath('../data/shadow_casting_lowres.csv')
shadowcasting_lattice.to_csv(csv_path)

# add a small number to make sure that it is showing all the voxels of the envelope
shadowcasting_lattice += 0.01

7. Interpolation

7.1. Import the high resolution lattice

# loading the lattice from csv
lattice_path = os.path.relpath('../data/voxelized_envelope_green_3.6_final.csv')
highres_env_lattice = tg.lattice_from_csv(lattice_path)

7.2. Interpolate the high res lattice with the shadow casting lattice csv

from scipy.interpolate import RegularGridInterpolator

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

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

    # high_res lattice
    full_lattice = env_lattice + 1

    # sample points
    sample_points = full_lattice.centroids

    # interpolation
    interpolated_values = interpolating_function(sample_points)

    # lattice construction
    shadowcasting_lattice = tg.to_lattice(interpolated_values.reshape(env_lattice.shape), env_lattice)

    # nulling the unavailable cells
    shadowcasting_lattice *= env_lattice

    return shadowcasting_lattice

highres_shadowcasting_lattice = interpolate(shadowcasting_lattice,highres_env_lattice)

7.3. Visualize the shadow casting lattice in high resolution

# 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

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

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

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

# Add the data values to the cell data
grid.point_arrays["Shadow Casting"] = base_lattice.flatten(order="F")  # Flatten the Lattice

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

# plotting
p.show(use_ipyvtk=True)

8. Save the Shadow casting Lattice in high resolution into a CSV

# save the shadow casting latice to csv
csv_path = os.path.relpath('../data/shadow_casting_highres.csv')
highres_shadowcasting_lattice.to_csv(csv_path)

9. Envelope selection

# extra import function
def lattice_from_csv(file_path):
    # read metadata
    meta_df = pd.read_csv(file_path, nrows=3)

    shape = np.array(meta_df['shape'])
    unit = np.array(meta_df['unit'])
    minbound = np.array(meta_df['minbound'])

    # read lattice
    lattice_df = pd.read_csv(file_path, skiprows=5)

    # create the buffer
    buffer = np.array(lattice_df['value']).reshape(shape)

    # create the lattice
    l = tg.to_lattice(buffer, minbound=minbound, unit=unit)

    return l

# loading the lattice from csv
shadow_casting_path = os.path.relpath('../data/shadow_casting_highres.csv')
shadow_casting_lattice = lattice_from_csv(shadow_casting_path)

9.1. Visualizing the selection

p = pv.Plotter(notebook=True)

base_lattice = shadow_casting_lattice

# Set the grid dimensions: shape + 1 because we want to inject our values on the CELL data
grid = pv.UniformGrid()
grid.dimensions = np.array(base_lattice.shape) + 1
# The bottom left corner of the data set
grid.origin = base_lattice.minbound - base_lattice.unit * 0.5
# These are the cell sizes along each axis
grid.spacing = base_lattice.unit 

# adding the meshes
p.add_mesh(tri_to_pv(context_mesh), style='surface')

def create_mesh(value):

    lattice = np.copy(shadow_casting_lattice)
    lattice[shadow_casting_lattice > value] *= 0.0
    # Add the data values to the cell data
    grid.cell_arrays["Agents"] = lattice.flatten(order="F")  # Flatten the array!
    # filtering the voxels
    threshed = grid.threshold([0.001, 1.0])
    # adding the voxels
    p.add_mesh(threshed, name='sphere', show_edges=True, opacity=1.0, show_scalar_bar=False, clim=[0.0, 1.0])

    return

p.add_slider_widget(create_mesh, [0, 1], title='', value=1.0, event_type="always", style="classic", pointa=(0.1, 0.1), pointb=(0.9, 0.1))
p.show(use_ipyvtk=True)

shadow_casting_lattice.max()

9.2. Generating an envelope based on the selection

lower_bound = 0.001
upper_bound = 0.471
lower_condition = shadow_casting_lattice > lower_bound
upper_condition = shadow_casting_lattice < upper_bound
new_avail_lattice = lower_condition * upper_condition

9.3. Visualize the new available lattice high resolution

p = pv.Plotter(notebook=True)

# adding the avilability lattice
new_avail_lattice.fast_vis(p)

p.show(use_ipyvtk=True)

9.4. Save new high resolution envelope to CSV

csv_path = os.path.relpath('../data/envelope_shadow_3.6.csv')
new_avail_lattice.to_csv(csv_path)

9.5. Low resolution selection

# loading the lattice from csv
shadow_casting_path = os.path.relpath('../data/shadow_casting_lowres.csv')
shadow_casting_lattice = lattice_from_csv(shadow_casting_path)

lower_bound = 0.001
upper_bound = 0.471
lower_condition = shadow_casting_lattice > lower_bound
upper_condition = shadow_casting_lattice < upper_bound
new_avail_lattice = lower_condition * upper_condition

9.6. Visualize the new available lattice low resolution

p = pv.Plotter(notebook=True)

# adding the avilability lattice
new_avail_lattice.fast_vis(p)

p.show(use_ipyvtk=True)

9.7. Save the new low resolution envelope to CSV

csv_path = os.path.relpath('../data/envelope_shadow_10.8.csv')
new_avail_lattice.to_csv(csv_path)

Credits

__author__ = "Shervin Azadi and Pirouz Nourian"
__changes_made_by__ = "Lotte Zwolsman"
__license__ = "MIT"
__version__ = "1.0"
__url__ = "https://github.com/frankvahstal/spatial_computing_workshops"
__summary__ = "Spatial Computing Design Studio Workshop on Solar Envelope"