oceanbox/shaver
6
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

feat: working boundary detection

Browse Source
This commit is contained in:
Ole Tytlandsvik
2024-09-06 16:08:07 +02:00
parent 8b146b858c
commit 5d855a7620
6 changed files with 271 additions and 87 deletions
Download Patch File Download Diff File Expand all files Collapse all files

3
.gitignore vendored
View File

@@ -8,4 +8,5 @@ obj/
*.txt
.direnv/*
*.sln.DotSettings.user
delaunay
delaunay
.pre-commit-config.yaml

56
Tests/Tests.fs
View File

@@ -45,20 +45,60 @@ let coastlineTests =
testList "Coastline Tests" [
testCase "Minimal grid separate coastline"
<| fun () ->
let neighborIndex = Grid.makeNeighborIndex minimalGrid
// Expect all nodes are "coastal nodes"
let expected = ([| 0; 1; 2; 3 |], [||])
let res = separateCoastline neighborIndex minimalGrid
let expected =
([| (0, 1); (1, 2); (2, 3); (0, 3) |] |> Set.ofArray, [||])
let res =
separateCoastline minimalGrid
|> fun (a, b) -> (a |> Set.ofArray, b)
Expect.equal res expected "all should be part of coastline"
testCase "Grid with inner points separate coastline"
<| fun () ->
let neighborIndex = Grid.makeNeighborIndex withInnerPoint
// Inner point should not be a coastal node
let expected = ([| 0; 1; 2; 3; 4; 5; 6; 7 |], [| 8 |])
let res = separateCoastline neighborIndex withInnerPoint
let expected =
([|
(0, 1)
(1, 2)
(2, 3)
(3, 4)
(4, 5)
(5, 6)
(6, 7)
(0, 7)
|]
|> Set.ofArray,
[| 8 |])
let res =
separateCoastline withInnerPoint
|> fun (a, b) -> (a |> Set.ofArray, b)
Expect.equal res expected "inner point should be separated"
]
[<Tests>]
let groupPolyTest =
testList "Group by polygon Tests" [
testCase "Minimal grid group polygon"
<| fun () ->
let coastEdges = [| (0, 1); (1, 2); (2, 3); (0, 3) |]
let expected = [| [| 0; 1; 2; 3 |] |]
let res = buildPolygons coastEdges
Expect.equal res expected "polygon should be equal to coast nodes"
testCase "Grid with inner point group polygon"
<| fun () ->
let coastEdges = [|
(0, 1)
(1, 2)
(2, 3)
(3, 4)
(4, 5)
(5, 6)
(6, 7)
(0, 7)
|]
let expected = [| [| 0; 1; 2; 3; 4; 5; 6; 7 |] |]
let res = buildPolygons coastEdges
Expect.equal res expected "polygon should be equal to coast nodes"
]
[<EntryPoint>]
let main args = runTestsInAssemblyWithCLIArgs [] args

93
flake.lock generated
View File

@@ -1,5 +1,21 @@
{
"nodes": {
"flake-compat": {
"flake": false,
"locked": {
"lastModified": 1696426674,
"narHash": "sha256-kvjfFW7WAETZlt09AgDn1MrtKzP7t90Vf7vypd3OL1U=",
"owner": "edolstra",
"repo": "flake-compat",
"rev": "0f9255e01c2351cc7d116c072cb317785dd33b33",
"type": "github"
},
"original": {
"owner": "edolstra",
"repo": "flake-compat",
"type": "github"
}
},
"flake-utils": {
"inputs": {
"systems": "systems"
@@ -18,6 +34,27 @@
"type": "github"
}
},
"gitignore": {
"inputs": {
"nixpkgs": [
"pre-commit-hooks",
"nixpkgs"
]
},
"locked": {
"lastModified": 1709087332,
"narHash": "sha256-HG2cCnktfHsKV0s4XW83gU3F57gaTljL9KNSuG6bnQs=",
"owner": "hercules-ci",
"repo": "gitignore.nix",
"rev": "637db329424fd7e46cf4185293b9cc8c88c95394",
"type": "github"
},
"original": {
"owner": "hercules-ci",
"repo": "gitignore.nix",
"type": "github"
}
},
"nixpkgs": {
"locked": {
"lastModified": 1723362943,
@@ -34,10 +71,64 @@
"type": "github"
}
},
"nixpkgs-stable": {
"locked": {
"lastModified": 1720386169,
"narHash": "sha256-NGKVY4PjzwAa4upkGtAMz1npHGoRzWotlSnVlqI40mo=",
"owner": "NixOS",
"repo": "nixpkgs",
"rev": "194846768975b7ad2c4988bdb82572c00222c0d7",
"type": "github"
},
"original": {
"owner": "NixOS",
"ref": "nixos-24.05",
"repo": "nixpkgs",
"type": "github"
}
},
"nixpkgs_2": {
"locked": {
"lastModified": 1719082008,
"narHash": "sha256-jHJSUH619zBQ6WdC21fFAlDxHErKVDJ5fpN0Hgx4sjs=",
"owner": "NixOS",
"repo": "nixpkgs",
"rev": "9693852a2070b398ee123a329e68f0dab5526681",
"type": "github"
},
"original": {
"owner": "NixOS",
"ref": "nixpkgs-unstable",
"repo": "nixpkgs",
"type": "github"
}
},
"pre-commit-hooks": {
"inputs": {
"flake-compat": "flake-compat",
"gitignore": "gitignore",
"nixpkgs": "nixpkgs_2",
"nixpkgs-stable": "nixpkgs-stable"
},
"locked": {
"lastModified": 1724440431,
"narHash": "sha256-9etXEOUtzeMgqg1u0wp+EdwG7RpmrAZ2yX516bMj2aE=",
"owner": "cachix",
"repo": "pre-commit-hooks.nix",
"rev": "c8a54057aae480c56e28ef3e14e4960628ac495b",
"type": "github"
},
"original": {
"owner": "cachix",
"repo": "pre-commit-hooks.nix",
"type": "github"
}
},
"root": {
"inputs": {
"flake-utils": "flake-utils",
"nixpkgs": "nixpkgs"
"nixpkgs": "nixpkgs",
"pre-commit-hooks": "pre-commit-hooks"
}
},
"systems": {

23
flake.nix
View File

@@ -5,12 +5,17 @@
url = "github:numtide/flake-utils";
inputs.nixpkgs.follows = "nixpkgs";
};
pre-commit-hooks = {
url = "github:cachix/pre-commit-hooks.nix";
# inputs.nixpkgs.follows = "nixpgs";
};
};
outputs =
{
self,
nixpkgs,
flake-utils,
pre-commit-hooks,
...
}:
flake-utils.lib.eachDefaultSystem (
@@ -20,11 +25,29 @@
in
with pkgs;
{
checks = {
# On every git commit, format the F# code first
pre-commit-check = pre-commit-hooks.lib.${system}.run {
src = ./.;
hooks = {
fantomas = {
enable = true;
name = "fantomas";
description = "Format your F# code with fantomas.";
entry = "dotnet fantomas";
files = "(\\.fs$)|(\\.fsx$)";
};
};
};
};
devShells.default = mkShell {
inherit (self.checks.${system}.pre-commit-check) shellHook;
nativeBuildInputs = [
fsautocomplete
netcdf
dotnet-sdk_8
gdb
rustup
];
LD_LIBRARY_PATH = lib.makeLibraryPath [ netcdf ];
};

39
src/CLI/Program.fs
View File

@@ -14,6 +14,7 @@ type Arguments =
| [<AltCommandLine("-n")>] No_Reduce
| [<AltCommandLine("-q")>] Quiet
| [<AltCommandLine("-v")>] Verbose
| [<AltCommandLine("-m")>] Makai
interface IArgParserTemplate with
member arg.Usage =
@@ -24,6 +25,8 @@ type Arguments =
| No_Reduce -> "do not perform reduction on the grid."
| Quiet -> "disable all logging."
| Verbose -> "enable verbose logging."
| Makai ->
"write boundaries to makai compatible format to 'makai.txt'"
let configureSerilog level =
LoggerConfiguration()
@@ -63,6 +66,7 @@ let main argv =
// Whether to perform reduction. Defaults to true
let shouldReduce = args.Contains(No_Reduce) |> not
// Default to a filepath, but don't write to file if flag is not passed
let outPath =
if args.Contains(Output_Grid) then
@@ -71,43 +75,46 @@ let main argv =
|> Some
else
None
// Whether to write boundaries to makai format
let toMakai = args.Contains(Makai)
// Path of input file. Main argument and mandatory
let path = args.GetResult(Input_Grid)
Log.Information $"Using input file: %s{path}"
let ds = DataSetUri.Create path |> DataSet.Open
let grd = ds |> Fvcom.getGrid
let neighborIndex = grd.ToGrid() |> Grid.makeNeighborIndex
let coastNodes, innerNodes = separateCoastline neighborIndex (grd.ToGrid())
coastNodes
|> Array.tryFind (fun n ->
Map.find n neighborIndex.NodesAroundNode |> Array.length > 4
)
|> function
| Some n -> Log.Verbose $"Found non edge node: {n}"
| None -> Log.Verbose "No non edge node found"
// let neighborIndex = grd.ToGrid() |> Grid.makeNeighborIndex
let coastEdges, innerNodes = separateCoastline (grd.ToGrid())
Log.Verbose
$"Coast nodes: {coastNodes.Length} \t Inner nodes: {innerNodes.Length}"
$"Coast edges: {coastEdges.Length} \t Inner nodes: {innerNodes.Length}"
Log.Verbose "Polygons: "
let polygons = groupByPolygon coastNodes neighborIndex
let polygons = buildPolygons coastEdges
polygons |> Array.iter (fun n -> Log.Verbose $"{n.Length}")
// Assume longest polygon is boundary
let boundary = polygons |> Array.maxBy (_.Length)
let rest = polygons |> Array.filter (fun p -> p <> boundary)
let boundaryStrings = getPolygonString boundary grd
let restStrings =
rest |> Array.map (fun p -> getPolygonString p grd) |> Array.concat
File.WriteAllLines("boundary.txt", boundaryStrings)
File.WriteAllLines("islands.txt", restStrings)
if toMakai then
Log.Verbose "Writing boundaries to Makai format..."
let makaiStrings =
polygons
|> Array.sortBy _.Length
|> Array.map (fun p -> getPolygonStringMakai p grd)
|> Array.concat
File.WriteAllLines("makai.txt", makaiStrings)
let innerNodes =
if shouldReduce then
Log.Verbose "Reducing inner nodes..."
@@ -128,6 +135,10 @@ let main argv =
File.WriteAllLines("points.txt", innerNodeStrings)
// Extract each node index from coast edges
let coastNodes = coastEdges |> Array.map fst
// Write to ncdf file if output path is specified
match outPath with
| None -> ()
| Some p ->

144
src/Gridding/Polygon.fs
View File

@@ -4,90 +4,98 @@ open Oceanbox.FvcomKit
/// <summary>
/// Separate coastline nodes from inner nodes.
/// Returns (coastline nodes * inner nodes).
/// Returns (coastline edges * inner nodes).
/// </summary>
/// <param name="nIdx"></param>
/// <param name="grid"></param>
let separateCoastline (nIdx: Grid.NeighborIndex) (grid: Grid.Grid) =
let separateCoastline (grid: Grid.Grid) =
let getEdges (a, b, c) =
[ (a, b); (b, a); (a, c); (c, a); (b, c); (c, b) ] |> Set.ofList
let getEdges (a, b, c) = [ (a, b); (b, c); (c, a) ]
// Find whether an elem has less than two elems that it
// shares an edge with
// FIXME: This one is wrong
let isEdgeElem (elemIdx: Grid.ElemIdx) : bool =
let ownEdges = grid.Elem[elemIdx] |> getEdges
let allEdges =
(Map.empty, grid.Elem)
||> Array.fold (fun allEdges elem ->
// Get all edges of the elem and sort the coordinates
let elemEdges =
elem
|> getEdges
|> List.map (fun (a, b) -> if a < b then (a, b) else (b, a))
Grid.getElemsSurroundingElem nIdx grid elemIdx
|> Array.filter (fun surroundingIdx ->
// Avoid including self
if surroundingIdx = elemIdx then
false
else
let edges =
grid.Elem[surroundingIdx]
|> getEdges
edges = ownEdges
// Add edge count to map state
(allEdges, elemEdges)
||> List.fold (fun acc edge ->
match acc |> Map.tryFind edge with
| Some n -> acc |> Map.add edge (n + 1)
| None -> acc |> Map.add edge 1
)
)
|> Array.isEmpty
|> not
let cNodes =
([||], nIdx.ElemsAroundNode)
||> Map.fold (fun cNodes idx neighborElems ->
neighborElems
|> Array.tryFind (fun eIdx -> eIdx |> isEdgeElem)
|> function
| Some _ -> Array.append cNodes [| idx |]
| None -> cNodes
// Coastal edges will only have one reference
let coastalEdges =
(Array.empty, allEdges)
||> Map.fold (fun acc edge count ->
match count with
| 1 -> Array.append acc [| edge |]
| _ -> acc
)
|> Set.ofArray
let iNodes =
[| 0 .. (grid.Nodes.Length - 1) |]
|> Set.ofArray
|> fun allNodes -> Set.difference allNodes cNodes
// Find all nodes that are not part of a coastal edge
let innerNodes =
[| 0 .. grid.Nodes.Length - 1 |]
|> Array.filter (fun node ->
coastalEdges
|> Array.exists (fun (a, b) -> a = node || b = node)
|> not
)
(cNodes |> Array.ofSeq, iNodes |> Array.ofSeq)
(coastalEdges |> Array.ofSeq, innerNodes |> Array.ofSeq)
/// <summary>
/// Group given edge nodes by polygon
/// Build closed polygons using the given boundary edges
/// </summary>
/// <param name="nodes">Edge node indices</param>
/// <param name="nIdx"></param>
let groupByPolygon (nodes: Grid.NodeIdx[]) (nIdx: Grid.NeighborIndex) =
let nMap = nIdx.NodesAroundNode
/// <param name="edges"></param>
let buildPolygons (edges: (Grid.NodeIdx * Grid.NodeIdx) array) =
let isEdgeNode (node: Grid.NodeIdx) = nodes |> Array.contains node
// Find a single polygon given an edge node
/// Find a single polygon given an edge. Assume all polygons are closed
let rec findPolygon
(node: Grid.NodeIdx)
(acc: Set<Grid.NodeIdx>)
: Set<Grid.NodeIdx> =
Map.find node nMap
|> Array.tryFind (fun n ->
(n |> isEdgeNode) && (not <| (Set.contains n acc))
(edge: Grid.NodeIdx * Grid.NodeIdx)
(acc: Grid.NodeIdx array)
=
let alignedEdge =
match acc with
| a when a |> Array.contains (snd edge) -> (snd edge, fst edge)
| _ -> edge
let polygon = [| fst alignedEdge; snd alignedEdge |] |> Array.append acc
edges
|> Array.tryFind (fun (a, b) ->
(a = snd alignedEdge || b = snd alignedEdge)
&& (polygon |> Array.contains a |> not
|| polygon |> Array.contains b |> not)
)
|> function
| Some n ->
// Search deep, then check for remaining nodes
findPolygon n (acc |> Set.add node) |> findPolygon n
| None -> acc
| Some nextEdge -> findPolygon nextEdge polygon
| None -> polygon
// Find all polygons given array of edge nodes
let rec getGroups (remaining: Set<Grid.NodeIdx>) (acc: Grid.NodeIdx[][]) =
match remaining.Count with
/// Given an array of edges, group all edge nodes into polygons
let rec findAllPolygons
(remaining: (Grid.NodeIdx * Grid.NodeIdx) array)
(acc: Grid.NodeIdx array array)
=
match remaining.Length with
| 0 -> acc
| _ ->
let polygon = findPolygon (remaining |> Set.minElement) Set.empty
let newRemaining = Set.difference remaining polygon
getGroups
let polygon = findPolygon remaining[0] [||]
let newRemaining =
remaining
|> Array.filter (fun (a, b) ->
(polygon |> Array.contains a |> not)
&& (polygon |> Array.contains b |> not)
)
findAllPolygons
newRemaining
(Array.append acc [| polygon |> Array.ofSeq |])
(acc |> Array.append [| polygon |> Array.ofSeq |])
getGroups (nodes |> Set.ofArray) [||]
findAllPolygons edges [||]
/// <summary>
/// Create an array of strings to represent a polygon to input into Radovan's Delaunay program.
@@ -95,8 +103,18 @@ let groupByPolygon (nodes: Grid.NodeIdx[]) (nIdx: Grid.NeighborIndex) =
/// x and y coordinates
/// </summary>
let getPolygonString (nodes: Grid.NodeIdx array) (grid: Fvcom.FvcomGrid) =
nodes
|> Array.map (fun nIdx -> grid.Nodes[nIdx] |> fun n -> $"{fst n} {snd n}")
|> Array.append [| $"{nodes.Length}" |]
/// <summary>
/// Create an array of strings to represent polygons to input into the Makai tool
/// for visualization.
/// </summary>
let getPolygonStringMakai (nodes: Grid.NodeIdx array) (grid: Fvcom.FvcomGrid) =
nodes
|> Array.map (fun nIdx ->
grid.Nodes[nIdx] |> fun n -> $"{fst n} {snd n} 0.25"
)
|> Array.append [| $"{nodes.Length}" |]
// NOTE: Maybe we need to write boundary instead of island when applicable
|> Array.append [| $"island {nodes.Length} 0" |]