Fill urban roughness gaps with the closure's bare-soil endpoint

src/Lands/Lands.jl

@@ -28,6 +28,7 @@ export AbstractLand,
        IsotropicFrontalArea, EmpiricalFrontalArea,
        UniformHeight, VariableHeight,
        urban_roughness, compute_aerodynamic_roughness!, aerodynamic_parameters,
+       fill_aerodynamic_roughness_gaps!,
        # Atmosphere-facing accessors
        surface_temperature, surface_saturation
 

src/Lands/roughness/urban_roughness_field.jl

@@ -35,7 +35,8 @@ $(TYPEDSIGNATURES)
 Momentum roughness length `ℓᵐ` and zero-plane displacement `d` (as `Field`s on the grid
 of `h`) from a mean building-height field `h` and a plan-area index field `λᵖ`,
 under `closure` (default [`MorphometricRoughness`](@ref)). Where `λᵖ → 0` the result
-reduces to a bare-soil roughness.
+reduces to a bare-soil roughness, and cells of invalid morphometry are `NaN`; fill them
+with [`fill_aerodynamic_roughness_gaps!`](@ref) before passing the pair to a flux closure.
 """
 function urban_roughness(h, λᵖ; closure = MorphometricRoughness(eltype(h.grid)))
     grid = h.grid
@@ -45,6 +46,55 @@ function urban_roughness(h, λᵖ; closure = MorphometricRoughness(eltype(h.grid
     return ℓᵐ, d
 end
 
+#####
+##### Gap filling
+#####
+
+# A similarity profile needs a finite, positive roughness length.
+@inline roughness_gap(ℓᵐ, d) = !(isfinite(ℓᵐ) & (ℓᵐ > 0)) | !isfinite(d)
+
+@kernel function _fill_aerodynamic_roughness_gaps!(ℓᵐ, d, ℓᵘ, dᵘ)
+    i, j = @index(Global, NTuple)
+    @inbounds begin
+        ℓᵐᵢⱼ = ℓᵐ[i, j, 1]
+        dᵢⱼ  = d[i, j, 1]
+        gap  = roughness_gap(ℓᵐᵢⱼ, dᵢⱼ)
+        ℓᵐ[i, j, 1] = ifelse(gap, ℓᵘ, ℓᵐᵢⱼ)
+        d[i, j, 1]  = ifelse(gap, dᵘ, dᵢⱼ)
+    end
+end
+
+"""
+$(TYPEDSIGNATURES)
+
+Replace, in place, the gaps left where `closure` could not evaluate a cell with the
+`unbuilt` pair `(ℓ, d)`, by default the closure's own bare-soil endpoint
+`aerodynamic_parameters(closure, 0, 0)`. A gap is a non-finite `ℓᵐ` or `d`, or a
+non-positive `ℓᵐ`, which a similarity profile cannot take either. The roughness length `ℓᵐ`
+and zero-plane displacement `d` are filled together, so a gap in either becomes unbuilt
+surface in both.
+
+GHSL publishes no tile for the all-ocean cells of its Mollweide grid, so the offshore part
+of a coastal window comes back as gaps, as does any no-data pixel inside a published tile.
+The fill makes those cells evaluable and no more: bare soil (`ℓᵐ = 0.03` m by default) is
+two orders of magnitude rougher than a water surface, and the atmosphere-land flux closure
+applies land similarity theory wherever the land component runs. Where the gaps are known
+to be open water, pass a water roughness such as `unbuilt = (1e-4, 0)`, or keep the land
+component off those cells. For gaps that fall inside built-up land, inpaint `λᵖ` and `h`
+with `DataWrangling.inpaint_mask!` before evaluating the closure instead, so the fill comes
+from horizontal neighbors rather than from a constant.
+"""
+function fill_aerodynamic_roughness_gaps!(ℓᵐ, d, closure::AbstractUrbanRoughness; unbuilt = aerodynamic_parameters(closure, 0, 0))
+    ℓᵘ, dᵘ = unbuilt
+    ℓᵘ > 0 || throw(ArgumentError("the unbuilt roughness length must be positive, got $ℓᵘ"))
+    dᵘ >= 0 || throw(ArgumentError("the unbuilt zero-plane displacement must be non-negative, got $dᵘ"))
+
+    grid = ℓᵐ.grid
+    FT = eltype(ℓᵐ)
+    launch!(architecture(grid), grid, :xy, _fill_aerodynamic_roughness_gaps!, ℓᵐ, d, convert(FT, ℓᵘ), convert(FT, dᵘ))
+    return ℓᵐ, d
+end
+
 #####
 ##### Measured-morphometry builder: per-cell σʰ, hᵐᵃˣ and λᶠ from a footprint-level dataset
 ##### instead of the closure's regressions and frontal-area estimator.
@@ -90,7 +140,8 @@ Momentum roughness length `ℓᵐ` and zero-plane displacement `d` (as `Field`s
 the fields a footprint-level dataset such as `GlobalBuildingFootprints3D` aggregates. The
 closure (default [`MorphometricRoughness`](@ref)) is fed the measured height heterogeneity
 in place of its frontal-area estimator and `σʰ`/`hᵐᵃˣ` regressions. Where `λᵖ → 0` the
-result reduces to a bare-soil roughness.
+result reduces to a bare-soil roughness. Cells the closure cannot evaluate are `NaN`; fill
+them with [`fill_aerodynamic_roughness_gaps!`](@ref) before passing the pair to a flux closure.
 """
 function urban_roughness(h, λᵖ, σʰ, hᵐᵃˣ, λᶠ; closure = MorphometricRoughness(eltype(h.grid)))
     grid = h.grid

src/NumericalEarth.jl

@@ -109,6 +109,7 @@ export
     IsotropicFrontalArea, EmpiricalFrontalArea,
     UniformHeight, VariableHeight,
     urban_roughness,
+    fill_aerodynamic_roughness_gaps!,
     surface_temperature,
     surface_layer_diagnostics,
     regrid_bathymetry,

test/test_urban_roughness.jl

@@ -362,3 +362,43 @@ for arch in test_architectures
         @test all(≈(dref), interior(d))
     end
 end
+
+@testset "Gap filling to an unbuilt surface" begin
+    closure = MorphometricRoughness()
+    grid = LatitudeLongitudeGrid(CPU(), Float64; size = (4, 4),
+                                 longitude = (-0.1, 0.1), latitude = (51.4, 51.6),
+                                 topology = (Bounded, Bounded, Flat))
+    Nx = size(grid, 1)
+    built_roughness, built_displacement = aerodynamic_parameters(closure, 0.3, 15.0)
+    ℓˢᵒⁱˡ, dˢᵒⁱˡ = aerodynamic_parameters(closure, 0, 0)
+
+    # A building dataset with a missing tile: the closure marks it NaN by design.
+    λᵖ = Field{Center, Center, Nothing}(grid)
+    h  = Field{Center, Center, Nothing}(grid)
+    set!(λᵖ, (λ, φ) -> ifelse(λ < 0, 0.3, NaN))
+    set!(h,  (λ, φ) -> ifelse(λ < 0, 15.0, NaN))
+    roughness_and_displacement() = urban_roughness(h, λᵖ; closure)
+
+    ℓᵐ, d = roughness_and_displacement()
+    @test any(isnan, interior(ℓᵐ))
+    fill_aerodynamic_roughness_gaps!(ℓᵐ, d, closure)
+    @test all(isfinite, interior(ℓᵐ)) && all(isfinite, interior(d))
+    @test ℓᵐ[1, 1, 1] ≈ built_roughness && d[1, 1, 1] ≈ built_displacement   # built cells are untouched
+    @test ℓᵐ[Nx, 1, 1] ≈ ℓˢᵒⁱˡ && d[Nx, 1, 1] ≈ dˢᵒⁱˡ                          # the gap becomes unbuilt surface
+
+    # The unbuilt surface can be prescribed, e.g. open water.
+    ℓᵐ, d = roughness_and_displacement()
+    fill_aerodynamic_roughness_gaps!(ℓᵐ, d, closure; unbuilt = (1e-4, 0))
+    @test ℓᵐ[Nx, 1, 1] ≈ 1e-4 && d[Nx, 1, 1] == 0
+    @test ℓᵐ[1, 1, 1] ≈ built_roughness
+
+    # A non-positive roughness is a gap too.
+    ℓᵐ, d = roughness_and_displacement()
+    ℓᵐ[1, 1, 1] = 0
+    fill_aerodynamic_roughness_gaps!(ℓᵐ, d, closure)
+    @test ℓᵐ[1, 1, 1] ≈ ℓˢᵒⁱˡ && d[1, 1, 1] == 0
+
+    # The unbuilt surface itself must be evaluable.
+    @test_throws ArgumentError fill_aerodynamic_roughness_gaps!(ℓᵐ, d, closure; unbuilt = (0, 0))
+    @test_throws ArgumentError fill_aerodynamic_roughness_gaps!(ℓᵐ, d, closure; unbuilt = (0.03, -1))
+end