23 - Irradiance maps¶
IrradianceMapSensor divides the scene horizontally into a regular grid and records the downward and upward radiant flux in each cell. It is used to analyse solar radiation over complex terrain, the distribution of photosynthetically active radiation under a canopy, and the exchange of radiation between the ground and the canopy.
What is recorded¶
Photons start at the light sources and travel through the scene. Each grid cell (pixel) is a vertical column with a downwelling and an upwelling accumulator:
- Every time a photon interacts with the scene, its current flux is added to the downwelling value of that pixel;
- when a photon travels from an interaction in one pixel to another pixel, the flux is added to the upwelling value of the first pixel; consecutive interactions within one pixel add nothing to the upwelling value;
- when a photon finally leaves the scene, its flux is added to the upwelling value of the pixel of its last interaction.
For statistical vegetation, a scattering event on a statistical leaf is counted in the pixel where it happens. A canopy can span any number of pixels.
Basic use¶
import numpy as np
import less
scene = less.Scene()
scene.size = 40.0
scene.terrain = less.Terrain(property=less.Lambertian(reflectance=0.15))
scene.illumination = less.Illumination(
source=less.Sun(zenith=30, azimuth=180), atmosphere=less.NoAtmosphere())
sensor = less.IrradianceMapSensor(
bands=[550, 850],
grid=(40, 40), # number of pixels; or resolution=1.0 (m)
photons_per_m2=400, # photons per square metre of scene
max_bounces=None, # limit on recorded interactions; None means no limit
)
product = scene.simulate(sensor)
product.downwelling # (rows, columns, bands) total flux per pixel column, W nm⁻¹
product.upwelling
product.downwelling_irradiance() # divided by the horizontal pixel area, W m⁻² nm⁻¹
product.wavelengths
product.cell_size # pixel size (dx, dy), m
| Parameter | Meaning |
|---|---|
bands |
Bands (nm) |
grid |
Number of pixels (nx, ny), default (100, 100) |
resolution |
Pixel size (m); when given, grid is ignored |
photons_per_m2 |
Photons per square metre, typically 400–4000 |
num_photons |
Total number of photons, instead of photons_per_m2 |
max_bounces |
Record only the first interactions; 1 records only directly arriving sun and sky radiation |
Converting to irradiance¶
downwelling and upwelling are the total flux of each pixel column. To convert to irradiance, choose an area:
surface= |
Divide by | Meaning |
|---|---|---|
"horizontal" (default) |
Horizontal pixel area | Irradiance on a horizontal surface, as in solar radiation maps |
"slope" |
Horizontal area / cos(slope) | Reading of a radiometer mounted on the slope |
"surface" |
Actual ground area inside the pixel | Normalised by the real surface area |
horizontal = product.downwelling_irradiance()
up = product.upwelling_irradiance("horizontal")
par = product.downwelling_integrated((400, 700), surface="horizontal") # broadband, W m⁻²
On uneven terrain the result also contains the slope terrain_slope (radians) and the surface area terrain_surface_area (m²) of every pixel; on flat ground they are None.
Check: flat Lambertian ground¶
Without skylight, the downwelling irradiance over flat ground equals the solar irradiance on a horizontal surface, and the upwelling equals that times the reflectance:
import math
horizontal_sun = np.array([1.6203, 0.7190]) # expected horizontal irradiance, W m⁻² nm⁻¹
scene.illumination = less.Illumination(
source=less.Sun(zenith=30, azimuth=180, wavelengths=[550, 850],
irradiance=horizontal_sun / math.cos(math.radians(30))),
atmosphere=less.NoAtmosphere())
check = scene.simulate(less.IrradianceMapSensor(bands=[550, 850], resolution=1.0, photons_per_m2=400))
print("down:", check.downwelling_irradiance().mean(axis=(0, 1))) # close to horizontal_sun
print("up: ", check.upwelling_irradiance().mean(axis=(0, 1))) # close to 0.15 × horizontal_sun
Sun.irradiance is the irradiance perpendicular to the sun's rays, hence the division by cos 30°.
Uneven terrain¶
from scipy.ndimage import gaussian_filter
rng = np.random.default_rng(42)
dem = np.zeros((128, 128))
for sigma, amplitude in [(18, 6), (9, 3), (4, 1.5), (2, 0.5)]:
dem += amplitude * gaussian_filter(rng.standard_normal((128, 128)), sigma=sigma)
dem = (dem - dem.min()) * 12.0 / (dem.max() - dem.min()) # 12 m of relief
hills = less.Scene()
hills.size = 100.0
hills.terrain = less.Terrain(property=less.Lambertian(reflectance=0.30), dem=dem.astype(np.float32))
hills.illumination = less.Illumination(
source=less.Sun(zenith=45, azimuth=180), atmosphere=less.NoAtmosphere())
result = hills.simulate(less.IrradianceMapSensor(bands=[550], resolution=0.78, photons_per_m2=2000))
on_horizontal = result.downwelling_irradiance("horizontal")[:, :, 0]
on_slope = result.downwelling_irradiance("slope")[:, :, 0]
on_surface = result.downwelling_irradiance("surface")[:, :, 0]
ratio = result.upwelling_irradiance("horizontal")[:, :, 0] / np.maximum(on_horizontal, 1e-10)
print("mean upwelling/downwelling:", np.nanmean(ratio))
Downwelling irradiance differs strongly between sunny slopes and shadows. Radiation intercepted again within a pixel changes the relation between upwelling and downwelling, so on uneven terrain and in canopies upwelling/downwelling is not the surface albedo, and downwelling minus upwelling is not the local absorption.
The sensor grid may differ from the elevation grid. The surface area inside a pixel is summed from the actual triangles, so the total area of a region does not depend on the grid. The slope is the arc cosine of horizontal projected area / 3-D area; pixels without ground have zero area and NaN slope.
Irradiance maps and the radiation field¶
| Irradiance map | Radiation field (chapter 31) | |
|---|---|---|
| Dimensions | 2-D grid | 3-D voxels or per scene element |
| Records | Flux entering and leaving each pixel column | Absorption in each voxel or element |
| Typical use | Ground radiation, terrain radiation | Canopy absorption, APAR, vertical profiles |
Companion script: scripts/23_irradiance_map.py