mibrammall.dev

Solar Shading from Lidar Data

General Setup

Setup of Libraries

WhiteBox Tools is a library used for geospatial data analysis. It can process lidar data, which has become popular due to self-driving cars, but we will use data obtained from aeroplane-mounted lidar. In what follows, I put the WhiteBox tools into the Lib folder of Python and everything worked just fine. The Python interface to WhiteBox Tools allows the executables to be separate and the location of the tools can be set from Python, but I didn't see a point to this for playing around with the toolbox.

We make use of pvlib to find the position of the sun in the sky at a given time. Also we will use pandas which is probably more familiar to Python users.

Setup of WhiteBox Tools in Python

This was super-simple, and we just need to create an instance of the class that we import from the library. We also need to set the working directory so that we can load our data. This has to be the complete path to the data directory, as it appears the tool does not accept relative paths, at least from my tests on Windows.

from WBT.whitebox_tools import WhiteboxTools
wbt = WhiteboxTools()

Since we're running a Jupyter notebook, we can't use the 'file' global to figure out what directory the script is running from, so we use path.abspath from Python's os standard library.

import os
d = os.path.abspath('')
wbt.set_working_dir(f'{d}/data')
wbt.set_verbose_mode(False)

The Shading Calculation

Generating a Digital Elevation Model

output_dem = "hyatt.tif"

The first step was to convert the lidar to a digital elevation model. This is simple using the WhiteBox Tools, and is a single function call. This, as all calls in the library, launches a subprocess through the Python interface. Managing this through Python is still nicer than manipulating the tool through the terminal, so I am grateful for the authors adding it.

wbt.lidar_idw_interpolation(
    i="hyatt.las",
    output=output_dem,
    parameter="elevation",
    returns="last",
    resolution=1.5,
    weight=2.0,
    radius=2.5,
)

And that's all there is to it

Calculating the Solar Position

import pandas as pd
from pvlib import solarposition

The solarposition function from pvlid needs some timezone and geographical location data. We could automate this to make less manual look-up necessary, but this is not the subject of this post and is quite involved. Instead, I manually searched and copied the data across.

Now we get a bunch of solar positions for a single day at the beginning of June. The WhiteBox tools writes data to disk and the images it produces can be quite large depending on the size of your input file, and a single day was enough to demonstrate the method.

tz = 'America/Kentucky/Louisville'
lat, lon = 35.5175, 86.5804

times = pd.date_range(
    '2019-06-01 00:00:00',
    '2019-06-02',
    closed='left',
    freq='H',
    tz=tz
)

solpos = solarposition.get_solarposition(
    times,
    lat,
    lon
)

There's little point in calculating the shading of something at nighttime, so we only use the times that the Sun is up, or above the horizon.

positions = solpos['apparent_elevation'] > 0]

Calculating the Hillshade

First, a wrapper function over WhiteBox Tools' hillshade function. It takes a filename to write the output to, and the azimuth and altitude of the Sun. The azimuth and altitude of the Sun at a given time has been calculated for us by pvlib in the previous section.

def draw_hillshade(
    output_file,
    azimuth,
    altitude
):
    wbt.hillshade(
        output_dem, 
        output_file,
        azimuth=azimuth,
        altitude=altitude
    )

def get_date_as_string(date):
    return date
        .tz_convert(None)
        .strftime("%Y-%b-%d-%H")

for date, row in positions.iterrows():
    azimuth = row["azimuth"]
    altitude = row["elevation"]
    date_str = get_date_as_string(date)
    output_file = f'{date_str}.tif'
    draw_hillshade(
        output_file,
        azimuth,
        altitude
    )

Making a GIF

Since we have a sequence of images, we could animate them to get a sense of how the light changes through the day as the Sun moves through the sky. We go ahead and import PIL and create an empty array to hold all the frames of the GIF that we will create.

We iterate over the same dataframe and read the images in to an image object using PIL, then we append to the array. Since the images are grayscale, it was necessary to convert the images into 'Pallette' mode to enable PIL to save the image as a GIF.

Finally, we save the image as a GIF by grabbing the first image in the array and using PIL from there to do the heavy lifting for us.

from PIL import Image
images = []

for date, row in positions.iterrows():
    date_str = get_date_as_string(date)
    file_name = f"{d}/data/{date_str}.tif"
    im = Image.open(file_name).convert('P')
    images.append(im)

images[0].save(
    f'{d}/daily.gif',
    save_all=True,
    append_images=images[1:],
    optimize=False,
    duration=360,
    loop=0
)