Ice can accumulate from air, while snow slides off slopes larger than 45-50°. In reality couloirs, north faces, wind exposure, surface roughness, terrain traps and nearby snow cover plays a role.
'Stetinden' is special in having just about only convex and flat faces, so we can probably get far with:
Use steepness as the only determining factor.
Affect reflection only, i.e. use the same hue for snow and stone.
In other words, we will keep func_directional_pallette and simply write a replacement for func_reflection_coefficient and its callee shade_exponent.
Previously, we adapted topo_relief. This time, we also need to adapt one of its callees (and that one's callee again). The original callees are:
"""
func_reflection_coefficient(; sun_deg = 202, light_elev_deg = [9, 30, 30, 30])
---> generic function (direction_no, z, n_ew, n_sn, n_up) ---> Float32
"""
function func_reflection_coefficient(; sun_deg = 202, light_elev_deg = [9, 30, 30, 30])
# Direction no: From sun, opposite sun, one side 60°, other side 60 °
Δazim_deg = [0, -180 - 60, - 180, - 180 + 60]
light_azim_deg = sun_deg .+ Δazim_deg
light_azim = light_azim_deg .* (π / 180)
light_elev = light_elev_deg .* (π / 180)
#
# Convert light azimuth and light elevation to x-y-z light unit vectors.
# Azimuth of 0 <=> Light from north, vector points north
# Azimuth of π / 2 <=> Light from east, vector points east
# Azimuth of π <=> Light from south, vector points south
# Elevation 0 <=> Light from horizon
# Elevation π / 2 <=> Light from above
#
# Unit vector components of all the light sources
l_ew = cos.(light_elev) .* sin.(light_azim)
l_sn = cos.(light_elev) .* cos.(light_azim)
l_up = sin.(light_elev)
# Closure on the light source vectors
f = let l_ew = l_ew, l_sn = l_sn, l_up = l_up
(dno, z, n_ew, n_sn, n_up) -> begin
@assert 1 <= dno <= 4
# The dot product of light and surface normal is the fraction of light
# reflected towards the observer. This is 'lambert shade'.
#
# We modify the lambert shading where there is snow.
# Lambert_reflection^shade_exponent(z, dno) naively does not exceed 1.0.
# Thrust, but check.
min(1.0f0, convert(Float32,
lambert_shade(n_ew, n_sn, n_up, l_ew[dno], l_sn[dno], l_up[dno] ^
shade_exponent(z, dno))))
end
end
end
"""
shade_exponent(z, direction_no)
---> Float64
Hard-coded values are:
Elevation z <= 400 --> 1.2
400 < z < 500 --> interpolated
500 <= z --> 0.4 for light sources 1 to 3
1.2 for light source 3 (opposite sun)
When used as an exponent to Lambertian shading, this means:
Elevation z <= 400 --> Slightly blank surface, reflects a little less all light angles considered
400 < z < 500 --> Transition
500 <= z --> For light source 1-3: Matte appearance, reflects more all light angles considered
For light source 3 (opposite sun), same for all elevations
"""
function shade_exponent(z, direction_no)
z1 = 400
z2 = 500
y1 = 1.2
y2 = direction_no == 3 ? y1 : 0.4
if z <= z1
y1
elseif z < z2
linterp(y1, y2, z1, z2, z)
else
y2
end
endWe define new functions outside of BitmapMaps, with new names and additional parameters.
The interesting change is in my_shade_exponent, but we need to change the callers too. For brevity, we drop repeating some commentary:
using BitmapMaps: func_directional_pallette, func_reflection_coefficient,
topo_relief, cell_to_utm_factor, full_folder_path, TOPORELIEF_FNAM,
lambert_shade, linterp
function my_steep_topo_relief(sb::SheetBuilder;
sun_deg = 156,
snow_lim_deg = 45,
exponent_factor = 15.0,
fna = "Toporelief_" * lpad(string(snow_lim_deg), 2, '0') * "°_" *
"exp_" * lpad(string(Int(round(exponent_factor))), 2, '0') *
".png")
#
ffna_source = joinpath(full_folder_path(sb), TOPORELIEF_FNAM)
if isfile(ffna_source)
rm(ffna_source; force = true)
end
f_hypso = func_directional_pallette()
f_reflect = my_func_refl(; sun_deg, snow_lim_deg, exponent_factor)
topo_relief(full_folder_path(sb), sb.cell_iter, cell_to_utm_factor(sb), f_hypso, f_reflect)
@show fna
ffna_dest = joinpath(full_folder_path(sb), fna)
cp(ffna_source, ffna_dest; force = true)
end
function my_func_refl(; sun_deg = 156,
light_elev_deg = [9, 30, 30, 30],
snow_lim_deg = 45,
exponent_factor = 15.0)
#
# Direction no: From sun, opposite sun, one side 60°, other side 60 °
Δazim_deg = [0, -180 - 60, - 180, - 180 + 60]
light_azim_deg = sun_deg .+ Δazim_deg
light_azim = light_azim_deg .* (π / 180)
light_elev = light_elev_deg .* (π / 180)
#
# Convert light azimuth and light elevation to x-y-z light unit vectors.
# Azimuth of 0 <=> Light from north, vector points north
# Azimuth of π / 2 <=> Light from east, vector points east
# Azimuth of π <=> Light from south, vector points south
# Elevation 0 <=> Light from horizon
# Elevation π / 2 <=> Light from above
#
# Unit vector components of all the light sources
l_ew = cos.(light_elev) .* sin.(light_azim)
l_sn = cos.(light_elev) .* cos.(light_azim)
l_up = sin.(light_elev)
# Save some time by doing this part just once
cos_lim = cos(snow_lim_deg * (π / 180))
# Closure on the light source vectors and snow limit
f = let l_ew = l_ew, l_sn = l_sn, l_up = l_up, cos_lim = cos_lim
(dno, z, n_ew, n_sn, n_up) -> begin
@assert 1 <= dno <= 4
# The dot product of light and surface normal is the fraction of light
# reflected towards the observer. This is 'lambert shade'.
#
# We modify the lambert shading where there is snow.
# Lambert_reflection^shade_exponent
# naively does not exceed 1.0.
# Thrust, but check.
min(1.0f0, convert(Float32,
lambert_shade(n_ew, n_sn, n_up, l_ew[dno], l_sn[dno], l_up[dno]) ^
my_shade_exponent(z, dno, n_up, cos_lim, exponent_factor)))
end
end
end
function my_shade_exponent(z, direction_no, n_up, cos_lim, exponent_factor)
z1 = 400
z2 = 500
y1 = 1.2
y2 = direction_no == 3 ? y1 : 0.4
if z <= z1
ex = y1
elseif z < z2
ex = linterp(y1, y2, z1, z2, z)
else
ex = y2
end
# Here comes the part about removing snow in steep places
if z > 450
# This is a faster way to say:
# acos(n_up) > snow_lim_deg * (π / 180)
if n_up < cos_lim
# A very high exponent makes light reflect in one direction only,
# in effect making the surface dark.
ex *= exponent_factor
end
end
ex
endLet's try that out with different steepness limits for snow cover:
for snow_lim_deg in 90:-1:0
my_steep_topo_relief(sb; snow_lim_deg)
end;
for exponent_factor in 1:25
my_steep_topo_relief(sb; exponent_factor)
end;
# Cleanup: restore the original `Toporelief.png`:
rm(joinpath(full_folder_path(sb), TOPORELIEF_FNAM))
topo_relief(sb);Here is a downscaled movie[1] of the images we just created, sorted by snow steepness angle. The exponent factor varies at snow steepness limit 45°.
Having selected snow steepness 45° and exponent factor 15, we conclude and show the composite map on the next page.
| [1] | Video made with Luxor / Cairo and FFMpeg. |