More work

Demeton.Console/Properties/launchSettings.json

@@ -2,7 +2,7 @@
   "profiles": {
     "Demeton.Console": {
       "commandName": "Project",
-      "commandLineArgs": "dem-with-water-bodies N45E005"
+      "commandLineArgs": "dem-with-water-bodies N46E006"
     }
   }
 }

Demeton.Tests/Commands tests/DemWithWaterBodiesCommand/Downsampling water bodies.fs

@@ -1,7 +1,7 @@
 module Demeton.Tests.Commands_tests.DemWithWaterBodiesCommand.Downsampling_water_bodies
 
 
-open Demeton.WorldCover.Funcs
+open Demeton.WaterBodies.Funcs
 open Xunit
 open Swensen.Unquote
 

Demeton.Tests/Water bodies/Load water bodies data.fs

@@ -15,6 +15,7 @@ open Demeton.Dem.Funcs
 open FileSys
 open Demeton.WorldCover.Fetch
 open Demeton.WorldCover.Funcs
+open Demeton.WaterBodies.Types
 open Demeton.WaterBodies.Png
 
 

Demeton/Bnry.fs

@@ -115,3 +115,8 @@ let readBigEndianUInt32 (stream: Stream) : uint32 =
     ||| (((uint32) (readByte stream)) <<< 16)
     ||| (((uint32) (readByte stream)) <<< 8)
     ||| (((uint32) (readByte stream)))
+
+let writeBigEndianInt16 (value: int16) (stream: Stream) : Stream =
+    stream
+    |> writeByte ((byte) (value >>> 8))
+    |> writeByte ((byte) value)

Demeton/Commands/DemWithWaterBodiesCommand.fs

@@ -8,13 +8,15 @@ open Demeton.Dem.Types
 open Demeton.Dem.Funcs
 open Demeton.Aw3d.Types
 open Demeton.Aw3d.Funcs
+open Demeton.WaterBodies.Funcs
 open Demeton.WorldCover.Types
 open Demeton.WorldCover.Fetch
 open Demeton.WorldCover.Funcs
 open Demeton.WorldCover.Coloring
 open Demeton.WorldCover.WaterBodiesFetch
 open Raster
-
+open FileSys
+open System.IO
 
 type Options =
     { TileId: DemTileId
@@ -118,6 +120,10 @@ let fetchAw3dTile
     |> Result.bind (fun _ -> (readAw3dTile localCacheDir tileId) |> Result.Ok)
 
 
+let downsampleAw3dTile demResolutionParameter heightsArray =
+    let downsamplingFactor = float demResolutionParameter / 90.
+
+    heightsArray |> downsampleHeightsArray downsamplingFactor
 
 
 // todo 50: in order to properly identify the water bodies, we need to
@@ -148,6 +154,18 @@ let fetchWorldCoverTile
     readWorldCoverTiffFile localCacheDir (Some tile1by1Rect) containingTileId
 
 
+let downsampleWaterBodiesHeightsArray demResolutionParameter heightsArray =
+    // first calculate how many pixels per degree will the final DEM have
+    let finalDemPixelsPerDegree =
+        float demResolutionParameter / 90. * (float Aw3dTileSize)
+
+    // now, based on that value calculate the downsampling factor needed to
+    // reduce WorldCover's 12000 degrees per degree to the final resolution
+    let downsamplingFactor =
+        finalDemPixelsPerDegree / (float WorldCoverCellsPerDegree)
+
+    heightsArray |> downsampleWaterBodiesHeightsArray downsamplingFactor
+
 let identifyAndSimplifyWaterBodies
     demResolutionParameter
     worldCoverHeightsArray
@@ -171,6 +189,51 @@ let identifyAndSimplifyWaterBodies
     worldCoverTargetResolution |> colorWaterBodies
 
 
+let mergeDemWithWaterBodies
+    (waterBodiesHeightsArray: HeightsArray)
+    (demHeightsArray: HeightsArray)
+    : HeightsArray =
+    if
+        waterBodiesHeightsArray.Width <> demHeightsArray.Width
+        || waterBodiesHeightsArray.Height <> demHeightsArray.Height
+    then
+        invalidArg
+            "waterBodiesHeightsArray"
+            (sprintf
+                "The water bodies array (%d,%d) must have the same dimensions as the DEM array (%d,%d)."
+                waterBodiesHeightsArray.Width
+                waterBodiesHeightsArray.Height
+                demHeightsArray.Width
+                demHeightsArray.Height)
+
+    for index in 0 .. waterBodiesHeightsArray.Cells.Length - 1 do
+        if waterBodiesHeightsArray.Cells.[index] = 1s then
+            demHeightsArray.Cells.[index] <- Int16.MinValue
+
+    demHeightsArray
+
+
+let writeHeightsArrayToHgtFile
+    (fileName: FileName)
+    (heightsArray: HeightsArray)
+    : FileName =
+    Log.debug "Writing the DEM tile to %s..." fileName
+
+    fileName |> Pth.directory |> FileSys.ensureDirectoryExists |> ignore
+
+    FileSys.openFileToWrite (fileName)
+    |> function
+        | Ok stream ->
+            use stream = stream
+
+            for height in heightsArray.Cells do
+                stream |> Bnry.writeBigEndianInt16 height |> ignore
+
+            stream.Close()
+
+            fileName
+        | Error error -> raise error.Exception
+
 let run (options: Options) : Result<unit, string> =
     fetchAw3dTile options.TileId options.LocalCacheDir
     |> Result.map (fun aw3dTile ->
@@ -182,25 +245,45 @@ let run (options: Options) : Result<unit, string> =
             |> listAllAvailableFiles FileSys.openFileToRead
             |> Set.ofSeq
 
-        let result =
-            loadWaterBodiesTileFromCache
-                options.LocalCacheDir
-                availableWorldCoverTiles
-                options.TileId
-            |> makeNoneFileIfNeeded options.LocalCacheDir
-            |> extractWaterBodiesTileFromWorldCoverTileIfNeeded
-                options.LocalCacheDir
-
-        // todo 10: continue with the command - simplify water bodies
-        //   and then merge it with the DEM and save it as DEM tile
+        let downsampledAw3dTile =
+            aw3dTile |> downsampleAw3dTile options.DemResolution
+
+        loadWaterBodiesTileFromCache
+            options.LocalCacheDir
+            availableWorldCoverTiles
+            options.TileId
+        |> makeNoneFileIfNeeded options.LocalCacheDir
+        |> extractWaterBodiesTileFromWorldCoverTileIfNeeded
+            options.LocalCacheDir
+        |> function
+            | CachedTileLoaded tileHeightsArray ->
+                // todo 20: simplify water bodies
+
+                tileHeightsArray
+                |> Option.map (fun tileHeightsArray ->
+                    let downsampledWaterBodiesTile =
+                        tileHeightsArray
+                        |> downsampleWaterBodiesHeightsArray
+                            options.DemResolution
+
+                    let hgtFileName =
+                        downsampledAw3dTile
+                        |> mergeDemWithWaterBodies downsampledWaterBodiesTile
+                        |> writeHeightsArrayToHgtFile (
+                            Path.Combine(
+                                options.OutputDir,
+                                toTileName options.TileId + ".hgt"
+                            )
+                        )
+
+                    hgtFileName)
+            | TileNeedsToBeDownloaded _ ->
+                invalidOp "Bug: this should never happen"
+            | TileDoesNotExistInWorldCover _ ->
+                NotImplementedException(
+                    "The case when there is no WorldCover tile"
+                )
+                |> raise
+        |> ignore
 
         ())
-// match result with
-// | CachedTileLoaded heightsArray ->
-//     test <@ heightsArray |> Option.isSome @>
-// | _ -> Should.fail "Unexpected result"
-//
-//
-// fetchWorldCoverTile options.TileId options.LocalCacheDir
-// |> identifyAndSimplifyWaterBodies options.DemResolution
-// |> ignore)

Demeton/Dem/Funcs.fs

@@ -547,3 +547,58 @@ let heightsArrayToGrayscale16BitImageData
         heightsArray.Width
         heightsArray.Height
         (Grayscale16Bit.ImageDataInitializer1D initializer)
+
+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
+
+                downsampledHeightsArray.setHeightAtLocal (x, y) downsampledValue
+
+    downsampledHeightsArray