BitmapMaps doc

No steep snow

  1. How we will remove snow
  2. Customizing functions from BitmapMaps
    1. Original
    2. Customized
  3. Optimize new parameters

How we will remove snow

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:

In other words, we will keep func_directional_pallette and simply write a replacement for func_reflection_coefficient and its callee shade_exponent.

Customizing functions from BitmapMaps

Original

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
end

Customized

We 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
end

Optimize new parameters

Let'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°.

Video: Varying two parameters

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.
Previous page Next page

CC BY-SA 4.0 hustf. Last modified: April 18, 2025. Website built with Franklin.jl and the Julia programming language.