oceanbox/Oceanbox.FvcomKit
6
0
Fork 0
Code Issues Pull Requests 1 Actions Packages Projects Releases 80 Wiki Activity

devel: format with fantomas

Browse Source
This commit is contained in:
Jonas Juselius
2025-03-06 07:33:13 +01:00
parent bf7ae5889b
commit 222b022662
13 changed files with 408 additions and 571 deletions
Download Patch File Download Diff File Expand all files Collapse all files

122
src/Adjoin.fs
View File

@@ -9,10 +9,10 @@ open Polygon
open KdTree // C# version
type CropBox = (float * float) * (float * float)
type Mask = bool [,]
type Mask = bool[,]
type Pos = float * float
type PosVec = Pos []
type BiPos = float [,] * float [,]
type PosVec = Pos[]
type BiPos = float[,] * float[,]
let private reProject (proj: IProj) ((lng, lat): BiPos) =
lng
@@ -20,7 +20,7 @@ let private reProject (proj: IProj) ((lng, lat): BiPos) =
|> Array2D.map proj.project
|> unzip2D
let private genCropIdx ((min, max): float * float) (m: float [,]) =
let private genCropIdx ((min, max): float * float) (m: float[,]) =
[|
for i = 0 to Array2D.length1 m - 1 do
for j = 0 to Array2D.length2 m - 1 do
@@ -44,22 +44,21 @@ let private mkCropBox (box: BBox) : CropBox =
// make an array of indeces of wet roms grid points inside the fvcom domain
let private genCullIdx (box: CropBox) (wet: Mask) (p: BiPos) =
genCropMask box p
|> Array.filter (fun (i, j) -> wet[i, j])
genCropMask box p |> Array.filter (fun (i, j) -> wet[i, j])
// crop roms grid coordinates based in index mask of active, overlapping points
let private cullCoords (p: BiPos) (idx: (int * int) []) =
let private cullCoords (p: BiPos) (idx: (int * int)[]) =
let x, y = p
let x' = idx |> Array.map (fun (i, j) -> x[i, j])
let y' = idx |> Array.map (fun (i, j) -> y[i, j])
Array.zip x' y'
let private createTree (points: Leaf<int> []) =
let tree = KdTree<float, int>(2, KdTree.Math.DoubleMath())
let private createTree (points: Leaf<int>[]) =
let tree = KdTree<float, int> (2, KdTree.Math.DoubleMath ())
points
|> Array.iter (fun a -> tree.Add([| fst a.Pos; snd a.Pos |], a.Data) |> ignore)
|> Array.iter (fun a -> tree.Add ([| fst a.Pos; snd a.Pos |], a.Data) |> ignore)
if points.Length > 0 then
tree.Balance()
tree.Balance ()
else
// Log.Warning $"Empty kd-tree"
()
@@ -75,7 +74,7 @@ let private buildTree (points: (float * float) array) =
// adjoinIdx: index to culled roms grid for each fvcom node
// cullIdx: index to roms points _actually_ in use
// oobIdx: out-of-bounds points, too far from a roms cell
type FvcomAdjoint = { adjoinIdx: int []; cullIdx: (int * int) []; oobIdx: int [] }
type FvcomAdjoint = { adjoinIdx: int[]; cullIdx: (int * int)[]; oobIdx: int[] }
let private dist (a: float * float) (b: float * float) =
let x0, y0 = a
@@ -107,14 +106,13 @@ let private distLngLat (proj: IProj) (a: float * float) (b: float * float) =
// (0, [])
// |> snd
// |> Array.ofList
let private genOobIdx (proj: IProj) (tree: KdTree<_, _>) (cullIdx: (int * int) []) (fPos: PosVec) (rPos: BiPos) =
let private genOobIdx (proj: IProj) (tree: KdTree<_, _>) (cullIdx: (int * int)[]) (fPos: PosVec) (rPos: BiPos) =
let rLon, rLat = rPos
let dMax =
distLngLat proj (rLon[0, 0], rLat[0, 0]) (rLon[0, 1], rLat[0, 1])
let dMax = distLngLat proj (rLon[0, 0], rLat[0, 0]) (rLon[0, 1], rLat[0, 1])
fPos
|> Array.fold
(fun (n, a) (x, y as p0) ->
tree.GetNearestNeighbours([| x; y |], 1)
tree.GetNearestNeighbours ([| x; y |], 1)
|> Array.head
|> fun p ->
let i0, i1 = cullIdx[p.Value]
@@ -148,15 +146,8 @@ let private mkFvcomAdjoint (proj: IProj) (fPos: PosVec) (bbox: BBox) ((rPos, wet
let tree = cullCoords pos' cullIdx |> buildTree
let nearest =
fPos
|> Array.map (fun (x, y) ->
tree.GetNearestNeighbours([| x; y |], 1)
|> Array.head
|> fun x -> x.Value)
{
adjoinIdx = nearest
cullIdx = cullIdx
oobIdx = genOobIdx proj tree cullIdx fPos rPos
}
|> Array.map (fun (x, y) -> tree.GetNearestNeighbours ([| x; y |], 1) |> Array.head |> (fun x -> x.Value))
{ adjoinIdx = nearest; cullIdx = cullIdx; oobIdx = genOobIdx proj tree cullIdx fPos rPos }
let inline float2 x = bimap float float x
@@ -169,19 +160,19 @@ let getNearestCell (proj: IProj) (fvcom: Grid) (roms: BiPos * Mask) =
let cells = Util.calcCentroids fvcom |> Array.map float2
mkFvcomAdjoint proj cells fvcom.BBox roms
let private createIdxTree (points: Leaf<int * int> []) =
let tree = KdTree<float, int * int>(2, KdTree.Math.DoubleMath())
let private createIdxTree (points: Leaf<int * int>[]) =
let tree = KdTree<float, int * int> (2, KdTree.Math.DoubleMath ())
points
|> Array.iter (fun a -> tree.Add([| fst a.Pos; snd a.Pos |], a.Data) |> ignore)
|> Array.iter (fun a -> tree.Add ([| fst a.Pos; snd a.Pos |], a.Data) |> ignore)
if points.Length > 0 then
tree.Balance()
tree.Balance ()
else
//Log.Warning $"Empty kd-tree"
()
tree
// rPos must be in the right projection!
let makeNestTree ((lng, lat): float [,] * float [,]) =
let makeNestTree ((lng, lat): float[,] * float[,]) =
[| 0 .. (Array2D.length1 lng) - 1 |]
|> Array.map (fun i ->
[| 0 .. (Array2D.length2 lng) - 1 |]
@@ -196,36 +187,31 @@ let makeNestTree ((lng, lat): float [,] * float [,]) =
let getNearestUpperLeft' (tree: KdTree<float, int * int>) (((ew, ns), mask): BiPos * Mask) ((x, y): float * float) =
let i, j =
// nearestNeighbor tree { X = x; Y = y }
tree.GetNearestNeighbours([| x; y |], 1)
|> fun x -> x[0].Value
let p1 =
[|
ew[i, j], ns[i, j]
ew[i + 1, j], ns[i + 1, j]
ew[i + 1, j - 1], ns[i + 1, j - 1]
ew[i, j - 1], ns[i, j - 1]
|] // lower right grid cell
let p2 =
[|
ew[i, j], ns[i, j]
ew[i + 1, j], ns[i + 1, j]
ew[i + 1, j + 1], ns[i + 1, j + 1]
ew[i, j + 1], ns[i, j + 1]
|] // upper right grid cell
let p3 =
[|
ew[i, j], ns[i, j]
ew[i - 1, j], ns[i - 1, j]
ew[i - 1, j - 1], ns[i - 1, j - 1]
ew[i, j - 1], ns[i, j - 1]
|] // lower left grid cell
let p4 =
[|
ew[i, j], ns[i, j]
ew[i - 1, j], ns[i - 1, j]
ew[i - 1, j + 1], ns[i - 1, j + 1]
ew[i, j + 1], ns[i, j + 1]
|] // upper left grid cell
tree.GetNearestNeighbours ([| x; y |], 1) |> fun x -> x[0].Value
let p1 = [|
ew[i, j], ns[i, j]
ew[i + 1, j], ns[i + 1, j]
ew[i + 1, j - 1], ns[i + 1, j - 1]
ew[i, j - 1], ns[i, j - 1]
|] // lower right grid cell
let p2 = [|
ew[i, j], ns[i, j]
ew[i + 1, j], ns[i + 1, j]
ew[i + 1, j + 1], ns[i + 1, j + 1]
ew[i, j + 1], ns[i, j + 1]
|] // upper right grid cell
let p3 = [|
ew[i, j], ns[i, j]
ew[i - 1, j], ns[i - 1, j]
ew[i - 1, j - 1], ns[i - 1, j - 1]
ew[i, j - 1], ns[i, j - 1]
|] // lower left grid cell
let p4 = [|
ew[i, j], ns[i, j]
ew[i - 1, j], ns[i - 1, j]
ew[i - 1, j + 1], ns[i - 1, j + 1]
ew[i, j + 1], ns[i, j + 1]
|] // upper left grid cell
let q =
if inpolygon p1 (x, y) then
[| i + 1, j - 1; i + 1, j; i, j; i, j - 1 |]
@@ -237,9 +223,7 @@ let getNearestUpperLeft' (tree: KdTree<float, int * int>) (((ew, ns), mask): BiP
[| i, j; i, j + 1; i - 1, j + 1; i - 1, j |]
else
[| 0, 0; 0, 0; 0, 0; 0, 0 |]
if (q
|> Array.map (fun (q1, q2) -> q1 + q2)
|> Array.sum) = 0 then
if (q |> Array.map (fun (q1, q2) -> q1 + q2) |> Array.sum) = 0 then
failwith "Surrounding grid cell not found"
let m = q |> Array.map (fun (n, m) -> mask[n, m])
Array.zip q m
@@ -299,24 +283,24 @@ let getNearestUpperLeft' (tree: KdTree<float, int * int>) (((ew, ns), mask): BiP
// let m = q |> Array.map (fun (n, m) -> mask[n, m])
// Array.zip q m
let getNearestCellCorner (coords, _ as grid: BiPos * Mask) (pos: (float * float) []) =
let getNearestCellCorner (coords, _ as grid: BiPos * Mask) (pos: (float * float)[]) =
let tree = makeNestTree coords
// pos |> Array.map (getNearestUpperLeft tree grid)
pos |> Array.map (getNearestUpperLeft' tree grid)
let getCellBox ((ew, ns): BiPos) (boxid: ((int * int) * bool) []) =
let getCellBox ((ew, ns): BiPos) (boxid: ((int * int) * bool)[]) =
boxid
|> Array.map (fun ((n, m), mask) -> (ew[n, m], ns[n, m]), mask)
|> Array.unzip
let getCellProps (prop: float [,]) ((n, m): int * int) =
let getCellProps (prop: float[,]) ((n, m): int * int) =
prop[n, m], prop[n + 1, m], prop[n + 1, m + 1], prop[n + 1, m]
// pick out elements actually in use
let inline private cullBiMatrix (culler: (int * int) []) (m: 'a [,]) =
let inline private cullBiMatrix (culler: (int * int)[]) (m: 'a[,]) =
culler |> Array.map (fun (i, j) -> m[i, j])
// adjoin fvcom and roms data based on nearest neighbours
let adjoinMatrix (adj: FvcomAdjoint) (m: 'a [,]) =
let adjoinMatrix (adj: FvcomAdjoint) (m: 'a[,]) =
let x = cullBiMatrix adj.cullIdx m
adj.adjoinIdx |> Array.map (fun n -> x[n])
adj.adjoinIdx |> Array.map (fun n -> x[n])

40
src/Arome.fs
View File

@@ -15,8 +15,7 @@ type SquareGrid = {
BBox: BBox
squareSize: float
projection: ProjNet.CoordinateSystems.CoordinateSystem
}
with
} with
member this.getBoundingBox() = this.BBox
static member empty = {
dimensions = 0, 0
@@ -40,16 +39,15 @@ let getBBox (xs: float array) (ys: float array) : BBox =
maxY = maxY
center = center
}
with
| e ->
with e ->
Log.Error $"{e}"
BBox.empty
let getGrid (ds: DataSet) : Result<SquareGrid, string> =
try
let dimensions = ds.Dimensions["x"].Length, ds.Dimensions["y"].Length
let xs = (ds["x"].GetData() :?> single []) |> Array.map float
let ys = (ds["y"].GetData() :?> single []) |> Array.map float
let xs = (ds["x"].GetData () :?> single[]) |> Array.map float
let ys = (ds["y"].GetData () :?> single[]) |> Array.map float
let bbox = getBBox xs ys
if xs.Length < 2 && ys.Length < 2 then
@@ -62,22 +60,26 @@ let getGrid (ds: DataSet) : Result<SquareGrid, string> =
let isSquare = lengthX = lengthY
if isSquare then
{ SquareGrid.empty with
dimensions = dimensions
BBox = bbox
squareSize = lengthX
}
{ SquareGrid.empty with dimensions = dimensions; BBox = bbox; squareSize = lengthX }
|> Ok
else
Log.Error("FvcomKit.Arome.getGrid grid is not square: {X1} - {X0} = {LengthX} = {LengthY} = {Y1} - {Y0}", x1, x0, lengthX, lengthY, y1, y0)
Log.Error (
"FvcomKit.Arome.getGrid grid is not square: {X1} - {X0} = {LengthX} = {LengthY} = {Y1} - {Y0}",
x1,
x0,
lengthX,
lengthY,
y1,
y0
)
Error "The given data set does not contain a grid made up of squares"
with exn ->
Log.Error(exn, "FvcomKit.Arome.getGrid exception")
Log.Error (exn, "FvcomKit.Arome.getGrid exception")
Error $"Error reading arome grid: {exn.Message}"
let readUV (ds: DataSet) t x y =
let u = ds[ "x_wind_10m" ].GetData([| t; y; x |], [| 1; 1; 1 |]) :?> single [,,]
let v = ds[ "y_wind_10m" ].GetData([| t; y; x |], [| 1; 1; 1 |]) :?> single [,,]
let u = ds["x_wind_10m"].GetData ([| t; y; x |], [| 1; 1; 1 |]) :?> single[,,]
let v = ds["y_wind_10m"].GetData ([| t; y; x |], [| 1; 1; 1 |]) :?> single[,,]
u[0, 0, 0], v[0, 0, 0]
/// Finds the index of a tile within a square grid, given its bounding box and square length
@@ -91,9 +93,13 @@ let tryFindIndex squareSize (wide, tall) (minX, minY) (maxX, maxY) (p0, p1) =
if xIdx < wide && yIdx < tall then
Some (xIdx, yIdx)
else
Log.Warning(
Log.Warning (
"Got wrong indices within the bounding box of the archive: min {@Min}, max {@Max}, point {@Point}, delta {@Delta}m, indices {@Indices}",
(minX, minY), (maxX, maxY), (p0, p1), (dx, dy), (xIdx, yIdx)
(minX, minY),
(maxX, maxY),
(p0, p1),
(dx, dy),
(xIdx, yIdx)
)
None
else

431
src/Fvcom.fs
View File

@@ -10,17 +10,16 @@ open Grid
open Types
type FvcomGrid =
{
Elem: Elem array
Nodes: Node array
BBox: BBox
Cells: Node array
Bathymetry: single []
Siglay: single [,]
SiglayCenter: single [,]
Siglev: single [,]
}
type FvcomGrid = {
Elem: Elem array
Nodes: Node array
BBox: BBox
Cells: Node array
Bathymetry: single[]
Siglay: single[,]
SiglayCenter: single[,]
Siglev: single[,]
} with
interface IGrid with
member x.getVertex n = x.Nodes[n]
member x.getCell n = x.Elem[n]
@@ -30,109 +29,90 @@ type FvcomGrid =
member x.getVertices() = x.Nodes
member x.getCells() = x.Elem
member x.getBoundingBox() = x.BBox
static member empty =
{
Elem = Array.empty
Nodes = Array.empty
BBox = BBox.empty
Cells = Array.empty
Bathymetry = Array.empty
Siglay = Array2D.zeroCreate 0 0
SiglayCenter = Array2D.zeroCreate 0 0
Siglev = Array2D.zeroCreate 0 0
}
member this.ToGrid() =
{
Elem = this.Elem
Nodes = this.Nodes
BBox = this.BBox
}
static member empty = {
Elem = Array.empty
Nodes = Array.empty
BBox = BBox.empty
Cells = Array.empty
Bathymetry = Array.empty
Siglay = Array2D.zeroCreate 0 0
SiglayCenter = Array2D.zeroCreate 0 0
Siglev = Array2D.zeroCreate 0 0
}
member this.ToGrid() = { Elem = this.Elem; Nodes = this.Nodes; BBox = this.BBox }
let getNumFrames (ds: DataSet) =
try
ds.Dimensions["time"].Length
with
| e ->
with e ->
Log.Error $"{e}"
0
let getNumSiglay (ds: DataSet) =
try
ds.Dimensions["siglay"].Length
with
| e ->
with e ->
Log.Error $"{e}"
0
let getNumSiglev (ds: DataSet) =
try
ds.Dimensions["siglev"].Length
with
| e ->
with e ->
Log.Error $"{e}"
0
let getTime (ds: DataSet) n =
try
let ts = ds[ "Times" ].GetData() :?> byte [,]
let ts = ds["Times"].GetData () :?> byte[,]
ts[n, *]
|> Array.map char
|> System.String
|> System.DateTime.Parse
|> Some
with
| e ->
ts[n, *] |> Array.map char |> System.String |> System.DateTime.Parse |> Some
with e ->
Log.Error $"getTime exception: {e.Message}"
None
let getTimeSpanSinceStart (ds: DataSet) n =
try
let days = ds[ "Itime" ].GetData() :?> int []
let msec = ds[ "Itime2" ].GetData() :?> int []
let days = ds["Itime"].GetData () :?> int[]
let msec = ds["Itime2"].GetData () :?> int[]
let t0 = TimeSpan.FromDays days[n]
let t1 = TimeSpan.FromMilliseconds msec[n]
t0 + t1 |> Some
with
| e ->
with e ->
Log.Error $"getTimeInDays exception: {e.Message}"
None
let readTauc (ds: DataSet) t =
try
let n = ds.Dimensions["nele"].Length
let tauc = ds[ "tauc" ].GetData([| t; 0 |], [| 1; n |]) :?> single [,]
let tauc = ds["tauc"].GetData ([| t; 0 |], [| 1; n |]) :?> single[,]
tauc[0, *]
with
| e ->
with e ->
Log.Error $"{e}"
Array.empty
let readUV (ds: DataSet) t l =
try
let n = ds.Dimensions["nele"].Length
let u = ds[ "u" ].GetData([| t; l; 0 |], [| 1; 1; n |]) :?> single [,,]
let v = ds[ "v" ].GetData([| t; l; 0 |], [| 1; 1; n |]) :?> single [,,]
let u = ds["u"].GetData ([| t; l; 0 |], [| 1; 1; n |]) :?> single[,,]
let v = ds["v"].GetData ([| t; l; 0 |], [| 1; 1; n |]) :?> single[,,]
Array.zip u[0, 0, *] v[0, 0, *]
with
| e ->
with e ->
Log.Error $"{e}"
Array.empty
let readUV' (ds: DataSet) t l =
readUV ds t l
|> Array.collect (fun (x, y) -> [| x; y |])
readUV ds t l |> Array.collect (fun (x, y) -> [| x; y |])
let readUVs (ds: DataSet) t l es =
try
let n = ds.Dimensions["nele"].Length
let u = ds[ "u" ].GetData([| t; l; 0 |], [| 1; 1; n |]) :?> single [,,]
let v = ds[ "v" ].GetData([| t; l; 0 |], [| 1; 1; n |]) :?> single [,,]
let u = ds["u"].GetData ([| t; l; 0 |], [| 1; 1; n |]) :?> single[,,]
let v = ds["v"].GetData ([| t; l; 0 |], [| 1; 1; n |]) :?> single[,,]
let u' = es |> Array.map (fun i -> u[0, 0, i])
let v' = es |> Array.map (fun i -> v[0, 0, i])
Array.zip u' v'
with
| e ->
with e ->
Log.Error $"{e}"
Array.empty
@@ -147,21 +127,26 @@ let readUVs (ds: DataSet) t l es =
/// <returns>Array of tuples, and empty on any errors.</returns>
let readUVRange (ds: DataSet) t l e0 en =
try
let u = ds[ "u" ].GetData([| t; l; e0 |], [| 1; 1; en |]) :?> single [,,]
let v = ds[ "v" ].GetData([| t; l; e0 |], [| 1; 1; en |]) :?> single [,,]
let u = ds["u"].GetData ([| t; l; e0 |], [| 1; 1; en |]) :?> single[,,]
let v = ds["v"].GetData ([| t; l; e0 |], [| 1; 1; en |]) :?> single[,,]
Array.zip u[0, 0, *] v[0, 0, *]
with
| e ->
Log.Error(e, "FvcomKit.Fvcom.readUVRange exception with input: time {Time}, layer {Layer}, starting elem {Elem0} and ending elem {ElemN}", t, l, e0, en)
with e ->
Log.Error (
e,
"FvcomKit.Fvcom.readUVRange exception with input: time {Time}, layer {Layer}, starting elem {Elem0} and ending elem {ElemN}",
t,
l,
e0,
en
)
Array.empty
let readOmega (ds: DataSet) t l =
try
let n = ds.Dimensions["node"].Length
let omega = ds[ "omega" ].GetData([| t; l; 0 |], [| 1; 1; n |]) :?> single [,,]
let omega = ds["omega"].GetData ([| t; l; 0 |], [| 1; 1; n |]) :?> single[,,]
omega[0, 0, *]
with
| e ->
with e ->
Log.Error $"{e}"
Array.empty
@@ -169,30 +154,27 @@ let readOmegaBlock (ds: DataSet) t =
try
let n = ds.Dimensions["node"].Length
let l = ds.Dimensions["siglev"].Length
let u = ds[ "omega" ].GetData([| t; 0; 0 |], [| 1; l; n |]) :?> single [,,]
let u = ds["omega"].GetData ([| t; 0; 0 |], [| 1; l; n |]) :?> single[,,]
u[0, *, *]
with
| e ->
with e ->
Log.Error $"{e}"
Array2D.zeroCreate 0 0
let readOmegas (ds: DataSet) t l es =
try
let n = ds.Dimensions["node"].Length
let omega = ds[ "omega" ].GetData([| t; l; 0 |], [| 1; 1; n |]) :?> single [,,]
let omega = ds["omega"].GetData ([| t; l; 0 |], [| 1; 1; n |]) :?> single[,,]
es |> Array.map (fun i -> omega[0, 0, i])
with
| e ->
with e ->
Log.Error $"{e}"
Array.empty
let readU (ds: DataSet) t l =
try
let n = ds.Dimensions["nele"].Length
let u = ds[ "u" ].GetData([| t; l; 0 |], [| 1; 1; n |]) :?> single [,,]
let u = ds["u"].GetData ([| t; l; 0 |], [| 1; 1; n |]) :?> single[,,]
u[0, 0, *]
with
| e ->
with e ->
Log.Error $"{e}"
Array.empty
@@ -200,20 +182,18 @@ let readUBlock (ds: DataSet) t =
try
let n = ds.Dimensions["nele"].Length
let l = ds.Dimensions["siglay"].Length
let u = ds[ "u" ].GetData([| t; 0; 0 |], [| 1; l; n |]) :?> single [,,]
let u = ds["u"].GetData ([| t; 0; 0 |], [| 1; l; n |]) :?> single[,,]
u[0, *, *]
with
| e ->
with e ->
Log.Error $"{e}"
Array2D.zeroCreate 0 0
let readV (ds: DataSet) t l =
try
let n = ds.Dimensions["nele"].Length
let v = ds[ "v" ].GetData([| t; l; 0 |], [| 1; 1; n |]) :?> single [,,]
let v = ds["v"].GetData ([| t; l; 0 |], [| 1; 1; n |]) :?> single[,,]
v[0, 0, *]
with
| e ->
with e ->
Log.Error $"{e}"
Array.empty
@@ -221,20 +201,18 @@ let readVBlock (ds: DataSet) t =
try
let n = ds.Dimensions["nele"].Length
let l = ds.Dimensions["siglay"].Length
let v = ds[ "v" ].GetData([| t; 0; 0 |], [| 1; l; n |]) :?> single [,,]
let v = ds["v"].GetData ([| t; 0; 0 |], [| 1; l; n |]) :?> single[,,]
v[0, *, *]
with
| e ->
with e ->
Log.Error $"{e}"
Array2D.zeroCreate 0 0
let readWw (ds: DataSet) t l =
try
let n = ds.Dimensions["nele"].Length
let ww = ds[ "ww" ].GetData([| t; l; 0 |], [| 1; 1; n |]) :?> single [,,]
let ww = ds["ww"].GetData ([| t; l; 0 |], [| 1; 1; n |]) :?> single[,,]
ww[0, 0, *]
with
| e ->
with e ->
Log.Error $"{e}"
Array.empty
@@ -242,69 +220,62 @@ let readWwBlock (ds: DataSet) t =
try
let n = ds.Dimensions["nele"].Length
let l = ds.Dimensions["siglay"].Length
let w = ds[ "ww" ].GetData([| t; 0; 0 |], [| 1; l; n |]) :?> single [,,]
let w = ds["ww"].GetData ([| t; 0; 0 |], [| 1; l; n |]) :?> single[,,]
w[0, *, *]
with
| e ->
with e ->
Log.Error $"{e}"
Array2D.zeroCreate 0 0
let readWws (ds: DataSet) t l es =
try
let n = ds.Dimensions["nele"].Length
let ww = ds[ "ww" ].GetData([| t; l; 0 |], [| 1; 1; n |]) :?> single [,,]
let ww = ds["ww"].GetData ([| t; l; 0 |], [| 1; 1; n |]) :?> single[,,]
es |> Array.map (fun i -> ww[0, 0, i])
with
| e ->
with e ->
Log.Error $"{e}"
Array.empty
let readUVW (ds: DataSet) t l =
try
let n = ds.Dimensions["nele"].Length
let u = ds[ "u" ].GetData([| t; l; 0 |], [| 1; 1; n |]) :?> single [,,]
let v = ds[ "v" ].GetData([| t; l; 0 |], [| 1; 1; n |]) :?> single [,,]
let w = ds[ "ww" ].GetData([| t; l; 0 |], [| 1; 1; n |]) :?> single [,,]
let u = ds["u"].GetData ([| t; l; 0 |], [| 1; 1; n |]) :?> single[,,]
let v = ds["v"].GetData ([| t; l; 0 |], [| 1; 1; n |]) :?> single[,,]
let w = ds["ww"].GetData ([| t; l; 0 |], [| 1; 1; n |]) :?> single[,,]
Array.zip3 u[0, 0, *] v[0, 0, *] w[0, 0, *]
with
| e ->
with e ->
Log.Error $"{e}"
Array.empty
let readUVWs (ds: DataSet) t l es =
try
let n = ds.Dimensions["nele"].Length
let u = ds[ "u" ].GetData([| t; l; 0 |], [| 1; 1; n |]) :?> single [,,]
let v = ds[ "v" ].GetData([| t; l; 0 |], [| 1; 1; n |]) :?> single [,,]
let w = ds[ "ww" ].GetData([| t; l; 0 |], [| 1; 1; n |]) :?> single [,,]
let u = ds["u"].GetData ([| t; l; 0 |], [| 1; 1; n |]) :?> single[,,]
let v = ds["v"].GetData ([| t; l; 0 |], [| 1; 1; n |]) :?> single[,,]
let w = ds["ww"].GetData ([| t; l; 0 |], [| 1; 1; n |]) :?> single[,,]
let u' = es |> Array.map (fun i -> u[0, 0, i])
let v' = es |> Array.map (fun i -> v[0, 0, i])
let w' = es |> Array.map (fun i -> w[0, 0, i])
Array.zip3 u' v' w'
with
| e ->
with e ->
Log.Error $"{e}"
Array.empty
let readUVWRange (ds: DataSet) t l e0 en =
try
let u = ds[ "u" ].GetData([| t; l; e0 |], [| 1; 1; en |]) :?> single [,,]
let v = ds[ "v" ].GetData([| t; l; e0 |], [| 1; 1; en |]) :?> single [,,]
let w = ds[ "ww" ].GetData([| t; l; e0 |], [| 1; 1; en |]) :?> single [,,]
let u = ds["u"].GetData ([| t; l; e0 |], [| 1; 1; en |]) :?> single[,,]
let v = ds["v"].GetData ([| t; l; e0 |], [| 1; 1; en |]) :?> single[,,]
let w = ds["ww"].GetData ([| t; l; e0 |], [| 1; 1; en |]) :?> single[,,]
Array.zip3 u[0, 0, *] v[0, 0, *] w[0, 0, *]
with
| e ->
with e ->
Log.Error $"{e}"
Array.empty
let readTemp (ds: DataSet) t l =
try
let n = ds.Dimensions["node"].Length
let temp =
ds[ "temp" ].GetData([| t; l; 0 |], [| 1; 1; n |]) :?> single [,,]
let temp = ds["temp"].GetData ([| t; l; 0 |], [| 1; 1; n |]) :?> single[,,]
temp[0, 0, *]
with
| e ->
with e ->
Log.Error $"{e}"
Array.empty
@@ -312,26 +283,18 @@ let readTempBlock (ds: DataSet) t =
try
let n = ds.Dimensions["node"].Length
let l = ds.Dimensions["siglay"].Length
let sal =
ds["temp"]
.GetData([| t; 0; 0 |], [| 1; l; n |])
:?> single [,,]
let sal = ds["temp"].GetData ([| t; 0; 0 |], [| 1; l; n |]) :?> single[,,]
sal[0, *, *]
with
| e ->
with e ->
Log.Error $"{e}"
Array2D.zeroCreate 0 0
let readSalinity (ds: DataSet) t l =
try
let n = ds.Dimensions["node"].Length
let sal =
ds["salinity"]
.GetData([| t; l; 0 |], [| 1; 1; n |])
:?> single [,,]
let sal = ds["salinity"].GetData ([| t; l; 0 |], [| 1; 1; n |]) :?> single[,,]
sal[0, 0, *]
with
| e ->
with e ->
Log.Error $"{e}"
Array.empty
@@ -339,48 +302,39 @@ let readSalinityBlock (ds: DataSet) t =
try
let n = ds.Dimensions["node"].Length
let l = ds.Dimensions["siglay"].Length
let sal =
ds["salinity"]
.GetData([| t; 0; 0 |], [| 1; l; n |])
:?> single [,,]
let sal = ds["salinity"].GetData ([| t; 0; 0 |], [| 1; l; n |]) :?> single[,,]
sal[0, *, *]
with
| e ->
with e ->
Log.Error $"{e}"
Array2D.zeroCreate 0 0
let readArt1 (ds: DataSet) =
try
ds[ "art1" ].GetData() :?> single []
with
| e ->
ds["art1"].GetData () :?> single[]
with e ->
Log.Error $"{e}"
Array.empty
let readZeta (ds: DataSet) t =
try
let n = ds.Dimensions["node"].Length
let zeta = ds[ "zeta" ].GetData([| t; 0 |], [| 1; n |]) :?> single [,]
let zeta = ds["zeta"].GetData ([| t; 0 |], [| 1; n |]) :?> single[,]
zeta[0, *]
with
| e ->
with e ->
Log.Error $"{e}"
Array.empty
let readBathymetry (ds: DataSet) =
try
ds[ "h" ].GetData() :?> single []
with
| e ->
ds["h"].GetData () :?> single[]
with e ->
Log.Error e.Message
Array.empty
let tryReadBathymetry (ds: DataSet) =
if ds.Variables.Contains "h" then
ds[ "h" ].GetData()
:?> single []
|> Some
ds["h"].GetData () :?> single[] |> Some
else
None
@@ -388,11 +342,10 @@ let readBathymetryAtCenters (ds: DataSet) =
let h = readBathymetry ds
let lc = ds.Dimensions["nele"].Length
try
let nv = ds[ "nv" ].GetData() :?> int [,] |> Array2D.map (fun n -> n - 1)
let nv = ds["nv"].GetData () :?> int[,] |> Array2D.map (fun n -> n - 1)
[| 0 .. lc - 1 |]
|> Array.map (fun i -> (h[nv[0,i]] + h[nv[1,i]] + h[nv[2,i]]) / 3.0f)
with
| e ->
|> Array.map (fun i -> (h[nv[0, i]] + h[nv[1, i]] + h[nv[2, i]]) / 3.0f)
with e ->
Log.Error e.Message
Array.empty
@@ -400,25 +353,23 @@ let readBathymetryAtCenter (ds: DataSet) e =
try
let h = readBathymetryAtCenters ds
h[e]
with
| e ->
with e ->
Log.Error e.Message
-1f
let readSiglev (ds: DataSet) n =
try
let l = ds.Dimensions["siglev"].Length
let siglev = ds[ "siglev" ].GetData([| 0; n |], [| l; 1 |]) :?> single [,]
let siglev = ds["siglev"].GetData ([| 0; n |], [| l; 1 |]) :?> single[,]
siglev[*, 0]
with
| err ->
with err ->
Log.Error $"{err}"
Array.empty
let tryReadSiglev (ds: DataSet) n =
if ds.Variables.Contains "siglev" then
let l = ds.Dimensions["siglev"].Length
let siglev = ds[ "siglev" ].GetData([| 0; n |], [| l; 1 |]) :?> single [,]
let siglev = ds["siglev"].GetData ([| 0; n |], [| l; 1 |]) :?> single[,]
Some siglev[*, 0]
else
None
@@ -426,55 +377,50 @@ let tryReadSiglev (ds: DataSet) n =
let readSiglay (ds: DataSet) n =
try
let l = ds.Dimensions["siglay"].Length
let siglay = ds[ "siglay" ].GetData([| 0; n |], [| l; 1 |]) :?> single [,]
let siglay = ds["siglay"].GetData ([| 0; n |], [| l; 1 |]) :?> single[,]
siglay[*, 0]
with
| err ->
with err ->
Log.Error $"{err}"
Array.empty
let tryReadSiglay (ds: DataSet) n =
if ds.Variables.Contains "siglay" then
let l = ds.Dimensions["siglay"].Length
let siglay = ds[ "siglay" ].GetData([| 0; n |], [| l; 1 |]) :?> single [,]
let siglay = ds["siglay"].GetData ([| 0; n |], [| l; 1 |]) :?> single[,]
Some siglay[*, 0]
else
None
let readSiglevAtCenter (ds: DataSet) e =
try
let nv = ds[ "nv" ].GetData() :?> int [,] |> Array2D.map (fun n -> n - 1)
let nv = ds["nv"].GetData () :?> int[,] |> Array2D.map (fun n -> n - 1)
let s1 = readSiglev ds nv[0, e]
let s2 = readSiglev ds nv[1, e]
let s3 = readSiglev ds nv[2, e]
[| 0 .. s1.Length - 1 |]
|> Array.map ( fun i -> (s1[i] + s2[i] + s3[i]) / 3.0f)
with
| e ->
[| 0 .. s1.Length - 1 |] |> Array.map (fun i -> (s1[i] + s2[i] + s3[i]) / 3.0f)
with e ->
Log.Error $"{e}"
Array.empty
let readSiglayAtCenter (ds: DataSet) e =
try
let nv = ds[ "nv" ].GetData() :?> int [,] |> Array2D.map (fun n -> n - 1)
let nv = ds["nv"].GetData () :?> int[,] |> Array2D.map (fun n -> n - 1)
let s1 = readSiglay ds nv[0, e]
let s2 = readSiglay ds nv[1, e]
let s3 = readSiglay ds nv[2, e]
[| 0 .. s1.Length - 1 |]
|> Array.map ( fun i -> (s1[i] + s2[i] + s3[i]) / 3.0f)
with
| e ->
[| 0 .. s1.Length - 1 |] |> Array.map (fun i -> (s1[i] + s2[i] + s3[i]) / 3.0f)
with e ->
Log.Error $"{e}"
Array.empty
let tryReadSiglayAtCenter (ds: DataSet) e =
if ds.Variables.Contains "siglay" then
let nv = ds[ "nv" ].GetData() :?> int [,] |> Array2D.map (fun n -> n - 1)
let nv = ds["nv"].GetData () :?> int[,] |> Array2D.map (fun n -> n - 1)
let s1 = readSiglay ds nv[0, e]
let s2 = readSiglay ds nv[1, e]
let s3 = readSiglay ds nv[2, e]
[| 0 .. s1.Length - 1 |]
|> Array.map ( fun i -> (s1[i] + s2[i] + s3[i]) / 3.0f)
|> Array.map (fun i -> (s1[i] + s2[i] + s3[i]) / 3.0f)
|> Some
else
None
@@ -485,49 +431,43 @@ let readSiglayCenter = readSiglayAtCenter
module Siglay =
let readSiglay (ds: DataSet) =
try
ds[ "siglay" ].GetData() :?> single [,]
with
| err ->
ds["siglay"].GetData () :?> single[,]
with err ->
Log.Error $"{err}"
Array2D.zeroCreate 0 0
let tryReadSiglay (ds: DataSet) =
if ds.Variables.Contains "siglay" then
ds[ "siglay" ].GetData()
:?> single [,]
|> Some
ds["siglay"].GetData () :?> single[,] |> Some
else
None
let readSiglayAtCenter (ds: DataSet) =
let siglay = readSiglay ds
try
let nv = ds[ "nv" ].GetData() :?> int [,] |> Array2D.map (fun n -> n - 1)
let nv = ds["nv"].GetData () :?> int[,] |> Array2D.map (fun n -> n - 1)
let l1 = Array2D.length1 siglay
let l2 = Array2D.length2 nv
let sc = Array2D.zeroCreate l1 l2
sc
|> Array2D.mapi (fun i j _ -> (siglay[i,nv[0,j]] + siglay[i,nv[1,j]] + siglay[i,nv[2,j]]) / 3.0f)
with
| e ->
|> Array2D.mapi (fun i j _ -> (siglay[i, nv[0, j]] + siglay[i, nv[1, j]] + siglay[i, nv[2, j]]) / 3.0f)
with e ->
Log.Error $"{e}"
Array2D.zeroCreate 0 0
let tryReadSiglayAtCenter (ds: DataSet) =
if ds.Variables.Contains "siglay" then
readSiglayAtCenter ds
|> Some
readSiglayAtCenter ds |> Some
else
None
let readUv (ds: DataSet) e t =
try
let l = ds.Dimensions["siglay"].Length
let u = ds[ "u" ].GetData([| t; 0; e |], [| 1; l; 1 |]) :?> single [,,]
let v = ds[ "v" ].GetData([| t; 0; e |], [| 1; l; 1 |]) :?> single [,,]
let u = ds["u"].GetData ([| t; 0; e |], [| 1; l; 1 |]) :?> single[,,]
let v = ds["v"].GetData ([| t; 0; e |], [| 1; l; 1 |]) :?> single[,,]
Array.zip u[0, *, 0] v[0, *, 0]
with
| err ->
with err ->
Log.Warning $"readUv {e} {t}"
Log.Error $"{err}"
Array.empty
@@ -535,12 +475,11 @@ module Siglay =
let readUvw (ds: DataSet) e t =
try
let l = ds.Dimensions["siglay"].Length
let u = ds[ "u" ].GetData([| t; 0; e |], [| 1; l; 1 |]) :?> single [,,]
let v = ds[ "v" ].GetData([| t; 0; e |], [| 1; l; 1 |]) :?> single [,,]
let w = ds[ "ww" ].GetData([| t; 0; e |], [| 1; l; 1 |]) :?> single [,,]
let u = ds["u"].GetData ([| t; 0; e |], [| 1; l; 1 |]) :?> single[,,]
let v = ds["v"].GetData ([| t; 0; e |], [| 1; l; 1 |]) :?> single[,,]
let w = ds["ww"].GetData ([| t; 0; e |], [| 1; l; 1 |]) :?> single[,,]
Array.zip3 u[0, *, 0] v[0, *, 0] w[0, *, 0]
with
| err ->
with err ->
Log.Warning $"readUv {e} {t}"
Log.Error $"{err}"
Array.empty
@@ -548,12 +487,10 @@ module Siglay =
let readUv' (ds: DataSet) e t =
try
let l = ds.Dimensions["siglay"].Length
let u = ds[ "u" ].GetData([| t; 0; e |], [| 1; l; 1 |]) :?> single [,,]
let v = ds[ "v" ].GetData([| t; 0; e |], [| 1; l; 1 |]) :?> single [,,]
Array.zip u[0, *, 0] v[0, *, 0]
|> Array.collect (fun (x, y) -> [| x; y |])
with
| err ->
let u = ds["u"].GetData ([| t; 0; e |], [| 1; l; 1 |]) :?> single[,,]
let v = ds["v"].GetData ([| t; 0; e |], [| 1; l; 1 |]) :?> single[,,]
Array.zip u[0, *, 0] v[0, *, 0] |> Array.collect (fun (x, y) -> [| x; y |])
with err ->
Log.Warning $"readUv' {e} {t}"
Log.Error $"{err}"
Array.empty
@@ -561,10 +498,9 @@ module Siglay =
let readTemp (ds: DataSet) n t =
try
let l = ds.Dimensions["siglay"].Length
let t = ds[ "temp" ].GetData([| t; 0; n |], [| 1; l; 1 |]) :?> single [,,]
let t = ds["temp"].GetData ([| t; 0; n |], [| 1; l; 1 |]) :?> single[,,]
t[0, *, 0]
with
| err ->
with err ->
Log.Warning $"readTemp {n} {t}"
Log.Error $"{err}"
Array.empty
@@ -572,13 +508,9 @@ module Siglay =
let readSalinity (ds: DataSet) n t =
try
let l = ds.Dimensions["siglay"].Length
let t =
ds["salinity"]
.GetData([| t; 0; n |], [| 1; l; 1 |])
:?> single [,,]
let t = ds["salinity"].GetData ([| t; 0; n |], [| 1; l; 1 |]) :?> single[,,]
t[0, *, 0]
with
| err ->
with err ->
Log.Warning $"readS {n} {t}"
Log.Error $"{err}"
Array.empty
@@ -586,49 +518,44 @@ module Siglay =
module Singular =
let readUv (ds: DataSet) e t d =
try
let u = ds[ "u" ].GetData([| t; d; e |], [| 1; 1; 1 |]) :?> single [,,]
let v = ds[ "v" ].GetData([| t; d; e |], [| 1; 1; 1 |]) :?> single [,,]
let u = ds["u"].GetData ([| t; d; e |], [| 1; 1; 1 |]) :?> single[,,]
let v = ds["v"].GetData ([| t; d; e |], [| 1; 1; 1 |]) :?> single[,,]
u[0, 0, 0], v[0, 0, 0]
with
| err ->
with err ->
Log.Warning $"readUv {e} {t}"
Log.Error $"{err}"
0f, 0f
let readTemp (ds: DataSet) n t d =
try
let t = ds[ "temp" ].GetData([| t; d; n |], [| 1; 1; 1 |]) :?> single [,,]
let t = ds["temp"].GetData ([| t; d; n |], [| 1; 1; 1 |]) :?> single[,,]
t[0, 0, 0]
with
| err ->
with err ->
Log.Warning $"readTemp {n} {t}"
Log.Error $"{err}"
0f
let readSalinity (ds: DataSet) n t d =
try
let t =
ds["salinity"].GetData([| t; d; n |], [| 1; 1; 1 |]) :?> single [,,]
let t = ds["salinity"].GetData ([| t; d; n |], [| 1; 1; 1 |]) :?> single[,,]
t[0, 0, 0]
with
| err ->
with err ->
Log.Warning $"readS {n} {t}"
Log.Error $"{err}"
0f
let readZeta (ds: DataSet) n t =
try
let zeta = ds[ "zeta" ].GetData([| t; n |], [| 1; 1 |]) :?> single [,]
let zeta = ds["zeta"].GetData ([| t; n |], [| 1; 1 |]) :?> single[,]
zeta[0, 0]
with
| e ->
with e ->
Log.Error $"{e}"
0f
let getBBox (ds: DataSet) : BBox =
try
let x = ds[ "x" ].GetData() :?> single []
let y = ds[ "y" ].GetData() :?> single []
let x = ds["x"].GetData () :?> single[]
let y = ds["y"].GetData () :?> single[]
let minX = Array.min x
let maxX = Array.max x
let minY = Array.min y
@@ -641,26 +568,17 @@ let getBBox (ds: DataSet) : BBox =
maxY = float maxY
center = center
}
with
| e ->
with e ->
Log.Error $"{e}"
BBox.empty
let getGrid (ds: DataSet) : FvcomGrid =
try
let x =
ds[ "x" ].GetData() :?> single []
|> Array.map float
let y =
ds[ "y" ].GetData() :?> single []
|> Array.map float
let xc =
ds[ "xc" ].GetData() :?> single []
|> Array.map float
let yc =
ds[ "yc" ].GetData() :?> single []
|> Array.map float
let nv = ds[ "nv" ].GetData() :?> int [,]
let x = ds["x"].GetData () :?> single[] |> Array.map float
let y = ds["y"].GetData () :?> single[] |> Array.map float
let xc = ds["xc"].GetData () :?> single[] |> Array.map float
let yc = ds["yc"].GetData () :?> single[] |> Array.map float
let nv = ds["nv"].GetData () :?> int[,]
let h =
match tryReadBathymetry ds with
@@ -676,7 +594,7 @@ let getGrid (ds: DataSet) : FvcomGrid =
| None -> Array2D.zeroCreate 0 0
let siglev =
if ds.Variables.Contains "siglev" then
ds[ "siglev" ].GetData() :?> single [,]
ds["siglev"].GetData () :?> single[,]
else
Array2D.zeroCreate 0 0
let elem =
@@ -694,26 +612,17 @@ let getGrid (ds: DataSet) : FvcomGrid =
SiglayCenter = siglay_c
Siglev = siglev
}
with
| e ->
with e ->
Log.Error $"{e}"
FvcomGrid.empty
let getGridLonLat (ds: DataSet) : FvcomGrid =
try
let x =
ds[ "lon" ].GetData() :?> single []
|> Array.map float
let y =
ds[ "lat" ].GetData() :?> single []
|> Array.map float
let xc =
ds[ "lonc" ].GetData() :?> single []
|> Array.map float
let yc =
ds[ "latc" ].GetData() :?> single []
|> Array.map float
let nv = ds[ "nv" ].GetData() :?> int [,]
let x = ds["lon"].GetData () :?> single[] |> Array.map float
let y = ds["lat"].GetData () :?> single[] |> Array.map float
let xc = ds["lonc"].GetData () :?> single[] |> Array.map float
let yc = ds["latc"].GetData () :?> single[] |> Array.map float
let nv = ds["nv"].GetData () :?> int[,]
let h =
match tryReadBathymetry ds with
@@ -729,7 +638,7 @@ let getGridLonLat (ds: DataSet) : FvcomGrid =
| None -> Array2D.zeroCreate 0 0
let siglev =
if ds.Variables.Contains "siglev" then
ds[ "siglev" ].GetData() :?> single [,]
ds["siglev"].GetData () :?> single[,]
else
Array2D.zeroCreate 0 0
let elem =
@@ -747,24 +656,20 @@ let getGridLonLat (ds: DataSet) : FvcomGrid =
SiglayCenter = siglay_c
Siglev = siglev
}
with
| e ->
with e ->
Log.Error $"{e}"
FvcomGrid.empty
let getNbve (ds: DataSet) =
let nbve = ds[ "nbve" ].GetData() :?> int [,]
let nbve = ds["nbve"].GetData () :?> int[,]
[|
for i = 0 to (Array2D.length2 nbve - 1) do
nbve[*, i]
|> Array.filter ((<>) 0)
|> Array.map (fun x -> x - 1)
|> Array.rev
nbve[*, i] |> Array.filter ((<>) 0) |> Array.map (fun x -> x - 1) |> Array.rev
|]
let projectFvcomGrid proj (grid: FvcomGrid) : FvcomGrid =
{ grid with
let projectFvcomGrid proj (grid: FvcomGrid) : FvcomGrid = {
grid with
Nodes = grid.Nodes |> Array.Parallel.map proj
BBox = projectBBox proj grid.BBox
Cells = grid.Cells |> Array.Parallel.map proj
}
}

8
src/Grid.fs
View File

@@ -525,7 +525,7 @@ type ExtendedGrid(grid: IGrid) =
let mutable elementTree: IdxTree option = None
let mutable neighborIndex: NeighborIndex option = None
let mutable centroids: Vertex [] option = None
let mutable gridHash: byte[] = [||]
let mutable gridSha1: byte[] = [||]
let getNeighborIdx () =
match neighborIndex with
@@ -614,7 +614,7 @@ type ExtendedGrid(grid: IGrid) =
Util.tryFindElement grid elementTree.Value p
member this.getSha1() =
if gridHash.Length = 0 then
if gridSha1.Length = 0 then
let bg : BinGrid = {
hash = [||]
vertices = (this :> IGrid).getVertices ()
@@ -622,8 +622,8 @@ type ExtendedGrid(grid: IGrid) =
}
let bytes = MessagePackSerializer.Serialize(bg)
let sha1 = System.Security.Cryptography.SHA1.Create()
gridHash <- sha1.ComputeHash bytes
gridHash
gridSha1 <- sha1.ComputeHash bytes
gridSha1
member this.getSha1String() =
this.getSha1() |> Convert.ToHexStringLower

76
src/Interpol.fs
View File

@@ -8,8 +8,8 @@ type InterpolCoefs = ((int * int) * (float * float)) array array
type DepthInterpolCoefs = { iRho: InterpolCoefs; iU: InterpolCoefs; iV: InterpolCoefs }
let private findNearestZ z0 (h: float []) =
let rec findNerest' z0 (h: float []) d n =
let private findNearestZ z0 (h: float[]) =
let rec findNerest' z0 (h: float[]) d n =
if n < h.Length - 1 then
let d' = abs (h[n] - z0)
if d' < d then findNerest' z0 h d' (n + 1) else n
@@ -40,42 +40,36 @@ let private calcInterpolationWeightedIdx (rn, fn) =
else
failwith "not reachable")
let private calcInterpMatrices (rz: float [] []) (fz: float [] []) =
Array.zip rz fz
|> Array.Parallel.map calcInterpolationWeightedIdx
let private calcInterpMatrices (rz: float[][]) (fz: float[][]) =
Array.zip rz fz |> Array.Parallel.map calcInterpolationWeightedIdx
let mkDepthInterpolCoefs (siglay: single [,]) (h: float []) (roms: float [] []) =
let s =
siglay
|> Array2D.map float
|> Matrix.Build.DenseOfArray
let mkDepthInterpolCoefs (siglay: single[,]) (h: float[]) (roms: float[][]) =
let s = siglay |> Array2D.map float |> Matrix.Build.DenseOfArray
let rescale (sigma: float Matrix) h =
let sh = h |> Array.map float
sigma.MapIndexed(fun _ j x -> -x * sh[j])
|> fun x -> x.ToColumnArrays()
sigma.MapIndexed (fun _ j x -> -x * sh[j]) |> fun x -> x.ToColumnArrays ()
let conv m =
m
|> Array.map (Array.map ((*) -1.0))
|> matrix
|> fun x -> x.ToColumnArrays()
|> fun x -> x.ToColumnArrays ()
let romz = conv roms
let z = rescale s h
calcInterpMatrices romz z
let zInterpolProp (iz: InterpolCoefs) (adjRomsProp: float [] []) =
let zInterpolProp (iz: InterpolCoefs) (adjRomsProp: float[][]) =
let pz = adjRomsProp |> matrix
iz
|> Array.Parallel.mapi (fun n x ->
let p = pz[ *, n ].ToArray() |> Array.rev
x
|> Array.map (fun ((i1, i2), (w1, w2)) -> p[i1] * w1 + p[i2] * w2))
let p = pz[*, n].ToArray () |> Array.rev
x |> Array.map (fun ((i1, i2), (w1, w2)) -> p[i1] * w1 + p[i2] * w2))
|> matrix
|> fun x -> x.Transpose().ToArray()
|> fun x -> x.Transpose().ToArray ()
// Trivially already satisfied by ajoint
let hInterpolNearestProp (adjRomsProp: float [] []) = adjRomsProp
let hInterpolNearestProp (adjRomsProp: float[][]) = adjRomsProp
type BiWght = float * float * float * float
type BiProp = float * float * float * float
@@ -92,47 +86,43 @@ let genBilinearInterpolationWgts
let CD = (x3 - x2) ** 2.0 + (y3 - y2) ** 2.0 |> sqrt
let DA = (x0 - x3) ** 2.0 + (y0 - y3) ** 2.0 |> sqrt
let h0 =
let a0 =
(x * (y0 - y3) + x0 * (y3 - y) + x3 * (y - y0))
* 0.5
let a0 = (x * (y0 - y3) + x0 * (y3 - y) + x3 * (y - y0)) * 0.5
2.0 * a0 / DA
let h1 =
let a1 =
(x * (y1 - y0) + x0 * (y - y1) + x1 * (y0 - y))
* 0.5
let a1 = (x * (y1 - y0) + x0 * (y - y1) + x1 * (y0 - y)) * 0.5
2.0 * a1 / AB
let w00 = (h0 / AB - 1.0) * (h1 / DA - 1.0)
let w10 = h0 / AB * (1.0 - h1 / DA)
let w11 = h1 / DA * h0 / CD
let w01 = h1 / DA * (1.0 - h0 / CD)
Some(w00, w10, w11, w01)
Some (w00, w10, w11, w01)
else
None
let makeBiWeights
(pos: (float * float) [])
(boxes: ((float * float) * (float * float) * (float * float) * (float * float)) [])
(mask: (bool * bool * bool * bool) [])
(pos: (float * float)[])
(boxes: ((float * float) * (float * float) * (float * float) * (float * float))[])
(mask: (bool * bool * bool * bool)[])
=
Array.zip3 pos boxes mask
|> Array.map (fun (p, b, m) -> genBilinearInterpolationWgts p b m)
let interpolateCells (wghts: BiWght option []) (prop: BiProp []) =
let interpolateCells (wghts: BiWght option[]) (prop: BiProp[]) =
Array.zip wghts prop
|> Array.Parallel.map (fun (w, (p00, p01, p11, p10)) ->
match w with
| Some w ->
let w00, w01, w11, w10 = w
Some(w00 * p00 + w01 * p01 + w11 * p11 + w10 * p10)
Some (w00 * p00 + w01 * p01 + w11 * p11 + w10 * p10)
| None -> None)
// |> Array.Parallel.map (fun ((w00, w01, w11, w10), (p00, p01, p11, p10)) ->
// w00 * p00 + w01 * p01 + w11 * p11 + w10 * p10)
let private valtest (intval: float option []) =
let private valtest (intval: float option[]) =
let somval = intval |> Array.filter (fun v -> Option.isSome v)
somval.Length = intval.Length
let private fillOutOfBounds (grid: Grid) (interpVal: float option []) =
let private fillOutOfBounds (grid: Grid) (interpVal: float option[]) =
let node = interpVal.Length = grid.Nodes.Length
let nbridx = makeNeighborIndex (grid :> IGrid)
let mutable allval = valtest interpVal
@@ -141,24 +131,20 @@ let private fillOutOfBounds (grid: Grid) (interpVal: float option []) =
if Option.isNone interpVal[n] then
let nb =
if node then
nbridx.NodesAroundNode[n]
|> Array.filter (fun k -> Option.isSome interpVal[k])
nbridx.NodesAroundNode[n] |> Array.filter (fun k -> Option.isSome interpVal[k])
else
getElemsSurroundingElem nbridx grid n
|> Array.filter (fun k -> Option.isSome interpVal[k])
if nb.Length > 0 then
let nval =
nb
|> Array.map (fun k -> interpVal[k].Value)
|> Array.average
let nval = nb |> Array.map (fun k -> interpVal[k].Value) |> Array.average
interpVal[n] <- Some nval
allval <- valtest interpVal
interpVal |> Array.map (fun v -> v.Value)
let interpolateProp
(coords, _ as grid: BiPos * Mask)
(pos: (float * float) [])
(prop: float [,])
(pos: (float * float)[])
(prop: float[,])
(fvgrid: Grid)
(oob: bool)
=
@@ -189,7 +175,7 @@ let interpolateProp
else
iprop |> Array.map (fun v -> v.Value)
let interpCoefs (coords, _ as grid: BiPos * Mask) (pos: (float * float) []) =
let interpCoefs (coords, _ as grid: BiPos * Mask) (pos: (float * float)[]) =
let cm = getNearestCellCorner grid pos
let boxes, mask =
cm
@@ -198,12 +184,10 @@ let interpCoefs (coords, _ as grid: BiPos * Mask) (pos: (float * float) []) =
|> Array.unzip
let wgths = makeBiWeights pos boxes mask
let corneridx =
cm
|> Array.map (fun c -> Array.unzip c)
|> Array.map (fun (i, _) -> i)
cm |> Array.map (fun c -> Array.unzip c) |> Array.map (fun (i, _) -> i)
corneridx, wgths
let interpProp (corneridx: (int * int) [] []) (wgths: BiWght option []) (prop: float [,]) (grid: Grid) (oob: bool) =
let interpProp (corneridx: (int * int)[][]) (wgths: BiWght option[]) (prop: float[,]) (grid: Grid) (oob: bool) =
let boxProps =
corneridx
|> Array.map (fun box -> box |> Array.map (fun (n, m) -> prop[n, m]))

8
src/NorKyst.fs
View File

@@ -10,8 +10,8 @@ let private getNorkystUrl (threddsUrl: string) (d: DateTime) =
$"{url}/fou-hi/new_norkyst800m/his/ocean_his.an.{fmt d}"
let private getArchiveUrls urlf (t: DateTime) =
let t = t.ToUniversalTime()
let now = DateTime.Now.ToUniversalTime()
let t = t.ToUniversalTime ()
let now = DateTime.Now.ToUniversalTime ()
let dDay = (now.Date - t.Date).Days
if dDay < 0 then // no data available
[]
@@ -21,6 +21,4 @@ let private getArchiveUrls urlf (t: DateTime) =
let tryGetArchive (threddsUrl: string) (t: DateTime) =
getArchiveUrls (getNorkystUrl threddsUrl) t
|> tryOpenThredds
|> Option.bind (fun (_, ds) ->
tryGetTimeIndex ds t
|> Option.map (fun idx -> (ds, idx)))
|> Option.bind (fun (_, ds) -> tryGetTimeIndex ds t |> Option.map (fun idx -> (ds, idx)))

11
src/NorShelf.fs
View File

@@ -13,8 +13,8 @@ let private getNorshelfUrl (threddsUrl: string) (kind: Kind) (mode: Mode) (date:
$"{url}/{date.Year}/%02d{date.Month}/norshelf_{kind}_{mode}_{fmt date}"
let private getArchiveUrls urlf (t: DateTime) =
let t = t.ToUniversalTime()
let now = DateTime.Now.ToUniversalTime()
let t = t.ToUniversalTime ()
let now = DateTime.Now.ToUniversalTime ()
let dDay = (now.Date - t.Date).Days
if dDay < -3 then // no data available
[]
@@ -29,9 +29,7 @@ let tryGetArchive threddsUrl avg (t: DateTime) =
let kind = if avg then Avg else Qck
getArchiveUrls (getNorshelfUrl threddsUrl kind) t
|> tryOpenThredds
|> Option.bind (fun (fname, ds) ->
tryGetTimeIndex ds t
|> Option.map (fun idx -> (fname, ds, idx)))
|> Option.bind (fun (fname, ds) -> tryGetTimeIndex ds t |> Option.map (fun idx -> (fname, ds, idx)))
let readArchive file reader =
tryOpenArchive file
@@ -40,8 +38,7 @@ let readArchive file reader =
let data = reader nc
Log.Information $"Read NorShelf data from {file}"
Some data
with
| e ->
with e ->
Log.Error e.Message
None)

2
src/Oceanbox.FvcomKit.fsproj
View File

@@ -37,6 +37,6 @@
<PackageReference Include="Serilog.Sinks.Console" Version="6.0.0" />
<PackageReference Include="Serilog.Sinks.Seq" Version="9.0.0" />
<PackageReference Include="Thoth.Json.Net" Version="12.0.0" />
<PackageReference Update="FSharp.Core" Version="9.0.200" />
<PackageReference Update="FSharp.Core" Version="9.0.201" />
</ItemGroup>
</Project>

10
src/Polygon.fs
View File

@@ -2,7 +2,7 @@
type Point = float * float
type Polygon = Point []
type Polygon = Point[]
let private maxdist (pol: Polygon) (pnt: Point) =
pol
@@ -28,10 +28,7 @@ let private onsegment (p: Polygon) =
let maxy = y[0..1] |> Array.max
let miny = y[0..1] |> Array.min
if x[2] >= minx
&& x[2] <= maxx
&& y[2] >= miny
&& y[2] <= maxy then
if x[2] >= minx && x[2] <= maxx && y[2] >= miny && y[2] <= maxy then
true
else
false
@@ -56,8 +53,7 @@ let private intersect (p1: Polygon) (p2: Polygon) =
let private isodd a = if a % 2 = 0 then false else true
let inpolygon (pol: Polygon) (point: Point) =
let (line: Polygon) =
[| fst point, snd point; fst point + (maxdist pol point), snd point |]
let (line: Polygon) = [| fst point, snd point; fst point + (maxdist pol point), snd point |]
let poly =
if pol[0] <> pol[pol.Length - 1] then
Array.append pol [| pol[0] |] |> Array.pairwise

170
src/ROMS.fs
View File

@@ -11,36 +11,33 @@ open Adjoin
open Interpol
type GridData<'a> = { rho: 'a; u: 'a; v: 'a }
type GridMask = GridData<bool [,]>
type GridMask = GridData<bool[,]>
type GridPos = GridData<BiPos>
type GridVars = GridData<float [,]>
type GridVars = GridData<float[,]>
type RomsGrid =
{
h: GridVars
z: GridVars []
pos: GridPos
wetMask: GridMask
angle: float [,]
}
type RomsGrid = {
h: GridVars
z: GridVars[]
pos: GridPos
wetMask: GridMask
angle: float[,]
}
type RomsProps =
{
salt: float [,] []
temp: float [,] []
zeta: float [,]
u: float [,] []
v: float [,] []
}
type RomsProps = {
salt: float[,][]
temp: float[,][]
zeta: float[,]
u: float[,][]
v: float[,][]
}
type AdjoinedProp =
{
salt: float [] []
temp: float [] []
zeta: float []
u: float [] []
v: float [] []
}
type AdjoinedProp = {
salt: float[][]
temp: float[][]
zeta: float[]
u: float[][]
v: float[][]
}
type FvcomGrid = Fvcom.FvcomGrid
@@ -49,58 +46,45 @@ let private layerZ csR (vars: GridVars) =
let rho = csR |> Array.map (mult vars.rho)
let u = csR |> Array.map (mult vars.u)
let v = csR |> Array.map (mult vars.v)
Array.zip3 rho u v
|> Array.map (fun (a, b, c) -> { rho = a; u = b; v = c })
Array.zip3 rho u v |> Array.map (fun (a, b, c) -> { rho = a; u = b; v = c })
// convert fractional ROMS landmask to ocean mask
let readWetMask (nc: DataSet) =
let isWet = Array2D.map ((=) 1.)
{
rho = nc[ "mask_rho" ].GetData() :?> float [,] |> isWet
u = nc[ "mask_u" ].GetData() :?> float [,] |> isWet
v = nc[ "mask_v" ].GetData() :?> float [,] |> isWet
rho = nc["mask_rho"].GetData () :?> float[,] |> isWet
u = nc["mask_u"].GetData () :?> float[,] |> isWet
v = nc["mask_v"].GetData () :?> float[,] |> isWet
}
let private readProp3 (nc: DataSet) (p: string) t (prop: string) =
let nEta = nc.Dimensions[$"eta_{p}"].Length
let nXi = nc.Dimensions[$"xi_{p}"].Length
let scaling =
nc.Variables[prop].Metadata["scale_factor"] :?> single
|> float
let offset =
nc.Variables[prop].Metadata["add_offset"] :?> single
|> float
let scaling = nc.Variables[prop].Metadata["scale_factor"] :?> single |> float
let offset = nc.Variables[prop].Metadata["add_offset"] :?> single |> float
let o = [| t; 0; 0 |]
let shp = [| 1; nEta; nXi |]
let convert (x: int16 [,,]) =
x[0, *, *]
|> Array2D.map (fun x -> float x * scaling + offset)
nc[ prop ].GetData(o, shp) :?> int16 [,,]
|> convert
let convert (x: int16[,,]) =
x[0, *, *] |> Array2D.map (fun x -> float x * scaling + offset)
nc[prop].GetData (o, shp) :?> int16[,,] |> convert
let private readProp4 (nc: DataSet) (p: string) t (prop: string) =
let nEta = nc.Dimensions[$"eta_{p}"].Length
let nXi = nc.Dimensions[$"xi_{p}"].Length
let nS = nc.Dimensions["s_rho"].Length
let scaling =
nc.Variables[prop].Metadata["scale_factor"] :?> single
|> float
let offset =
nc.Variables[prop].Metadata["add_offset"] :?> single
|> float
let scaling = nc.Variables[prop].Metadata["scale_factor"] :?> single |> float
let offset = nc.Variables[prop].Metadata["add_offset"] :?> single |> float
let o = [| t; 0; 0; 0 |]
let shp = [| 1; nS; nEta; nXi |]
let convert (x: int16 [,,,]) =
let convert (x: int16[,,,]) =
x[0, *, *, *]
|> Array3D.map (fun x -> float x * scaling + offset)
|> fun m ->
[|
for z = 0 to nS - 1 do
m[z, *, *]
|]
|> fun m -> [|
for z = 0 to nS - 1 do
m[z, *, *]
|]
nc[ prop ].GetData(o, shp) :?> int16 [,,,]
|> convert
nc[prop].GetData (o, shp) :?> int16[,,,] |> convert
let readProps (nc: DataSet) t : RomsProps =
Log.Information "Reading ROMS props..."
@@ -119,29 +103,26 @@ let readVerticalGrid path =
match Thredds.tryOpenArchive path with
| Some ds -> ds
| None -> failwith "open vertical grid failed"
nc[ "s_rho" ].GetData() :?> float [] |> Array.rev
nc["s_rho"].GetData () :?> float[] |> Array.rev
let readGrid (nc: DataSet) =
Log.Information "Reading ROMS grid..."
let hRho =
nc[ "h" ].GetData() :?> float [,]
|> Matrix.Build.DenseOfArray
let hRho = nc["h"].GetData () :?> float[,] |> Matrix.Build.DenseOfArray
let srho = nc[ "s_rho" ].GetData() :?> float [] |> Array.rev
let srho = nc["s_rho"].GetData () :?> float[] |> Array.rev
let e1, e2 = hRho.ColumnCount - 2, hRho.RowCount - 2
let hU = (hRho[*, 1..] + hRho[*, ..e1]) / 2.
let hV = (hRho[1.., *] + hRho[..e2, *]) / 2.
let pos: GridPos =
{
rho = nc[ "lon_rho" ].GetData() :?> float [,], nc[ "lat_rho" ].GetData() :?> float [,]
u = nc[ "lon_u" ].GetData() :?> float [,], nc[ "lat_u" ].GetData() :?> float [,]
v = nc[ "lon_v" ].GetData() :?> float [,], nc[ "lat_v" ].GetData() :?> float [,]
}
let angle = nc[ "angle" ].GetData() :?> float [,]
let h = { rho = hRho.ToArray(); u = hU.ToArray(); v = hV.ToArray() }
let pos: GridPos = {
rho = nc["lon_rho"].GetData () :?> float[,], nc["lat_rho"].GetData () :?> float[,]
u = nc["lon_u"].GetData () :?> float[,], nc["lat_u"].GetData () :?> float[,]
v = nc["lon_v"].GetData () :?> float[,], nc["lat_v"].GetData () :?> float[,]
}
let angle = nc["angle"].GetData () :?> float[,]
let h = { rho = hRho.ToArray (); u = hU.ToArray (); v = hV.ToArray () }
{
h = h
z = layerZ srho h
@@ -152,31 +133,28 @@ let readGrid (nc: DataSet) =
type IRomsToFvcom =
abstract uv: bool
abstract rhoToFvcom: 'a [,] -> 'a []
abstract uToFvcom: 'a [,] -> 'a []
abstract vToFvcom: 'a [,] -> 'a []
abstract rhoToFvcom: 'a[,] -> 'a[]
abstract uToFvcom: 'a[,] -> 'a[]
abstract vToFvcom: 'a[,] -> 'a[]
abstract rhoAdjoint: FvcomAdjoint
abstract uAdjoint: FvcomAdjoint
abstract vAdjoint: FvcomAdjoint
type AdjoinedGrid =
{
h: float []
zRho: float [] []
zU: float [] []
zV: float [] []
angle: float []
u: float []
v: float []
}
type AdjoinedGrid = {
h: float[]
zRho: float[][]
zU: float[][]
zV: float[][]
angle: float[]
u: float[]
v: float[]
}
let adjoinGirds (r2f: IRomsToFvcom) (roms: RomsGrid) =
Log.Information "Adjoining grids."
{
h = r2f.rhoToFvcom roms.h.rho
zRho =
roms.z
|> Array.map (fun h -> r2f.rhoToFvcom h.rho)
zRho = roms.z |> Array.map (fun h -> r2f.rhoToFvcom h.rho)
zU = roms.z |> Array.map (fun h -> r2f.uToFvcom h.u)
zV = roms.z |> Array.map (fun h -> r2f.vToFvcom h.v)
angle = r2f.rhoToFvcom roms.angle
@@ -216,8 +194,7 @@ let adjoinRomsToFvcom (proj: IProj) (fvcom: Grid) (roms: RomsGrid) (props: RomsP
let computeAngles (grid: FvcomGrid) (adj: AdjoinedGrid) =
let a = adj.angle
grid.Elem
|> Array.map (fun (i, j, k) -> (a[i] + a[j] + a[k]) / 3.0)
grid.Elem |> Array.map (fun (i, j, k) -> (a[i] + a[j] + a[k]) / 3.0)
let calcUVBar (fvcom: FvcomGrid) (u, v) =
let h = fvcom.Bathymetry
@@ -231,23 +208,18 @@ let calcUVBar (fvcom: FvcomGrid) (u, v) =
[|
for z = 0 to nZ - 1 do
fvcom.Elem
|> Array.mapi (fun n (i, j, k) ->
hc[n] * (s[z, i] + s[z, j] + s[z, k] |> float)
/ 3.0)
|> Array.mapi (fun n (i, j, k) -> hc[n] * (s[z, i] + s[z, j] + s[z, k] |> float) / 3.0)
|]
|> matrix
let dz =
sigz.ToColumnArrays()
|> Array.map (
Array.pairwise
>> Array.map (fun (x0, x1) -> x1 - x0)
)
sigz.ToColumnArrays ()
|> Array.map (Array.pairwise >> Array.map (fun (x0, x1) -> x1 - x0))
|> matrix
|> fun x -> x.PointwiseAbs().Transpose()
|> fun x -> x.PointwiseAbs().Transpose ()
let bar (m: Matrix<float>) =
m.PointwiseMultiply dz
|> fun x ->
let s = x.ColumnSums()
let s = x.ColumnSums ()
s.PointwiseDivide hc
let u' = Matrix.Build.DenseOfRowArrays u
let v' = Matrix.Build.DenseOfRowArrays v
@@ -256,9 +228,7 @@ let calcUVBar (fvcom: FvcomGrid) (u, v) =
// Linear interpolation to nearest vertical neighbor
let mkAllDepthInterpolCoefs (fvcom: Fvcom.FvcomGrid) (roms: AdjoinedGrid) =
Log.Information "Computing interpolation coefficients."
let uv =
Array.zip roms.u roms.v
|> Array.map (fun (u, v) -> (u + v) * 0.5)
let uv = Array.zip roms.u roms.v |> Array.map (fun (u, v) -> (u + v) * 0.5)
{
iRho = mkDepthInterpolCoefs fvcom.Siglay roms.h roms.zRho
iU = mkDepthInterpolCoefs fvcom.SiglayCenter uv roms.zU

50
src/Smoothing.fs
View File

@@ -7,61 +7,61 @@ open ProjNet.FSharp
let private OOBVAL = -100000f
let smooth (prop: single []) source =
let smooth (prop: single[]) source =
source
|> Array.fold (fun a x -> a + prop[x]) 0f
|> fun x -> x / single source.Length
let smooth' (prop: single []) source =
let smooth' (prop: single[]) source =
source
|> Array.fold (fun (n, a) x -> if prop[x] > OOBVAL then n + 1, a + prop[x] else n, a) (0, 0f)
|> fun (n, a) -> a / single n
let smoothNodes (idx: NeighborIndex) (prop: single []) =
idx.NodesAroundNode.Values
|> Seq.toArray
|> Array.Parallel.map (smooth prop)
let smoothNodes (idx: NeighborIndex) (prop: single[]) =
idx.NodesAroundNode.Values |> Seq.toArray |> Array.Parallel.map (smooth prop)
let smoothNodes2D (idx: NeighborIndex) (prop: single [,]) =
let smoothNodes2D (idx: NeighborIndex) (prop: single[,]) =
for k = 0 to Array2D.length1 prop - 1 do
Log.Debug $"2D smoothing nodes z = {k}"
prop[k, *] <- idx.NodesAroundNode.Values
|> Seq.toArray
|> Array.Parallel.map (smooth prop[k, *])
prop[k, *] <-
idx.NodesAroundNode.Values
|> Seq.toArray
|> Array.Parallel.map (smooth prop[k, *])
prop
let smoothNodes3D (idx: NeighborIndex) (prop: single [,,]) =
let smoothNodes3D (idx: NeighborIndex) (prop: single[,,]) =
for k = 0 to Array3D.length1 prop - 1 do
for n = 0 to Array3D.length2 prop - 1 do
Log.Debug $"3D smoothing nodes z = {n}"
prop[k, n, *] <- idx.NodesAroundNode.Values
|> Seq.toArray
|> Array.Parallel.map (smooth prop[k, n, *])
prop[k, n, *] <-
idx.NodesAroundNode.Values
|> Seq.toArray
|> Array.Parallel.map (smooth prop[k, n, *])
prop
let smoothElements (idx: NeighborIndex) (grid: Grid) (prop: single []) =
let smoothElements (idx: NeighborIndex) (grid: Grid) (prop: single[]) =
grid.Elem
|> Array.Parallel.mapi (fun n _ ->
let elems = getElemsSurroundingElem idx grid n
smooth prop elems)
let smoothElements2D (idx: NeighborIndex) (grid: Grid) (prop: single [,]) =
let smoothElements2D (idx: NeighborIndex) (grid: Grid) (prop: single[,]) =
for k = 0 to Array2D.length1 prop - 1 do
Log.Debug $"2D smoothing elements z = {k}"
prop[k, *] <- smoothElements idx grid prop[k, *]
prop
let smoothElements3D (idx: NeighborIndex) (grid: Grid) (prop: single [,,]) =
let smoothElements3D (idx: NeighborIndex) (grid: Grid) (prop: single[,,]) =
for k = 0 to Array3D.length1 prop - 1 do
for n = 0 to Array3D.length2 prop - 1 do
Log.Debug $"3D smoothing elements z = {n}"
prop[k, n, *] <- smoothElements idx grid prop[k, n, *]
prop
let private rectifyOob (getSurrounding: int -> int []) (oob: int []) (prop: single []) =
let private rectifyOob (getSurrounding: int -> int[]) (oob: int[]) (prop: single[]) =
for i in oob do
prop[i] <- OOBVAL // reset
let rec rectify (oob: int []) =
let rec rectify (oob: int[]) =
Log.Debug $"rectify oob remaining {oob.Length}"
for n in oob do
let ns = getSurrounding n
@@ -73,34 +73,34 @@ let private rectifyOob (getSurrounding: int -> int []) (oob: int []) (prop: sing
rectify oob
prop
let rectifyOutOfBoundsNodes (idx: NeighborIndex) (oob: int []) (prop: single []) =
let rectifyOutOfBoundsNodes (idx: NeighborIndex) (oob: int[]) (prop: single[]) =
let f = getNodesSurroundingNode idx
rectifyOob f oob prop
let rectifyOutOfBoundsNodes2D (idx: NeighborIndex) (oob: int []) (prop: single [,]) =
let rectifyOutOfBoundsNodes2D (idx: NeighborIndex) (oob: int[]) (prop: single[,]) =
let f = getNodesSurroundingNode idx
for k = 0 to Array2D.length1 prop - 1 do
prop[k, *] <- rectifyOob f oob prop[k, *]
prop
let rectifyOutOfBoundsNodes3D (idx: NeighborIndex) (oob: int []) (prop: single [,,]) =
let rectifyOutOfBoundsNodes3D (idx: NeighborIndex) (oob: int[]) (prop: single[,,]) =
let f = getNodesSurroundingNode idx
for k = 0 to Array3D.length1 prop - 1 do
for n = 0 to Array3D.length2 prop - 1 do
prop[k, n, *] <- rectifyOob f oob prop[k, n, *]
prop
let rectifyOutOfBoundsElements (idx: NeighborIndex) (grid: Grid) (oob: int []) (prop: single []) =
let rectifyOutOfBoundsElements (idx: NeighborIndex) (grid: Grid) (oob: int[]) (prop: single[]) =
let f = getElemsSurroundingElem idx grid
rectifyOob f oob prop
let rectifyOutOfBoundsElements2D (idx: NeighborIndex) (grid: Grid) (oob: int []) (prop: single [,]) =
let rectifyOutOfBoundsElements2D (idx: NeighborIndex) (grid: Grid) (oob: int[]) (prop: single[,]) =
let f = getElemsSurroundingElem idx grid
for k = 0 to Array2D.length1 prop - 1 do
prop[k, *] <- rectifyOob f oob prop[k, *]
prop
let rectifyOutOfBoundsElements3D (idx: NeighborIndex) (grid: Grid) (oob: int []) (prop: single [,,]) =
let rectifyOutOfBoundsElements3D (idx: NeighborIndex) (grid: Grid) (oob: int[]) (prop: single[,,]) =
let f = getElemsSurroundingElem idx grid
for k = 0 to Array3D.length1 prop - 1 do
for n = 0 to Array3D.length2 prop - 1 do

17
src/Thredds.fs
View File

@@ -23,18 +23,17 @@ type Kind =
let dateRange (start: DateTime) days =
if days < 0 then
List.unfold (fun d -> if d < days then None else Some(start.AddDays d, d - 1)) 0
List.unfold (fun d -> if d < days then None else Some (start.AddDays d, d - 1)) 0
else
List.unfold (fun d -> if d > days then None else Some(start.AddDays d, d + 1)) 0
List.unfold (fun d -> if d > days then None else Some (start.AddDays d, d + 1)) 0
let tryOpenArchive url =
let uri = NetCDFUri()
let uri = NetCDFUri ()
uri.Url <- url
try
let ds = NetCDFDataSet.Open uri
Some ds
with
| e ->
with e ->
Log.Debug e.Message
None
@@ -45,13 +44,13 @@ let rec tryOpenThredds =
match tryOpenArchive x with
| Some ds ->
Log.Debug $"thredds: {x}"
Some(x, ds)
Some (x, ds)
| None -> tryOpenThredds xs
let tryGetTimeIndex (ds: DataSet) (t: DateTime) =
let ot = ds[ "ocean_time" ].GetData() :?> double []
let t0 = DateTimeOffset.FromUnixTimeSeconds(ot[0] |> int64)
let t1 = DateTimeOffset.FromUnixTimeSeconds(ot[^0] |> int64)
let ot = ds["ocean_time"].GetData () :?> double[]
let t0 = DateTimeOffset.FromUnixTimeSeconds (ot[0] |> int64)
let t1 = DateTimeOffset.FromUnixTimeSeconds (ot[^0] |> int64)
Log.Debug $"t={t} t0={t0} tn={t1}"
if t < t0.DateTime || t > t1.DateTime then
Log.Error "time is out of bounds"

34
src/Types.fs
View File

@@ -4,25 +4,23 @@ open System
type Vertex = float * float
type BBox =
{
minX: float
maxX: float
minY: float
maxY: float
center: float * float
type BBox = {
minX: float
maxX: float
minY: float
maxY: float
center: float * float
} with
static member empty = {
minX = Double.MaxValue
maxX = Double.MinValue
minY = Double.MaxValue
maxY = Double.MinValue
center = 0, 0
}
static member empty =
{
minX = Double.MaxValue
maxX = Double.MinValue
minY = Double.MaxValue
maxY = Double.MinValue
center = 0, 0
}
let unzip2D (array: _ [,]) =
let unzip2D (array: _[,]) =
let len1 = Array2D.length1 array
let len2 = Array2D.length2 array
let res1 = Array2D.zeroCreate len1 len2
@@ -34,7 +32,7 @@ let unzip2D (array: _ [,]) =
res2[i, j] <- y
res1, res2
let inline propTo3D (p: 'a [,]) =
let inline propTo3D (p: 'a[,]) =
let z = Array2D.length1 p
let n = Array2D.length2 p
let m = Array3D.zeroCreate 1 z n
@@ -43,7 +41,7 @@ let inline propTo3D (p: 'a [,]) =
m[0, i, j] <- p[i, j] |> single
m
let inline propTo2D (p: 'a []) =
let inline propTo2D (p: 'a[]) =
let n = p.Length
let m = Array2D.zeroCreate 1 n
for i = 0 to n - 1 do