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

feat: organize into lib and cli

Browse Source
This commit is contained in:
Ole Tytlandsvik
2024-08-13 11:10:47 +02:00
parent ea56af6edc
commit 6200c17249
11 changed files with 388 additions and 254 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
.envrc Normal file
View File

@@ -0,0 +1 @@
use flake

2
.gitignore vendored
View File

@@ -6,4 +6,6 @@ obj/
.vscode
*.nc
*.txt
.direnv/*
*.sln.DotSettings.user
delaunay

7
Gridding.sln
View File

@@ -5,6 +5,8 @@ Project("{2150E333-8FDC-42A3-9474-1A3956D46DE8}") = "src", "src", "{A42886E3-363
EndProject
Project("{F2A71F9B-5D33-465A-A702-920D77279786}") = "CLI", "src\CLI\CLI.fsproj", "{80DF5F74-4690-4D8C-AA45-9E8758705192}"
EndProject
Project("{F2A71F9B-5D33-465A-A702-920D77279786}") = "Gridding", "src\Gridding\Gridding.fsproj", "{05D9B830-0CD6-4E59-9AD3-A196B82D60A8}"
EndProject
Global
GlobalSection(SolutionConfigurationPlatforms) = preSolution
Debug|Any CPU = Debug|Any CPU
@@ -15,8 +17,13 @@ Global
{80DF5F74-4690-4D8C-AA45-9E8758705192}.Debug|Any CPU.Build.0 = Debug|Any CPU
{80DF5F74-4690-4D8C-AA45-9E8758705192}.Release|Any CPU.ActiveCfg = Release|Any CPU
{80DF5F74-4690-4D8C-AA45-9E8758705192}.Release|Any CPU.Build.0 = Release|Any CPU
{05D9B830-0CD6-4E59-9AD3-A196B82D60A8}.Debug|Any CPU.ActiveCfg = Debug|Any CPU
{05D9B830-0CD6-4E59-9AD3-A196B82D60A8}.Debug|Any CPU.Build.0 = Debug|Any CPU
{05D9B830-0CD6-4E59-9AD3-A196B82D60A8}.Release|Any CPU.ActiveCfg = Release|Any CPU
{05D9B830-0CD6-4E59-9AD3-A196B82D60A8}.Release|Any CPU.Build.0 = Release|Any CPU
EndGlobalSection
GlobalSection(NestedProjects) = preSolution
{80DF5F74-4690-4D8C-AA45-9E8758705192} = {A42886E3-3635-45C7-AFB9-3271A447B3A0}
{05D9B830-0CD6-4E59-9AD3-A196B82D60A8} = {A42886E3-3635-45C7-AFB9-3271A447B3A0}
EndGlobalSection
EndGlobal

61
flake.lock generated Normal file
View File

@@ -0,0 +1,61 @@
{
"nodes": {
"flake-utils": {
"inputs": {
"systems": "systems"
},
"locked": {
"lastModified": 1710146030,
"narHash": "sha256-SZ5L6eA7HJ/nmkzGG7/ISclqe6oZdOZTNoesiInkXPQ=",
"owner": "numtide",
"repo": "flake-utils",
"rev": "b1d9ab70662946ef0850d488da1c9019f3a9752a",
"type": "github"
},
"original": {
"owner": "numtide",
"repo": "flake-utils",
"type": "github"
}
},
"nixpkgs": {
"locked": {
"lastModified": 1723362943,
"narHash": "sha256-dFZRVSgmJkyM0bkPpaYRtG/kRMRTorUIDj8BxoOt1T4=",
"owner": "NixOS",
"repo": "nixpkgs",
"rev": "a58bc8ad779655e790115244571758e8de055e3d",
"type": "github"
},
"original": {
"owner": "NixOS",
"ref": "nixos-unstable",
"repo": "nixpkgs",
"type": "github"
}
},
"root": {
"inputs": {
"flake-utils": "flake-utils",
"nixpkgs": "nixpkgs"
}
},
"systems": {
"locked": {
"lastModified": 1681028828,
"narHash": "sha256-Vy1rq5AaRuLzOxct8nz4T6wlgyUR7zLU309k9mBC768=",
"owner": "nix-systems",
"repo": "default",
"rev": "da67096a3b9bf56a91d16901293e51ba5b49a27e",
"type": "github"
},
"original": {
"owner": "nix-systems",
"repo": "default",
"type": "github"
}
}
},
"root": "root",
"version": 7
}

32
flake.nix Normal file
View File

@@ -0,0 +1,32 @@
{
inputs = {
nixpkgs.url = "github:NixOS/nixpkgs/nixos-unstable";
flake-utils = {
url = "github:numtide/flake-utils";
inputs.nixpkgs.follows = "nixpkgs";
};
};
outputs =
{
self,
nixpkgs,
flake-utils,
...
}:
flake-utils.lib.eachDefaultSystem (
system:
let
pkgs = import nixpkgs { inherit system; };
in
with pkgs;
{
devShells.default = mkShell {
nativeBuildInputs = [
netcdf
dotnet-sdk_8
];
LD_LIBRARY_PATH = lib.makeLibraryPath [netcdf];
};
}
);
}

4
src/CLI/CLI.fsproj
View File

@@ -13,5 +13,9 @@
<PackageReference Include="Oceanbox.FvcomKit" Version="5.5.3" />
<PackageReference Include="SDSLite-O" Version="2.7.2" />
</ItemGroup>
<ItemGroup>
<ProjectReference Include="..\Gridding\Gridding.fsproj" />
</ItemGroup>
</Project>

259
src/CLI/Program.fs
View File

@@ -1,260 +1,11 @@
open System
open System.IO
open System.IO
open Oceanbox.FvcomKit
open Microsoft.Research.Science.Data
open ProjNet.FSharp
open Gridding.Polygon
[<Literal>]
let filePath = "testfile.nc"
let private projection = CoordSys.UTMn 33
let private utm2ll =
makeTransform projection CoordSys.WGS84
/// <summary>
/// Separate coastline nodes from inner nodes.
/// Returns (coastline nodes * inner nodes).
/// </summary>
/// <param name="nIdx"></param>
/// <param name="grid"></param>
let separateCoastline (nIdx: Grid.NeighborIndex) (grid: Grid.Grid) =
let getEdges (a, b, c) =
[(a, b); (b, a); (a, c); (c, a); (b, c); (c, b)]
|> Set.ofList
// Find whether an elem has less than two elems that it
// shares an edge with
let isEdgeElem (elemIdx: Grid.ElemIdx): bool =
let ownEdges = grid.Elem[elemIdx] |> getEdges
Grid.getElemsSurroundingElem nIdx grid elemIdx
|> Array.filter (fun surroundingIdx ->
// Avoid including self
if surroundingIdx = elemIdx then
false
else
grid.Elem[surroundingIdx]
|> getEdges
|> Set.intersect ownEdges
|> Set.count > 0
)
|> Array.length <= 2
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)
|> Set.ofArray
let iNodes =
[| 0 .. (grid.Nodes.Length - 1) |]
|> Set.ofArray
|> fun allNodes -> Set.difference allNodes cNodes
(cNodes |> Array.ofSeq, iNodes |> Array.ofSeq)
/// <summary>
/// Group given edge nodes by polygon
/// </summary>
/// <param name="nodes">Edge node indices</param>
/// <param name="nIdx"></param>
let groupByPolygon (nodes: Grid.NodeIdx []) (nIdx: Grid.NeighborIndex) =
let nMap = nIdx.NodesAroundNode
let isEdgeNode (node: Grid.NodeIdx) =
nodes
|> Array.contains node
// Find a single polygon given an edge node
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)))
|> function
| Some n ->
// Search deep, then check for remaining nodes
findPolygon n (acc |> Set.add node)
|> findPolygon n
| None -> acc
// Find all polygons given array of edge nodes
let rec getGroups (remaining: Set<Grid.NodeIdx>) (acc: Grid.NodeIdx[][]) =
match remaining.Count with
| 0 -> acc
| _ ->
let polygon =
findPolygon (remaining |> Set.minElement) Set.empty
let newRemaining =
Set.difference remaining polygon
getGroups
(newRemaining)
(Array.append acc [| polygon |> Array.ofSeq |])
getGroups (nodes |> Set.ofArray) [||]
/// <summary>
/// Return a subset of the given node indices
/// </summary>
/// <param name="nodes"></param>
let reduceNodes (nodes: Grid.NodeIdx []) =
nodes
|> Array.choose (fun x ->
// Randomly remove half of array
let rnd = Random ()
rnd.Next 2
|> function
| 0 -> Some x
| _ -> None)
/// <summary>
/// Write variables in Fvcom grid type to file. These include
/// x, y, lon, lat, h
/// </summary>
/// <param name="nodes">Node indices from grid to write</param>
/// <param name="grid">Original grid</param>
/// <param name="file">Open NetCDFDataSet file handle</param>
let writeFGrid (nodes: Grid.NodeIdx []) (grid: Fvcom.FvcomGrid) (file: NetCDF4.NetCDFDataSet) =
let x, y =
nodes
|> Array.map (fun idx ->
grid.Nodes.[idx]
|> (fun n -> float32 (fst n), float32 (snd n)))
|> Array.unzip
let lon, lat = Array.zip x y |> Array.map utm2ll.project |> Array.unzip
let h =
nodes
|> Array.map (fun idx ->
grid.Bathymetry.[idx])
printfn $"node dimension length: {nodes.Length}"
printfn $"x length: {x.Length}, y Length: {y.Length}"
file.AddVariable<float32>("x", x, [|"node"|]) |> ignore
file.AddVariable<float32>("y", y, [|"node"|]) |> ignore
file.AddVariable<float32>("lon", lon, [|"node"|]) |> ignore
file.AddVariable<float32>("lat", lat, [|"node"|]) |> ignore
file.AddVariable<float32>("h", h, [|"node"|]) |> ignore
file
/// <summary>
/// Return an array containing only the items on the
/// given indices
/// </summary>
/// <param name="nodes">Indices to keep</param>
/// <param name="arr">Array to reduce</param>
let private reduce (nodes: Grid.NodeIdx []) (arr: 'T []) =
nodes |> Array.map (fun idx -> arr[idx])
/// <summary>
/// Write variables into given ncdf file.
/// Does not write the following variables:
/// x, y, lon, lat, h
/// </summary>
/// <param name="nodes"></param>
/// <param name="ds"></param>
/// <param name="file"></param>
let writeVariables (nodes: Grid.NodeIdx []) (ds: DataSet) (file: NetCDF4.NetCDFDataSet) =
let timeFrames = Fvcom.getNumFrames ds
file.CreateDimension("time", timeFrames)
let siglay = Fvcom.getNumSiglay ds
// Don't create dimension, use only z dimension when writing new file
// Read zeta values for each time frame,
// add readings to a 2d array
let rec readZetas (t: int) (acc: single [][]) =
match t with
| x when x >= timeFrames -> acc
| _ ->
Fvcom.readZeta ds t
|> Array.map single
|> reduce nodes
|> fun a -> [|a|]
|> Array.append acc
|> readZetas (t + 1)
let reducedZetas =
readZetas 0 [||]
// |> fun a ->
// printfn $"{a.Length}"
// a
|> fun a -> Array2D.init timeFrames nodes.Length (fun i j -> a[i][j])
// |> fun a ->
// printfn $"{a.Length}"
// a
file.AddVariable<single>("zeta", Array.empty, [|"time"; "node"|]) |> ignore
file.Variables[ "zeta" ].PutData([| 0; 0 |], reducedZetas)
// Read salinity for each timeframe,
// use top and bottom layers only
let rec readSalinity (t: int) (acc: single [][][]) =
// Helper to read at given layer
let readAtLayer (l: int) =
Fvcom.readSalinity ds t l
|> Array.map single
|> reduce nodes
|> fun a -> [|a|]
match t with
| x when x >= timeFrames -> acc
| _ ->
let bottomLayer = readAtLayer 0
let topLayer = readAtLayer siglay
Array.append bottomLayer topLayer
|> fun a -> [|a|]
|> Array.append acc
|> readSalinity (t + 1)
let reducedSalinity =
readSalinity 0 [||]
|> fun a -> Array3D.init timeFrames 2 nodes.Length (fun i j k -> a.[i].[j].[k])
file.AddVariable<single>("salinity", Array.empty, [|"time"; "z"; "node"|]) |> ignore
file.Variables[ "salinity" ].PutData([| 0; 0; 0 |], reducedSalinity)
/// <summary>
/// Write a new grid based on an original grid using only data from the
/// node indices given
/// </summary>
/// <param name="nodes">Remaining nodes after reduction</param>
/// <param name="grid">Grid structure read from the original grid</param>
/// <param name="ds">The original dataset</param>
/// <param name="filename">Path of output file</param>
let writeReducedGrid (nodes: Grid.NodeIdx []) (grid: Fvcom.FvcomGrid) (ds: DataSet) (filename: string) =
let file = new NetCDF4.NetCDFDataSet(filename)
file.CreateDimension("node", nodes.Length)
file.CreateDimension("z", 2)
file
|> writeFGrid nodes grid
|> writeVariables nodes ds
file.Commit ()
file.Dispose ()
/// <summary>
/// Create an array of strings to represent a polygon to input into Radovan's Delaunay program.
/// First prints the length of the number of nodes in the polygon before printing each node's
/// 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}" |]
[<EntryPoint>]
let main args =
@@ -304,11 +55,11 @@ let main args =
File.WriteAllLines("boundary.txt", boundaryStrings)
File.WriteAllLines("islands.txt", restStrings)
let reducedInnerNodes =
innerNodes |> reduceNodes
// let reducedInnerNodes =
// innerNodes |> reduceNodes
let innerNodeStrings =
reducedInnerNodes
innerNodes
|> Array.map (fun nIdx ->
grd.Nodes[nIdx]
|> fun n -> $"{fst n} {snd n}")

139
src/Gridding/Fvcom.fs Normal file
View File

@@ -0,0 +1,139 @@
module Gridding.Fvcom
open Oceanbox.FvcomKit
open Microsoft.Research.Science.Data
open ProjNet.FSharp
let private projection = CoordSys.UTMn 33
let private utm2ll =
makeTransform projection CoordSys.WGS84
/// <summary>
/// Return an array containing only the items on the
/// given indices
/// </summary>
/// <param name="nodes">Indices to keep</param>
/// <param name="arr">Array to reduce</param>
let private reduce (nodes: Grid.NodeIdx []) (arr: 'T []) =
nodes |> Array.map (fun idx -> arr[idx])
/// <summary>
/// Write variables in Fvcom grid type to file. These include
/// x, y, lon, lat, h
/// </summary>
/// <param name="nodes">Node indices from grid to write</param>
/// <param name="grid">Original grid</param>
/// <param name="file">Open NetCDFDataSet file handle</param>
let private writeFGrid (nodes: Grid.NodeIdx []) (grid: Fvcom.FvcomGrid) (file: NetCDF4.NetCDFDataSet) =
let x, y =
nodes
|> Array.map (fun idx ->
grid.Nodes[idx]
|> (fun n -> float32 (fst n), float32 (snd n)))
|> Array.unzip
let lon, lat = Array.zip x y |> Array.map utm2ll.project |> Array.unzip
let h =
nodes
|> Array.map (fun idx ->
grid.Bathymetry.[idx])
printfn $"node dimension length: {nodes.Length}"
printfn $"x length: {x.Length}, y Length: {y.Length}"
file.AddVariable<float32>("x", x, [|"node"|]) |> ignore
file.AddVariable<float32>("y", y, [|"node"|]) |> ignore
file.AddVariable<float32>("lon", lon, [|"node"|]) |> ignore
file.AddVariable<float32>("lat", lat, [|"node"|]) |> ignore
file.AddVariable<float32>("h", h, [|"node"|]) |> ignore
file
/// <summary>
/// Write variables into given ncdf file.
/// Does not write the following variables:
/// x, y, lon, lat, h
/// </summary>
/// <param name="nodes"></param>
/// <param name="ds"></param>
/// <param name="file"></param>
let private writeVariables (nodes: Grid.NodeIdx []) (ds: DataSet) (file: NetCDF4.NetCDFDataSet) =
let timeFrames = Fvcom.getNumFrames ds
file.CreateDimension("time", timeFrames)
let siglay = Fvcom.getNumSiglay ds
// Don't create dimension, use only z dimension when writing new file
// Read zeta values for each time frame,
// add readings to a 2d array
let rec readZetas (t: int) (acc: single [][]) =
match t with
| x when x >= timeFrames -> acc
| _ ->
Fvcom.readZeta ds t
|> Array.map single
|> reduce nodes
|> fun a -> [|a|]
|> Array.append acc
|> readZetas (t + 1)
let reducedZetas =
readZetas 0 [||]
// |> fun a ->
// printfn $"{a.Length}"
// a
|> fun a -> Array2D.init timeFrames nodes.Length (fun i j -> a[i][j])
// |> fun a ->
// printfn $"{a.Length}"
// a
file.AddVariable<single>("zeta", Array.empty, [|"time"; "node"|]) |> ignore
file.Variables[ "zeta" ].PutData([| 0; 0 |], reducedZetas)
// Read salinity for each timeframe,
// use top and bottom layers only
let rec readSalinity (t: int) (acc: single [][][]) =
// Helper to read at given layer
let readAtLayer (l: int) =
Fvcom.readSalinity ds t l
|> Array.map single
|> reduce nodes
|> fun a -> [|a|]
match t with
| x when x >= timeFrames -> acc
| _ ->
let bottomLayer = readAtLayer 0
let topLayer = readAtLayer siglay
Array.append bottomLayer topLayer
|> fun a -> [|a|]
|> Array.append acc
|> readSalinity (t + 1)
let reducedSalinity =
readSalinity 0 [||]
|> fun a -> Array3D.init timeFrames 2 nodes.Length (fun i j k -> a.[i].[j].[k])
file.AddVariable<single>("salinity", Array.empty, [|"time"; "z"; "node"|]) |> ignore
file.Variables[ "salinity" ].PutData([| 0; 0; 0 |], reducedSalinity)
/// <summary>
/// Write a new grid based on an original grid using only data from the
/// node indices given
/// </summary>
/// <param name="nodes">Remaining nodes after reduction</param>
/// <param name="grid">Grid structure read from the original grid</param>
/// <param name="ds">The original dataset</param>
/// <param name="filename">Path of output file</param>
let writeReducedGrid (nodes: Grid.NodeIdx []) (grid: Fvcom.FvcomGrid) (ds: DataSet) (filename: string) =
let file = new NetCDF4.NetCDFDataSet(filename)
file.CreateDimension("node", nodes.Length)
file.CreateDimension("z", 2)
file
|> writeFGrid nodes grid
|> writeVariables nodes ds
file.Commit ()
file.Dispose ()

20
src/Gridding/Gridding.fsproj Normal file
View File

@@ -0,0 +1,20 @@
<Project Sdk="Microsoft.NET.Sdk">
<PropertyGroup>
<OutputType>Library</OutputType>
<TargetFramework>net8.0</TargetFramework>
<IsPackable>true</IsPackable>
</PropertyGroup>
<ItemGroup>
<Compile Include="Polygon.fs" />
<Compile Include="Fvcom.fs" />
<Compile Include="Reduce.fs" />
</ItemGroup>
<ItemGroup>
<PackageReference Include="Oceanbox.FvcomKit" Version="5.5.3" />
<PackageReference Include="SDSLite-O" Version="2.7.2" />
</ItemGroup>
</Project>

99
src/Gridding/Polygon.fs Normal file
View File

@@ -0,0 +1,99 @@
module Gridding.Polygon
open Oceanbox.FvcomKit
/// <summary>
/// Separate coastline nodes from inner nodes.
/// Returns (coastline nodes * inner nodes).
/// </summary>
/// <param name="nIdx"></param>
/// <param name="grid"></param>
let separateCoastline (nIdx: Grid.NeighborIndex) (grid: Grid.Grid) =
let getEdges (a, b, c) =
[(a, b); (b, a); (a, c); (c, a); (b, c); (c, b)]
|> Set.ofList
// Find whether an elem has less than two elems that it
// shares an edge with
let isEdgeElem (elemIdx: Grid.ElemIdx): bool =
let ownEdges = grid.Elem[elemIdx] |> getEdges
Grid.getElemsSurroundingElem nIdx grid elemIdx
|> Array.filter (fun surroundingIdx ->
// Avoid including self
if surroundingIdx = elemIdx then
false
else
grid.Elem[surroundingIdx]
|> getEdges
|> Set.intersect ownEdges
|> Set.count > 0
)
|> Array.length <= 2
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)
|> Set.ofArray
let iNodes =
[| 0 .. (grid.Nodes.Length - 1) |]
|> Set.ofArray
|> fun allNodes -> Set.difference allNodes cNodes
(cNodes |> Array.ofSeq, iNodes |> Array.ofSeq)
/// <summary>
/// Group given edge nodes by polygon
/// </summary>
/// <param name="nodes">Edge node indices</param>
/// <param name="nIdx"></param>
let groupByPolygon (nodes: Grid.NodeIdx []) (nIdx: Grid.NeighborIndex) =
let nMap = nIdx.NodesAroundNode
let isEdgeNode (node: Grid.NodeIdx) =
nodes
|> Array.contains node
// Find a single polygon given an edge node
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)))
|> function
| Some n ->
// Search deep, then check for remaining nodes
findPolygon n (acc |> Set.add node)
|> findPolygon n
| None -> acc
// Find all polygons given array of edge nodes
let rec getGroups (remaining: Set<Grid.NodeIdx>) (acc: Grid.NodeIdx[][]) =
match remaining.Count with
| 0 -> acc
| _ ->
let polygon =
findPolygon (remaining |> Set.minElement) Set.empty
let newRemaining =
Set.difference remaining polygon
getGroups
(newRemaining)
(Array.append acc [| polygon |> Array.ofSeq |])
getGroups (nodes |> Set.ofArray) [||]
/// <summary>
/// Create an array of strings to represent a polygon to input into Radovan's Delaunay program.
/// First prints the length of the number of nodes in the polygon before printing each node's
/// 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}" |]

18
src/Gridding/Reduce.fs Normal file
View File

@@ -0,0 +1,18 @@
module Gridding.Reduce
open System
open Oceanbox.FvcomKit
/// <summary>
/// Return a subset of the given node indices
/// </summary>
/// <param name="nodes"></param>
let reduceNodes (nodes: Grid.NodeIdx []) =
nodes
|> Array.choose (fun x ->
// Randomly remove half of array
let rnd = Random ()
rnd.Next 2
|> function
| 0 -> Some x
| _ -> None)