Solar Envelope, Solar accessibility and Sky view factor computations

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
import pickle
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/optional_envelope.obj")
context_path = os.path.relpath("../data/immediate_context2.obj")
skydome_path = os.path.relpath("../data/skydome.obj")

# load the mesh from file
envelope_mesh = tm.load(envelope_path)
context_mesh = tm.load(context_path)
skydome_mesh = tm.load(skydome_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')

#trying the skydome
# p.add_mesh(tri_to_pv(skydome_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/lattice1620.csv")
envelope_lattice = tg.lattice_from_csv(lattice_path)

# creating the full lattice
full_lattice = envelope_lattice * 0 + 1

2.2. Visualize the Context Mesh + Envelope Lattice

pv.set_plot_theme("document")

# 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

# Visualize the mesh using pyvista plotter

# 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), opacity=0.1, style='wireframe')

# plotting
cpos = [(785.8704805788776, 708.4540755788776, 741.8613927288776),
 (65.08283250000001, -12.333572500000002, 21.07374465),
 (0.0, 0.0, 1.0)]
p.camera_position = cpos
p.window_size = 2000, 2000
p.show(use_ipyvtk=True)
# p.screenshot("solar_1")
print(p.camera_position)

3. Sun Vectors

3.1. Compute Sun Vectors

# initiate sunpath
sp = Sunpath(longitude=4.4777, latitude=51.9244)

# define sun hours : A list of hours of the year for each sun vector
# there are 8760 hours in a year, so the following integers refer to specific hours throughout the year
hoys = []
sun_vectors = []
day_multiple = 30

# for every day in 365 days of a year
for d in range(365):
    # filter days
    if d % day_multiple == 0:
        # for every hour in 24 hours of a day
        for h in range(24):
            # reconstruct the hoy
            hoy = d*24 + h
            # compute the sun object
            sun = sp.calculate_sun_from_hoy(hoy)
            # extract the sun vector
            sun_vector = sun.sun_vector.to_array()
            # remove sun vectors beneath the horizon
            if sun_vector[2] < 0.0:
                # add the hoy to the list
                hoys.append(hoy)
                # add the sun vetor to the list
                sun_vectors.append(sun_vector)

# convert sun_vectors list to numpy array                
sun_vectors = np.array(sun_vectors)
# compute the rotation matrix : (angle, axis)
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)

#checking how many rays to calculate
print(len(sun_vectors))

3.2. Computing SVF vectors with skydome, (exporting normals to list)

#adding normals of the dome to a numpy list
svf_vectors = np.array(skydome_mesh.face_normals)

#checking how many rays to calculate
print(len(svf_vectors))

3.3. Visualization

# Solar envelope and solar accessibility
# 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

# Visualize the mesh using pyvista plotter

# 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), opacity=0.1, style='wireframe')

# add the sun locations, color orange
for loc in (- sun_vectors * 250):
    p.add_mesh(pv.Sphere(radius= 2, center = loc), color="#ffa500", show_edges=False, lighting = False)

# plotting
cpos = [(785.8704805788776, 708.4540755788776, 741.8613927288776),
 (65.08283250000001, -12.333572500000002, 21.07374465),
 (0.0, 0.0, 1.0)]
p.camera_position = cpos
p.window_size = 2000, 2000
p.show(use_ipyvtk=True)
# p.screenshot("solar_2")
print(p.camera_position)

# SVF
# 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

# Visualize the mesh using pyvista plotter

# 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), opacity=0.1, style='wireframe')

# add the sun locations, color orange
for loc in (svf_vectors * 250):
    p.add_mesh(pv.Sphere(radius= 2, center = loc), color="#ffa500", show_edges=False, lighting = False)

# plotting
cpos = [(785.8704805788776, 708.4540755788776, 741.8613927288776),
 (65.08283250000001, -12.333572500000002, 21.07374465),
 (0.0, 0.0, 1.0)]
p.camera_position = cpos
p.window_size = 2000, 2000
p.show(use_ipyvtk=True)
# p.screenshot("solar_3")
print(p.camera_position)

4. Compute Intersection of Sun Rays with Context Mesh

4.1. Preparing the List of Ray Directions and Origins

# # Coreating list of directions and rays for solar accesibility and envelope

# constructing the sun direction from the sun vectors in a numpy array
converted_sun_vectors = sun_vectors.astype('float16')

sun_dirs = -np.array(converted_sun_vectors)

# exract the centroids of the envelope voxels
vox_cens = full_lattice.centroids

# next step we need to shoot in all of the sun directions from all of the voxels, todo so, 
# we need repeat the sun direction for the number of voxels to construct the ray_dir 
# (which is the list of all ray directions).

# creating a second ray_dir and ray_src list, also keeping track
# of which voxels need to shoot a reversed ray with an index

ray_dir = []
ray_src = []
ray_dir2 = []
ray_src2 = []
ray_dir2_index = []

for v_cen in vox_cens:
        for s_dir in sun_dirs:
            ray_src.append(v_cen)
            ray_dir.append(s_dir)

# checking if the ray intersects with any geometry, if not, appending the reversed ray to a second src snd dir list
# Also creating an index list so that the voxel id can be correctly reinserted later on

for v_cen in vox_cens:
        for s_dir in sun_dirs:
            if context_mesh.ray.intersects_any(ray_origins=[v_cen], ray_directions=[s_dir]):
                ray_dir2_index.append(0)
                continue
            else:
                ray_src2.append(v_cen)
                ray_dir2.append(-s_dir) # adding the reversed ray
                ray_dir2_index.append(1)


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

# Creating list of directions and rays for the SVF

converted_svf_vectors = svf_vectors.astype('float16')
sun_dirs2 = np.array(converted_svf_vectors)
vox_cens = full_lattice.centroids

ray_dir3 = []
ray_src3 = []

for v_cen in vox_cens:
        for s_dir in sun_dirs2:
            ray_src3.append(v_cen)
            ray_dir3.append(s_dir)

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

# checking how many calculations to do for solar accesibility and solar envelope
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 shot :",ray_src.shape)

# checking how many calculations to do for SVF
print("number of voxels to shoot rays from :",vox_cens.shape)
print("number of rays per each voxel :",sun_dirs2.shape)
print("number of rays to be shot :",ray_src3.shape)

4.2. Computing the Intersection

# computing the intersections of rays with the context mesh for the Solar accesibility and solar envelope
tri_id, ray_id = context_mesh.ray.intersects_id(ray_origins=ray_src, ray_directions=ray_dir, multiple_hits=False)

# computing the intersections of inversed rays with the context mesh for the solar envelope
tri_id2, ray_id2 = context_mesh.ray.intersects_id(ray_origins=ray_src2, ray_directions=ray_dir2, multiple_hits=False)

print("This way of computing saved: " + str(len(ray_dir) - len(ray_dir2)) + " calculations for the inverted ray directions, saving about 3 or 4 miniutes of calculations :)")

# computing the intersections of rays with the context mesh for the SVF
tri_id3, ray_id3 = context_mesh.ray.intersects_id(ray_origins=ray_src3, ray_directions=ray_dir3, 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

hits = [0]*len(ray_dir) # sun_acc
hits2 = [0]*len(ray_dir2) # sun_block

hits3 = [] # sun_block with correct length

# setting the rays that had an intersection to 1
for id in ray_id:
    hits[id] = 1 # sun_acc

for id in ray_id2:
    hits2[id] = 1 # sun_block

# fixing the length of the hits list and setting 'ignored' non intersections to 0
hcounter = 0
for i in ray_dir2_index:
    if i == 1:
        hits3.append(hits2[hcounter])
        hcounter += 1
    else:
        hits3.append(0)

sun_count = len(sun_dirs)
vox_count = len(vox_cens)

# initiating the list of ratio
vox_sun_acc = []
vox_sun_blockage = []

# iterate over the voxels
for v_id in range(vox_count):
    # counter for the intersection
    int_count = 0 # sun_acc
    int_count2 = 0 # sun_block
    # iterate over the sun rays
    for s_id in range(sun_count):
        # computing the ray id from voxel id and sun id
        r_id = s_id + v_id * sun_count

        # summing the intersections
        int_count += hits[r_id]
        int_count2 += hits3[r_id]


    # computing the percentage of the rays that DID NOT have 
    # an intersection (aka could see the sun / did not block sunlight)
    sun_access = 1.0 - int_count / sun_count
    sun_blockage = 1.0 - int_count2 / sun_count

    # add the ratio to list
    vox_sun_acc.append(sun_access)
    vox_sun_blockage.append(sun_blockage)

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

# Same calculations, but for the SVF
# initializing the hits list full of zeros

hits4 = [0]*len(ray_dir3)

for id in ray_id3:
    hits4[id] = 1 

sun_count = len(sun_dirs2)
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(sun_count):
        # computing the ray id from voxel id and sun id
        r_id = s_id + v_id * sun_count

        # summing the intersections
        int_count += hits4[r_id]

    # computing the percentage of the rays that DID NOT have 
    # an intersection (aka could see the skydome)
    sky_access = 1.0 - int_count / sun_count

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

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

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 / blockage
all_vox_sun_acc = []
all_vox_sun_block = []

# 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 sun access and sun blockage and append it to the list of all voxel sun access / sun blockage
        all_vox_sun_acc.append(vox_sun_acc[v_id])
        all_vox_sun_block.append(vox_sun_blockage[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 sun access / sun blockage
        all_vox_sun_acc.append(0.0)
        all_vox_sun_block.append(0.0)

# convert to array
sunacc_array = np.array(all_vox_sun_acc)
sunblock_array = np.array(all_vox_sun_block)

# reshape to lattice shape
sunacc_array = sunacc_array.reshape(full_lattice.shape)
sunblock_array = sunblock_array.reshape(full_lattice.shape)

# convert to lattice
sunacc_lattice = tg.to_lattice(sunacc_array, full_lattice)
sunblock_lattice = tg.to_lattice(sunblock_array, full_lattice)

# Same calculation, but for SVF
# getting the condition of all voxels: are they inside the envelop or not
env_all_vox = full_lattice.flatten()

# all voxels sky 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 the voxel was outside...
    if vox_in == True:
        # read its value of sky acces and append it to the list of all voxel sky 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
    # if the voxel was not inside... 
    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)

# Further info: for vectorized version of this code check: https://github.com/shervinazadi/spatial_computing_workshops/blob/master/notebooks/w2_solar_envelope.ipynb

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

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

5.3.1. Visualize the sun access 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 = sunacc_lattice.shape
# The bottom left corner of the data set
grid.origin = sunacc_lattice.minbound
# These are the cell sizes along each axis
grid.spacing = sunacc_lattice.unit

# Add the data values to the cell data
grid.point_arrays["Sun acces"] = sunacc_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=False)

# plotting
cpos = [(785.8704805788776, 708.4540755788776, 741.8613927288776),
 (65.08283250000001, -12.333572500000002, 21.07374465),
 (0.0, 0.0, 1.0)]
p.camera_position = cpos
p.window_size = 2000, 2000
p.show(use_ipyvtk=True)
# p.screenshot("solar_4")
print(p.camera_position)

5.3.2. Visualize solar envelope

# 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 = sunblock_lattice.shape
# The bottom left corner of the data set
grid.origin = sunblock_lattice.minbound
# These are the cell sizes along each axis
grid.spacing = sunblock_lattice.unit

# Add the data values to the cell data
grid.point_arrays["Solar envelope (where 0 is blocking sunlight for neighbours)"] = sunblock_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=False)

# plotting
cpos = [(785.8704805788776, 708.4540755788776, 741.8613927288776),
 (65.08283250000001, -12.333572500000002, 21.07374465),
 (0.0, 0.0, 1.0)]
p.camera_position = cpos
p.window_size = 2000, 2000
p.show(use_ipyvtk=True)
# p.screenshot("solar_5")
print(p.camera_position)

5.3.3. Visualize the sky view factor

# 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 = sky_acc_lattice.shape
# The bottom left corner of the data set
grid.origin = sky_acc_lattice.minbound
# These are the cell sizes along each axis
grid.spacing = sky_acc_lattice.unit

# Add the data values to the cell data
grid.point_arrays["Sky view factor"] = sky_acc_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=False)

# plotting
cpos = [(785.8704805788776, 708.4540755788776, 741.8613927288776),
 (65.08283250000001, -12.333572500000002, 21.07374465),
 (0.0, 0.0, 1.0)]
p.camera_position = cpos
p.window_size = 2000, 2000
p.show(use_ipyvtk=True)
# p.screenshot("solar_6")
print(p.camera_position)

6. Interpolate lowres lattices to highres

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

#loading the interpolation function

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 point
    sample_points = full_lattice.centroids

    # interpolation
    interpolated_values = interpolating_function(sample_points)

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

    # nulling the unavailable cells
    sunacc_lattice *= env_lattice

    return sunacc_lattice

6.1. Interpolate sun acces lattice

highres_sunacc_lattice = interpolate(sunacc_lattice,highres_env_lattice)

6.2. Interpolate solar envelope

highres_sunblock_lattice = interpolate(sunblock_lattice,highres_env_lattice)

6.3. Interpolate SVF

highres_sky_acc_lattice = interpolate(sky_acc_lattice,highres_env_lattice)

7. Setting the solar envelope to be a boolean lattice

highres_solar_envelope_flat = highres_sunblock_lattice.flatten()
full_highres_lattice = highres_env_lattice * 0 + 1

for i in range(len(highres_solar_envelope_flat)):
    if highres_solar_envelope_flat[i] < 0.75: #vary this value to increase or decrease the envelope size
        highres_solar_envelope_flat[i] = False
    else: 
        highres_solar_envelope_flat[i] = True

highres_solar_envelope_bool = (highres_solar_envelope_flat > 0).tolist() 

# The soon to be exported csv file now no longer contains "0" and "1", but "False" and "True" as boolean values, 
# if this is not case, topogenesis on windows messes up  loading of the file, 
# so this needs to be done in order to open the envelope in a different jupyter notebook.

highres_solar_envelope_bool = np.array(highres_solar_envelope_bool) #array conversion

highres_solar_envelope_bool = highres_solar_envelope_bool.reshape(full_highres_lattice.shape) #lattice shape

highres_solar_envelope_minmax = tg.to_lattice(highres_solar_envelope_bool, full_highres_lattice) #lattice creation

#Visualizing the final envelope

# 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

# Visualize the mesh using pyvista plotter
#######

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

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

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

# plotting
cpos = [(785.8704805788776, 708.4540755788776, 741.8613927288776),
 (65.08283250000001, -12.333572500000002, 21.07374465),
 (0.0, 0.0, 1.0)]
p.camera_position = cpos
p.window_size = 2000, 2000
p.show(use_ipyvtk=True)
# p.screenshot("solar_7")
print(p.camera_position)

7.2.1. Visualize sun acces

# 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_sunacc_lattice * highres_solar_envelope_minmax
# 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["Sun Access"] = base_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])
p.add_volume(grid, cmap="coolwarm", clim=[0.5, 1.0],opacity=opacity, shade=False)

# plotting
cpos = [(785.8704805788776, 708.4540755788776, 741.8613927288776),
 (65.08283250000001, -12.333572500000002, 21.07374465),
 (0.0, 0.0, 1.0)]
p.camera_position = cpos
p.window_size = 2000, 2000
p.show(use_ipyvtk=True)
# p.screenshot("solar_8")
print(p.camera_position)

7.2.2. Visualize solar envelope

# 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_sunblock_lattice * highres_solar_envelope_minmax

# 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["Solar envelope (where 0 is blocking sunlight for neighbours)"] = base_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])
p.add_volume(grid, cmap="coolwarm", clim=[0.75, 1.0],opacity=opacity, shade=False)

# plotting
cpos = [(785.8704805788776, 708.4540755788776, 741.8613927288776),
 (65.08283250000001, -12.333572500000002, 21.07374465),
 (0.0, 0.0, 1.0)]
p.camera_position = cpos
p.window_size = 2000, 2000
p.show(use_ipyvtk=True)
p.screenshot("solar_9")
print(p.camera_position)

7.2.2. Visualize SVF

# 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_sky_acc_lattice * highres_solar_envelope_minmax
# 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["Sky view factor"] = base_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])
p.add_volume(grid, cmap="coolwarm", clim=[0.5, 1.0],opacity=opacity, shade=False)

# plotting
cpos = [(785.8704805788776, 708.4540755788776, 741.8613927288776),
 (65.08283250000001, -12.333572500000002, 21.07374465),
 (0.0, 0.0, 1.0)]
p.camera_position = cpos
p.window_size = 2000, 2000
p.show(use_ipyvtk=True)
p.screenshot("solar_10")
print(p.camera_position)

7. Save Sun Access Lattice into a CSV

# save the sun access latice to csv

csv_path = os.path.relpath('../data/solar_accessibility_324.csv')
highres_sunacc_lattice.to_csv(csv_path)
csv_path = os.path.relpath('../data/solar_envelope_324.csv')
highres_solar_envelope_minmax.to_csv(csv_path)

# save the SVF latice to csv

csv_path = os.path.relpath('../data/SVF_324.csv')
highres_sky_acc_lattice.to_csv(csv_path)

Credits

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