Misc

Demeton.Tests/Aw3d/AW3D experiments.fs

@@ -52,7 +52,7 @@ let xcTracerHillshader
         | true -> Rgba8Bit.TransparentColor
         | false ->
             let aspectDiff =
-                Demeton.Geometry.Common.differenceBetweenAngles
+                differenceBetweenAngles
                     aspect
                     parameters.SunAzimuth
                     (Math.PI * 2.)

Demeton.Tests/Commands tests/TileShadeCommand/Calculating geo area needed.fs

@@ -0,0 +1,31 @@
+module Tests.Commands_tests.TileShadeCommand.Calculating_geo_area_needed
+
+open Demeton.Commands
+
+open Xunit
+open Swensen.Unquote
+open TestHelp
+
+[<Fact>]
+let ``Geo area needed is calculated correctly``() =
+    let options: TileShadeCommand.Options = {
+        TileWidth = 800
+        TileHeight = 600
+        TileCenter = (4, 46)
+        PixelSize = None
+        MapScale = Some 250000
+        Dpi = 245
+        WaterBodiesColor = "#49C8FF" |> Png.Rgba8Bit.parseColorHexValue
+    }
+
+    match
+        TileShadeCommand.createProjection options
+        |> Result.bind (TileShadeCommand.calculateGeoAreaNeeded options) with
+    | Ok geoAreaNeeded ->
+        test <@ geoAreaNeeded.MinLon |> isApproxEqualTo 3.866 (Decimals 3)  @>
+        test <@ geoAreaNeeded.MinLat |> isApproxEqualTo 45.930 (Decimals 3)  @>
+        test <@ geoAreaNeeded.MaxLon |> isApproxEqualTo 4.134 (Decimals 3)  @>
+        test <@ geoAreaNeeded.MaxLat |> isApproxEqualTo 46.069 (Decimals 3)  @>
+
+    | Error error -> invalidOp error
+

Demeton.Tests/Commands tests/TileShadeCommand/Creating map projection.fs

@@ -0,0 +1,23 @@
+module Tests.Commands_tests.TileShadeCommand.Creating_map_projection
+
+open Demeton.Commands
+
+open Xunit
+
+[<Fact>]
+let ``Projection is created``() =
+    let options: TileShadeCommand.Options = {
+        TileWidth= 100
+        TileHeight = 200
+        TileCenter = (10, 20)
+        PixelSize= None
+        MapScale=  Some 250000
+        Dpi=245
+        WaterBodiesColor = "#49C8FF" |> Png.Rgba8Bit.parseColorHexValue
+    }
+
+    let projection = TileShadeCommand.createProjection options
+
+    match projection with
+    | Ok projection -> ()
+    | Error error -> invalidOp error

Demeton.Tests/Demeton.Tests.fsproj

@@ -102,13 +102,15 @@
         <Compile Include="Vectorization tests\Steps to moves transformation.fs"/>
         <Compile Include="Svg\SVG path tests.fs"/>
         <Compile Include="Svg\SVG XML writing tests.fs"/>
+        <Compile Include="Commands tests\Command line parameters parsing.fs"/>
         <Compile Include="Commands tests\ShadeCommand\Command line parsing.fs"/>
         <Compile Include="Commands tests\ShadeCommand\Running the command.fs"/>
         <Compile Include="Commands tests\ShadeCommand\Splitting into intervals.fs"/>
         <Compile Include="Commands tests\ShadeCommand\Generating shaded tile.fs"/>
         <Compile Include="Commands tests\ShadeCommand\Shading the raster.fs"/>
         <Compile Include="Commands tests\ShadeCommand\Saving the tile.fs"/>
-        <Compile Include="Commands tests\Command line parameters parsing.fs"/>
+        <Compile Include="Commands tests\TileShadeCommand\Creating map projection.fs" />
+        <Compile Include="Commands tests\TileShadeCommand\Calculating geo area needed.fs" />
         <Compile Include="Commands tests\ImportSrtmTilesCommand tests.fs"/>
         <Compile Include="Acceptance tests\Generating hillshading rasters.fs"/>
         <Compile Include="FParsec exploratory tests.fs"/>

Demeton/Aw3d/Funcs.fs

@@ -140,6 +140,13 @@ let ensureAw3dTiles
     : Result<DemTileId list, string> =
     Log.info "Ensuring all needed AW3D tiles are there..."
 
+    Log.info
+        "Geo area needed: minLon: %f, minLat: %f, maxLon: %f, maxLat: %f"
+        bounds.MinLon
+        bounds.MinLat
+        bounds.MaxLon
+        bounds.MaxLat
+
     let aw3dTilesNeeded = bounds |> boundsToAw3dTiles |> Seq.toList
 
     let aw3dTileResults =

Demeton/Commands/ShadeCommand.fs

@@ -327,6 +327,82 @@ let splitIntoIntervals minValue maxValue intervalSize =
         let intervalMaxValue = min (intervalMinValue + intervalSize) maxValue
         (intervalIndex, intervalMinValue, intervalMaxValue))
 
+
+/// <summary>
+/// Based on the coverage points specified in the options, calculate the
+/// minimum bounding rectangle of the points projected onto the map.
+/// This effectively calculates the total raster area that needs to be shaded.
+/// </summary>
+let calculateRasterMbr mapProjection options =
+    // project each coverage point
+    let projectedPoints =
+        options.CoveragePoints
+        |> List.map (fun (lonDegrees, latDegrees) ->
+            mapProjection.Proj (degToRad lonDegrees) (degToRad latDegrees))
+        |> List.filter Option.isSome
+        |> List.map (Option.get >> (fun (x, y) -> (x, -y)))
+
+    // calculate the minimum bounding rectangle of all the projected points
+    let rasterMbr = Bounds.mbrOf projectedPoints
+
+    // round off the raster so we work with integer coordinates
+    Rect.asMinMax
+        (int (floor rasterMbr.MinX))
+        (int (floor rasterMbr.MinY))
+        (int (ceil rasterMbr.MaxX))
+        (int (ceil rasterMbr.MaxY))
+
+
+/// <summary>
+/// Based on the map projection and the raster MBR, calculate the SRTM level
+/// needed for the shading process.
+/// </summary>
+let calculateSrtmLevelNeeded mapProjection rasterMbr =
+    minLonLatDelta rasterMbr mapProjection.Invert
+    |> lonLatDeltaToSrtmLevel 3600
+
+
+/// <summary>
+/// Calculates the list of raster tiles to generate based on the raster MBR and
+/// options.
+/// Also returns the maximum tile index (either on X or Y axes, whichever is
+/// larger).
+/// </summary>
+let constructRasterTilesList options (rasterMbr: Rect)  =
+    let tilesToGenerate =
+        [ for yIndex, tileMinY, tileMaxY in
+              splitIntoIntervals
+                  rasterMbr.MinY
+                  rasterMbr.MaxY
+                  options.TileSize do
+              for xIndex, tileMinX, tileMaxX in
+                  splitIntoIntervals
+                      rasterMbr.MinX
+                      rasterMbr.MaxX
+                      options.TileSize do
+                  let tileBounds =
+                      Rect.asMinMax tileMinX tileMinY tileMaxX tileMaxY
+
+                  yield (xIndex, yIndex, tileBounds) ]
+
+    let maxTileIndexX =
+        tilesToGenerate |> List.map (fun (xIndex, _, _) -> xIndex) |> List.max
+
+    let maxTileIndexY =
+        tilesToGenerate |> List.map (fun (_, yIndex, _) -> yIndex) |> List.max
+
+    let maxTileIndex = max maxTileIndexX maxTileIndexY
+
+    Log.info
+        "NOTE: The command will generate a total raster size of %dx%d pixels (%dx%d tiles)."
+        rasterMbr.Width
+        rasterMbr.Height
+        (maxTileIndexX + 1)
+        (maxTileIndexY + 1)
+
+    tilesToGenerate, maxTileIndex
+
+
 /// <summary>
 /// A function the generates a shaded raster tile.
 /// </summary>
@@ -490,96 +566,54 @@ let saveShadedRasterTile
             tilePngFileName)
         |> Result.mapError fileSysErrorMessage
 
+
+
+
 let runWithProjection
     mapProjection
     (options: Options)
     (generateTile: ShadedRasterTileGenerator)
     (saveTile: ShadedRasterTileSaver)
     : Result<unit, string> =
-
-    // project each coverage point
-    let projectedPoints =
-        options.CoveragePoints
-        |> List.map (fun (lonDegrees, latDegrees) ->
-            mapProjection.Proj (degToRad lonDegrees) (degToRad latDegrees))
-        |> List.filter Option.isSome
-        |> List.map (Option.get >> (fun (x, y) -> (x, -y)))
-
-    // calculate the minimum bounding rectangle of all the projected points
-    let rasterMbr = Bounds.mbrOf projectedPoints
-
-    // round off the raster so we work with integer coordinates
-    let rasterMbrRounded =
-        Rect.asMinMax
-            (int (floor rasterMbr.MinX))
-            (int (floor rasterMbr.MinY))
-            (int (ceil rasterMbr.MaxX))
-            (int (ceil rasterMbr.MaxY))
+    let rasterMbr = calculateRasterMbr mapProjection options
 
     // calculate SRTM level needed
-    let srtmLevel =
-        minLonLatDelta rasterMbrRounded mapProjection.Invert
-        |> lonLatDeltaToSrtmLevel 3600
+    let srtmLevel = calculateSrtmLevelNeeded mapProjection rasterMbr
 
     // then split it up into tiles
-    let tilesToGenerate =
-        [ for yIndex, tileMinY, tileMaxY in
-              splitIntoIntervals
-                  rasterMbrRounded.MinY
-                  rasterMbrRounded.MaxY
-                  options.TileSize do
-              for xIndex, tileMinX, tileMaxX in
-                  splitIntoIntervals
-                      rasterMbrRounded.MinX
-                      rasterMbrRounded.MaxX
-                      options.TileSize do
-                  let tileBounds =
-                      Rect.asMinMax tileMinX tileMinY tileMaxX tileMaxY
-
-                  yield (xIndex, yIndex, tileBounds) ]
-
-    let maxTileIndexX =
-        tilesToGenerate |> List.map (fun (xIndex, _, _) -> xIndex) |> List.max
+    let tilesToGenerate, maxTileIndex = constructRasterTilesList options rasterMbr
+
+    let generateAndSaveHillshadingTile
+        xIndex yIndex tileBounds tilesGeneratedSoFar =
+        Log.info $"Generating a shade tile %d{xIndex}/%d{yIndex}..."
+
+        generateTile
+            srtmLevel
+            tileBounds
+            options.RootShadingStep
+            mapProjection
+        |> Result.map (fun maybeGeneratedTile ->
+            match maybeGeneratedTile with
+            | Some imageData ->
+                Log.info "Saving the shade tile..."
+
+                saveTile
+                    options
+                    maxTileIndex
+                    (xIndex, yIndex)
+                    tileBounds
+                    imageData
+                |> ignore
 
-    let maxTileIndexY =
-        tilesToGenerate |> List.map (fun (_, yIndex, _) -> yIndex) |> List.max
+                tilesGeneratedSoFar + 1
+            | None -> tilesGeneratedSoFar)
 
-    let maxTileIndex = max maxTileIndexX maxTileIndexY
-
-    Log.info
-        "NOTE: The command will generate a total raster size of %dx%d pixels (%dx%d tiles)."
-        rasterMbrRounded.Width
-        rasterMbrRounded.Height
-        (maxTileIndexX + 1)
-        (maxTileIndexY + 1)
 
     tilesToGenerate
     |> List.fold
         (fun state (xIndex, yIndex, tileBounds) ->
             state
-            |> Result.bind (fun tilesGeneratedSoFar ->
-                Log.info $"Generating a shade tile %d{xIndex}/%d{yIndex}..."
-
-                generateTile
-                    srtmLevel
-                    tileBounds
-                    options.RootShadingStep
-                    mapProjection
-                |> Result.map (fun maybeGeneratedTile ->
-                    match maybeGeneratedTile with
-                    | Some imageData ->
-                        Log.info "Saving the shade tile..."
-
-                        saveTile
-                            options
-                            maxTileIndex
-                            (xIndex, yIndex)
-                            tileBounds
-                            imageData
-                        |> ignore
-
-                        tilesGeneratedSoFar + 1
-                    | None -> tilesGeneratedSoFar)))
+            |> Result.bind (generateAndSaveHillshadingTile xIndex yIndex tileBounds))
         (Ok 0)
     |> Result.map (fun actualTilesGenerated ->
         match actualTilesGenerated with