Fixed the downsampling bug

Demeton.Tests/Commands tests/DemWithWaterBodiesCommand/Dealing with DEM.fs

@@ -40,7 +40,7 @@ let ``Encoding water bodies info into DEM`` () =
         dem
         |> DemWithWaterBodiesCommand.encodeWaterBodiesInfoIntoDem waterBodies
 
-    let wb = Int16.MinValue
+    let wb = DemHeightNone
 
     test
         <@ encoded.Cells = [| 0s; 1s; 100s; 101s; -100s; -99s; -2s; -1s; 0s |] @>
@@ -56,5 +56,5 @@ let ``Extending heights array for SRTM HGT compatibility`` () =
         DemWithWaterBodiesCommand.extendHeightsArrayWithAdditionalRowAndColumn
             dem
 
-    let no = Int16.MinValue
+    let no = DemHeightNone
     test <@ extendedDem.Cells = [| 1s; 2s; no; 3s; 4s; no; no; no; no |] @>

Demeton.Tests/Commands tests/TileShadeCommand/Integration tests.fs

@@ -14,9 +14,9 @@ let ``run command`` () =
     let options: TileShadeCommand.Options =
         { TileWidth = 800
           TileHeight = 600
-          TileCenter = (7.17, 46)
-          HgtSize = 1200
-          MapScale = Some (TileShadeCommand.MapScaleOf 600000)
+          TileCenter = (7.17, 46.10)
+          HgtSize = 3600
+          MapScale = Some (TileShadeCommand.MapScaleOf 500000)
           Dpi = TileShadeCommand.DefaultDpi
           LambertHillshadingIntensity = 1.
           IgorHillshadingIntensity = 1.

Demeton.Tests/Dem/HeightArraysScenes.fs

@@ -42,7 +42,7 @@ let renderScene (heightsArray: HeightsArray) : string =
         for x in 0 .. heightsArray.Width - 1 do
             let cellChar =
                 match heightsArray.heightAtLocal (x, y) with
-                | Int16.MinValue -> " "
+                | DemHeightNone -> " "
                 | value -> ((value |> char) + '0' |> string)
 
             sb |> append cellChar |> ignore
@@ -67,7 +67,7 @@ let compareScenes (actualScene: string) (expectedScene: string) : string =
     else
         // Construct a new heights array that will contain just the cells
         // that differ between the two arrays. All other cells will have
-        // value Int16.MinValue. Also, record the information if there were
+        // value DemHeightNone. Also, record the information if there were
         // any actual differences.
 
         let anyDifferences = ref false
@@ -88,7 +88,7 @@ let compareScenes (actualScene: string) (expectedScene: string) : string =
                         anyDifferences := true
                         actualHeight
                     else
-                        DemHeight(Int16.MinValue))
+                        DemHeightNone)
             )
 
         // if there were no differences...

Demeton.Tests/Dem/Reading and writing HGT files.fs

@@ -109,7 +109,6 @@ let ``Can create heights array from SRTM heights sequence`` () =
     test <@ heights.Cells.[tileSize - 1] = bottomRightHeight @>
 
 [<Fact>]
-// [<Trait("Category", "slow")>]
 let ``Can read HGT file`` () =
     let tileSize = Demeton.Srtm.Funcs.SrtmTileSize
 
@@ -143,7 +142,7 @@ let ``Can read HGT file`` () =
 
 [<Fact>]
 let ``Can write and then read HGT data`` () =
-    let no = Int16.MinValue
+    let no = DemHeightNone
 
     let cells =
         [| 1s; 2s; 3s; no; 4s; 5s; 6s; no; 7s; 8s; 9s; no; no; no; no; no |]

Demeton.Tests/Shaders/Hillshading/Calculating aspect and slope.fs

@@ -38,7 +38,7 @@ let slopeAndAspectForHeights (heights3by3: float option [])  =
            Some(lon, lat - dx)
            Some(lon + dx, lat - dx) |]
 
-    match Hillshading.calculatePQ coords (Some heightsFlippedByYAxis) with
+    match Hillshading.calculatePQ coords heightsFlippedByYAxis with
     | Some(p, q) -> Hillshading.calculateSlopeAndAspect p q
     | _ -> invalidOp "bug"
 

Demeton.Tests/todo.fs

@@ -2,11 +2,6 @@
 
 module Demeton.Tests.todo
 
-// todo 5: the command does not work when --hgt-size is 1200
-//   - somewhere in the workflow 3600 is hardcoded
-
-// todo 10: fix the problem with the white cross in the hillshading
-
 // todo sometime 100: update shading docs now that we have added an array of fetchers
 
 // todo sometime 100: implement support for variable-resolution height arrays, not just

Demeton/Commands/DemWithWaterBodiesCommand.fs

@@ -201,14 +201,14 @@ let extendHeightsArrayWithAdditionalRowAndColumn (heightsArray: HeightsArray) =
 
         if (sourceIndex + 1) % oldWidth = 0 then
             // fill the final column with the "missing value" marker
-            extendedCells.[destIndex + 1] <- Int16.MinValue
+            extendedCells.[destIndex + 1] <- DemHeightNone
             destIndex <- destIndex + 2
         else
             destIndex <- destIndex + 1
 
     // fill the final row with the "missing value" marker
     while destIndex < extendedCells.Length do
-        extendedCells.[destIndex] <- Int16.MinValue
+        extendedCells.[destIndex] <- DemHeightNone
         destIndex <- destIndex + 1
 
     HeightsArray(
@@ -229,7 +229,7 @@ let ensureHgtFile cacheDir hgtSize tileId : HeightsArray option =
         )
 
     if FileSys.fileExists hgtFileName then
-        Log.info "The HGT file already exists at '%s'." hgtFileName
+        // Log.info "The HGT file already exists at '%s'." hgtFileName
 
         FileSys.openFileToRead hgtFileName
         |> function

Demeton/Dem/Funcs.fs

@@ -382,7 +382,7 @@ let merge
 /// <summary>
 /// Extracts a sub-array from the given array. If the sub-array is outside
 /// the bounds of the given array, the missing values are filled with
-/// Int16.MinValue.
+/// DemHeightNone.
 /// </summary>
 let extract (extractBounds: Rect) (heightArray: HeightsArray) : HeightsArray =
 
@@ -397,7 +397,7 @@ let extract (extractBounds: Rect) (heightArray: HeightsArray) : HeightsArray =
             heightArray.heightAt (x, y)
         else
             // this is not a water body pixel
-            Int16.MinValue
+            DemHeightNone
 
     HeightsArray(
         extractBounds.MinX,
@@ -554,58 +554,63 @@ let downsampleHeightsArray
     (factor: float)
     (heightsArray: HeightsArray)
     : HeightsArray =
-    let downsampledMinX = int (float heightsArray.MinX * factor)
-    let downsampledMinY = int (float heightsArray.MinY * factor)
-    let downsampledWidth = int (float heightsArray.Width * factor)
-    let downsampledHeight = int (float heightsArray.Height * factor)
-
-    let downsampledHeightsArray =
-        HeightsArray(
-            downsampledMinX,
-            downsampledMinY,
-            downsampledWidth,
-            downsampledHeight,
-            EmptyHeightsArray
-        )
-
-    for y in 0 .. downsampledHeight - 1 do
-        for x in 0 .. downsampledWidth - 1 do
-            // Calculate the exact original coordinates this new pixel maps to
-            let origStartX = (float x / factor)
-
-            let origEndX =
-                min (float (x + 1) / factor) (float (heightsArray.Width - 1))
-
-            let origStartY = (float y / factor)
-
-            let origEndY =
-                min (float (y + 1) / factor) (float (heightsArray.Height - 1))
-
-            let mutable weightedSum = 0.0
-            let mutable totalWeight = 0.0
-
-            // Iterate through the original pixels that overlap with the new pixel
-            for oy in int origStartY .. int origEndY do
-                for ox in int origStartX .. int origEndX do
-                    // Calculate the overlap area (weight) for each original pixel
-                    let left = max origStartX (float ox)
-                    let right = min origEndX (float ox + 1.0)
-                    let top = max origStartY (float oy)
-                    let bottom = min origEndY (float oy + 1.0)
-
-                    let weight = (right - left) * (bottom - top)
-
-                    let height = heightsArray.heightAtLocal (ox, oy)
-
-                    weightedSum <- weightedSum + (float height) * weight
-
-                    totalWeight <- totalWeight + weight
-
-            // Set the new pixel value based on the weighted average
-            if totalWeight > 0.0 then
-                let downsampledValue =
-                    Math.Round(weightedSum / totalWeight) |> int16
+    match factor with
+    | 1.0 -> heightsArray
+    | _ ->
+        let downsampledMinX = int (float heightsArray.MinX * factor)
+        let downsampledMinY = int (float heightsArray.MinY * factor)
+        let downsampledWidth = int (float heightsArray.Width * factor)
+        let downsampledHeight = int (float heightsArray.Height * factor)
 
-                downsampledHeightsArray.setHeightAtLocal (x, y) downsampledValue
+        let downsampledHeightsArray =
+            HeightsArray(
+                downsampledMinX,
+                downsampledMinY,
+                downsampledWidth,
+                downsampledHeight,
+                EmptyHeightsArray
+            )
 
-    downsampledHeightsArray
+        for y in 0 .. downsampledHeight - 1 do
+            for x in 0 .. downsampledWidth - 1 do
+                // Calculate the exact original coordinates this new pixel maps to
+                let origStartX = float x / factor
+                let origEndX = float (x + 1) / factor
+                let origStartY = float y / factor
+                let origEndY = float (y + 1) / factor
+
+                let mutable weightedSum = 0.0
+                let mutable totalWeight = 0.0
+
+                // Iterate through the original pixels that overlap with the new pixel
+                for oy in int origStartY .. int origEndY do
+                    if oy < heightsArray.Height then
+                        for ox in int origStartX .. int origEndX do
+                            if ox < heightsArray.Width then
+                                // Calculate the overlap area (weight) for each original pixel
+                                let left = max origStartX (float ox)
+                                let right = min origEndX (float ox + 1.0)
+                                let top = max origStartY (float oy)
+                                let bottom = min origEndY (float oy + 1.0)
+
+                                let weight = (right - left) * (bottom - top)
+
+                                if weight <> 0.0 then
+                                    let height =
+                                        heightsArray.heightAtLocal (ox, oy)
+
+                                    weightedSum <-
+                                        weightedSum + (float height) * weight
+
+                                    totalWeight <- totalWeight + weight
+
+                // Set the new pixel value based on the weighted average
+                if totalWeight > 0.0 then
+                    let downsampledValue =
+                        Math.Round(weightedSum / totalWeight) |> int16
+
+                    downsampledHeightsArray.setHeightAtLocal
+                        (x, y)
+                        downsampledValue
+
+        downsampledHeightsArray

Demeton/Shaders/Hillshading.fs

@@ -37,33 +37,27 @@ let gridSize (coords: LonLat option[]) =
 
     (width, height)
 
-let calculatePQ
-    (coords: LonLat option[])
-    (heightsMaybe: float option[] option)
-    =
-    match heightsMaybe with
-    | None -> None
-    | Some heights ->
-        let allHeightsAreAvailable = heights |> Array.forall Option.isSome
-
-        match allHeightsAreAvailable with
-        | false -> None
-        | true ->
-            let heights = heightsMaybe |> Option.get |> Array.map Option.get
-
-            let gridWidth, gridHeight = gridSize coords
-
-            let p =
-                ((heights.[8] + 2. * heights.[5] + heights.[2])
-                 - (heights.[6] + 2. * heights.[3] + heights.[0]))
-                / (8. * gridWidth)
-
-            let q =
-                ((heights.[8] + 2. * heights.[7] + heights.[6])
-                 - (heights.[2] + 2. * heights.[1] + heights.[0]))
-                / (8. * gridHeight)
-
-            Some(p, q)
+let calculatePQ (coords: LonLat option[]) (heights: float option[]) =
+    let allHeightsAreAvailable = heights |> Array.forall Option.isSome
+
+    match allHeightsAreAvailable with
+    | false -> None
+    | true ->
+        let heights = heights |> Array.map Option.get
+
+        let gridWidth, gridHeight = gridSize coords
+
+        let p =
+            ((heights.[8] + 2. * heights.[5] + heights.[2])
+             - (heights.[6] + 2. * heights.[3] + heights.[0]))
+            / (8. * gridWidth)
+
+        let q =
+            ((heights.[8] + 2. * heights.[7] + heights.[6])
+             - (heights.[2] + 2. * heights.[1] + heights.[0]))
+            / (8. * gridHeight)
+
+        Some(p, q)
 
 type SlopeAndAspect = float * float
 
@@ -140,10 +134,11 @@ let shadeRaster
 
                 let heights = neighborHeights neighborCoords
 
-                let pqMaybe = calculatePQ neighborCoords heights
+                heights
+                |> Option.bind (fun heights ->
+                    calculatePQ neighborCoords heights)
+                |> Option.map (fun (p, q) ->
 
-                match pqMaybe with
-                | Some(p, q) ->
                     let slope, aspect = calculateSlopeAndAspect p q
                     let height = (Option.get heights).[4] |> Option.get
 
@@ -156,9 +151,8 @@ let shadeRaster
                         // we flip the Y coordinate since DEM heights array
                         // is flipped vertically compared to the bitmap
                         (bitmapRect.MaxY - y - 1)
-                        pixelValue
-
-                | None -> ()
+                        pixelValue)
+                |> ignore
 
         Parallel.For(bitmapRect.MinY, bitmapRect.MaxY, processRasterLine)
         |> ignore

Demeton/Shaders/IgorHillshader.fs

@@ -28,7 +28,8 @@ let defaultParameters =
 let shadePixel parameters : Hillshading.PixelHillshader =
     fun _ slope aspect ->
         match Double.IsNaN(aspect) with
-        | true -> Rgba8Bit.TransparentColor
+        | true ->
+            Rgba8Bit.TransparentColor
         | false ->
             let aspectDiff =
                 differenceBetweenAngles