09 - Terrain and surface reflectance¶
This chapter covers the three kinds of terrain (flat ground, raster elevation and triangle-mesh terrain) and the common surface reflectance models: Lambertian, RPV bidirectional reflectance and water.
Flat ground¶
Without elevation data the terrain is a plane at z = 0 covering scene.size:
import less
scene = less.Scene()
scene.size = 50.0
scene.terrain = less.Terrain(property=less.Lambertian(reflectance=0.15))
Raster elevation (DEM)¶
dem can be a GeoTIFF path or a 2-D NumPy array in metres. The elevation grid spans the horizontal extent of scene.size: the first row of the array is the north edge and the first column the west edge. The lowest elevation is subtracted on import, so the lowest point of the terrain is at z = 0.
Here an array describes a small hill:
import numpy as np
x = np.linspace(0, 100, 101)
gx, gy = np.meshgrid(x, x)
hill = 15.0 * np.exp(-((gx - 50) ** 2 + (gy - 50) ** 2) / (2 * 20.0 ** 2))
scene = less.Scene()
scene.size = 100.0
scene.terrain = less.Terrain(dem=hill, property=less.Lambertian(reflectance=0.2))
scene.illumination = less.Illumination(
source=less.Sun(zenith=40, azimuth=135), atmosphere=less.NoAtmosphere())
image = scene.simulate(less.OpticalImager(
less.Orthographic(image_size=128), bands=[550], quality=32))
With a file, write less.Terrain(dem="elevation.tif", property=...). A few points:
- the grid needs at least two rows and two columns of finite values; fill or crop no-data areas before import;
subsample=Nkeeps every N-th sample, which reduces the resolution of large DEMs;- each grid cell is split into two triangles, and
height_at()interpolates on those triangles, exactly matching the simulated surface.
print(scene.terrain.height_at(50.0, 50.0)) # hilltop
print(scene.terrain.height_at([10.0, 50.0], [10.0, 50.0])) # arrays work too
To place objects on a DEM, give only (x, y); each object is put on the ground at that point:
tree = less.Object("ash", mesh=less.asset_path("FREX.obj"))
tree.set_property("leaves", less.Prospect(cab=40))
tree.set_property("stem_branch", less.Lambertian(reflectance=0.1))
scene.add(tree, positions=[[30.0, 50.0], [50.0, 50.0], [70.0, 50.0]])
Triangle-mesh terrain¶
Complex terrain, cliffs and man-made surfaces can be described by an OBJ mesh:
scene = less.Scene()
scene.size = (100, 100)
scene.terrain = less.Terrain(
mesh=less.asset_path("tutorial_terrain.obj"),
property=less.Lambertian(reflectance=0.2),
)
scene.illumination = less.Illumination(
source=less.Sun(zenith=30), atmosphere=less.NoAtmosphere())
image = scene.simulate(less.OpticalImager(
less.Orthographic(image_size=128), bands=[550, 850], quality=32))
Mesh terrain is read with the OBJ coordinate convention (see chapter 5) and keeps its size, position and elevation; it is not stretched to scene.size, so make sure the mesh lies within the area you observe. scene.terrain.height_at(x, y) returns the height of the topmost surface at that point.
Use either dem or mesh, not both. To change the terrain, assign scene.terrain again. Uneven terrain (DEM or mesh) cannot be combined with periodic boundaries (scene.repetitive).
Lambertian reflection¶
A Lambertian surface reflects the same radiance in every direction; it is the most common model for soils and man-made surfaces:
soil = less.Lambertian(reflectance=0.15)
soil = less.Lambertian(
wavelengths=[400, 500, 600, 700, 800, 900, 1000],
reflectance=[0.03, 0.05, 0.08, 0.10, 0.25, 0.30, 0.32],
)
With transmittance it describes a translucent sheet, and with back_reflectance the two faces reflect differently, which is enough for simple leaves (more realistic leaf models are in chapter 10):
RPV bidirectional reflectance¶
The RPV (Rahman–Pinty–Verstraete) model describes how reflectance changes with viewing direction. It is useful when a parameterised surface stands in for a canopy or rough soil:
| Parameter | Meaning |
|---|---|
rho0 |
Overall level of reflectance |
k |
Below 1, oblique directions are brighter (bowl shape); above 1, nadir is brighter (bell shape) |
theta |
Below 0, backscattering is stronger |
rhoc |
Hotspot strength |
Each parameter may be a spectral array. Chapter 21 shows how to simulate and plot bidirectional reflectance.
Water¶
WaterSurface treats the terrain as the mean water surface and describes wind-driven waves with the Cox–Munk slope distribution. It computes Fresnel specular reflection and suits sun glint and multi-angle water reflectance:
water = less.WaterSurface(
roughness=less.CoxMunk(wind_speed=5.0, wind_azimuth=90.0),
refractive_index=1.34,
)
scene = less.Scene()
scene.size = 100.0
scene.terrain = less.Terrain(property=water)
scene.illumination = less.Illumination(
source=less.Sun(zenith=30, azimuth=180), atmosphere=less.NoAtmosphere())
brf = scene.simulate(less.BRFSensor(
bands=[443, 560, 665],
directions=[(0, 0), (30, 0), (30, 90), (30, 180)],
num_photons=200_000,
))
print(brf.brf) # one row per direction, one column per band
wind_speed is the wind speed at 10 m, valid from 1 to 14 m s⁻¹; wind_azimuth is measured clockwise from north.
To include light coming out of the water, add a water-leaving term. With known remote-sensing reflectance Rrs (sr⁻¹):
water = less.WaterSurface(
roughness=less.CoxMunk(wind_speed=5.0),
water_leaving=less.PrescribedRrs(
wavelengths=[443, 560, 665], values=[0.0060, 0.0030, 0.0008]),
)
With known absorption and backscattering coefficients (m⁻¹):
water = less.WaterSurface(
roughness=less.CoxMunk(wind_speed=5.0),
water_leaving=less.DeepWaterIOP(
wavelengths=[443, 560, 665],
absorption=[0.025, 0.070, 0.430],
backscattering=[0.006, 0.004, 0.002]),
)
The water model is for optically deep water: it does not trace light below the surface and does not build wave geometry, so it cannot be used to study shallow bottoms, submerged targets or individual wave crests. It works with optical imaging, multi-angle reflectance and the 3-D radiation field.