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

Extracting and remapping nv to local indices.

Browse Source
This commit is contained in:
Isa Rosso
2025-08-15 12:02:27 +00:00
committed by Moritz Jörg
parent 55592f88c4
commit 6eabb35142
3 changed files with 105 additions and 42 deletions
Download Patch File Download Diff File Expand all files Collapse all files

9
src/Grid.fs
View File

@@ -17,6 +17,7 @@ type FvcomGrd = {
Node: int
Nele: int
nSiglays: int
nSiglevs: int
x: float[]
y: float[]
xc: float[]
@@ -26,7 +27,9 @@ type FvcomGrd = {
lonc: float[]
latc: float[]
siglay: float[,]
siglev: float[,]
h: float[]
nv: int[,]
//hc: float[]
}
@@ -39,6 +42,8 @@ let readFvcomGrid (ncfile: string) =
let node = dims["node"].Length
let nele = dims["nele"].Length
let nSiglays = dims["siglay"].Length
let nSiglevs = dims["siglev"].Length
let x = nc["x"].GetData () :?> float32[] |> Array.map float
let y = nc["y"].GetData () :?> float32[] |> Array.map float
let xc = nc["xc"].GetData () :?> float32[] |> Array.map float
@@ -51,6 +56,7 @@ let readFvcomGrid (ncfile: string) =
let h = nc["h"].GetData () :?> float32[] |> Array.map float
let siglay = nc["siglay"].GetData () :?> float32[,] |> Array2D.map float
let siglev = nc["siglev"].GetData () :?> float32[,] |> Array2D.map float
let nv = nc["nv"].GetData () :?> int[,]
let inds = [| for ind in 0..nele -> ind |]
//let hc = inds |> Array.map (fun i -> (h[nv[i, 0]] + h[nv[i, 1]] + h[nv[i, 2]]) / 3.)
@@ -61,6 +67,7 @@ let readFvcomGrid (ncfile: string) =
Node = node
Nele = nele
nSiglays = nSiglays
nSiglevs = nSiglevs
x = x
y = y
xc = xc
@@ -70,7 +77,9 @@ let readFvcomGrid (ncfile: string) =
lonc = lonc
latc = latc
siglay = siglay
siglev = siglev
h = h
nv = nv
//hc = hc
}

46
src/Main.fs
View File

@@ -219,8 +219,8 @@ let main argv =
//printfn "%A" nodeInds
let cellInds = findNearestCellInds grd stationListXY
//let cellInds = [|926882|]
createOutputFile outfile grd nodeInds cellInds |> ignore
exportToNc outfile fileArray nodeInds cellInds grd.nSiglays |> ignore
createOutputFile outfile grd nodeInds cellInds (None) |> ignore
exportToNc outfile fileArray nodeInds cellInds grd.nSiglays (None) |> ignore
elif args.Contains A then
let aArgs = args.GetResult A
@@ -244,12 +244,14 @@ let main argv =
let fileArray = createFileArray fileListCropped
//let fileArray = createFileArray fileListCropped
let grd = readFvcomGrid fileListCropped.Path[0]
let _, _, nodeInds = findNodesInside grd lonlims latlims
let _, _, nodeIndsPre = findNodesInside grd lonlims latlims
let _, _, cellInds = findCellsInside grd lonlims latlims
printfn "%d %d" nodeInds.Length cellInds.Length
printfn "%d %d" nodeIndsPre.Length cellInds.Length
let nodeInds, nv = findCellsNv fileListCropped.Path[0] grd cellInds
printfn "After checking nv, the length of cellInds is %d and of nodeInds is %d " cellInds.Length nodeInds.Length
//printfn "%A" fileListCropped
createOutputFile outfile grd nodeInds cellInds |> ignore
exportToNc outfile fileArray nodeInds cellInds grd.nSiglays |> ignore
createOutputFile outfile grd nodeInds cellInds (Some nv) |> ignore
exportToNc outfile fileArray nodeInds cellInds grd.nSiglays (Some nv) |> ignore
//exportTimestepToNc outfile fileListCropped nodeInds cellInds grd.nSiglays |> ignore
@@ -275,16 +277,19 @@ let main argv =
let fileArray = createFileArray fileListCropped
//let fileArray = createFileArray fileListCropped
let grd = readFvcomGrid fileListCropped.Path[0]
let _, _, nodeInds = findNodesInsideRadius grd center radius
let _, _, nodeIndsPre = findNodesInsideRadius grd center radius
let _, _, cellInds = findCellsInsideRadius grd center radius
printfn "%d %d" nodeIndsPre.Length cellInds.Length
let nodeInds, nv = findCellsNv fileListCropped.Path[0] grd cellInds
printfn "%d %d" nodeInds.Length cellInds.Length
//printfn "%A" fileListCropped
createOutputFile outfile grd nodeInds cellInds |> ignore
exportToNc outfile fileArray nodeInds cellInds grd.nSiglays |> ignore
//exportTimestepToNc outfile fileListCropped nodeInds cellInds grd.nSiglays |> ignore
elif args.Contains S then
createOutputFile outfile grd nodeInds cellInds (Some nv) |> ignore
exportToNc outfile fileArray nodeInds cellInds grd.nSiglays (Some nv) |> ignore
elif args.Contains S then
let sArgs = args.GetResult S
if sArgs.Contains ProjectionS then
printfn "%s" "Change projection"
projection <-
@@ -307,12 +312,12 @@ let main argv =
let grd = readFvcomGrid fileListCropped.Path[0]
let nodeInds = findNearestTransectNodes grd startPointXY endPointXY
let cellInds = findNearestTransectCells grd startPointXY endPointXY
createOutputFile outfile grd nodeInds cellInds |> ignore
exportToNc outfile fileArray nodeInds cellInds grd.nSiglays |> ignore
createOutputFile outfile grd nodeInds cellInds (None) |> ignore
exportToNc outfile fileArray nodeInds cellInds grd.nSiglays (None) |> ignore
elif args.Contains STAT then
let statArgs = args.GetResult STAT
if statArgs.Contains ProjectionStat then
printfn "%s" "Change projection"
projection <-
@@ -339,20 +344,7 @@ let main argv =
elif args.Contains GRD then
let grdArgs = args.GetResult GRD
if grdArgs.Contains ProjectionGrd then
printfn "%s" "Change projection"
projection <-
grdArgs.GetResult ProjectionGrd
|> ProjNet.FSharp.Projections.Projection.FromString
|> ProjNet.FSharp.Projections.Projection.ToCoordinateSystem
let inFile = grdArgs.GetResult InFile
let outfile = grdArgs.GetResult FileGrd
printfn "%s" outfile
writeGrid2Nc inFile outfile |> ignore
elif args.Contains GRD then
let grdArgs = args.GetResult GRD
if grdArgs.Contains ProjectionGrd then
printfn "%s" "Change projection"
projection <-

92
src/NetCDF.fs
View File

@@ -4,8 +4,10 @@ open System
//open Excavator.FileList
//open Excavator.Grid
open Excavator.FileList
open FSharpPlus
open FileList
open Grid
open Microsoft.Research.Science.Data.Imperative
open Microsoft.Research.Science.Data.NetCDF4
@@ -32,16 +34,64 @@ let readNcArray2D (nc: NetCDFDataSet) (var: string) (inds: int[]) =
let v = nc[var].GetData () :?> float32[,] |> Array2D.map float
readArray2D v inds
let createOutputFile (name: string) (grd: FvcomGrd) (nodeInds: int[]) (cellInds: int[]) =
let readArray2Dint (a: int[,]) (inds: int[]) =
let d = Array2D.zeroCreate (Array2D.length1 a) inds.Length
inds |> Array.iteri (fun i ind -> d[*, i] <- a[*, i])
d
let uniqueValuesSeq (arr: 'T[,]) =
arr
|> Seq.cast<'T> // flatten 2D array
|> Seq.distinct // filter unique values
|> toArray
let findCellsNv (fname: string) (grd: FvcomGrd) (cellInds: int[]) =
printfn "%s: checking cellInds from nv" fname
// FVCOM's starting index is 1-based, F#'s is 0-based
let nvFull =
let minNv = grd.nv |> Seq.cast<int> |> Seq.min
if minNv = 1 then
printfn "Detected 1-based indexing in grd.nv — converting to 0-based."
grd.nv |> Array2D.map (fun i -> i - 1)
else
printfn "grd.nv is already 0-based."
grd.nv
// Extract the nv subset
let nvSubset = Array2D.init 3 cellInds.Length (fun i j -> nvFull[i, cellInds[j]])
// Extract used node indices
let usedNodeInds = uniqueValuesSeq nvSubset
// Create global -> local mapping from used nodes only
let globalToLocal =
usedNodeInds
|> Array.mapi (fun localIdx globalIdx -> globalIdx, localIdx)
|> dict
// Remap nv using new mapping
let newNv =
Array2D.init 3 cellInds.Length (fun i j ->
let globalIdx = nvSubset[i, j]
globalToLocal[globalIdx])
let newNodeInds = usedNodeInds
newNodeInds, newNv
let createOutputFile (name: string) (grd: FvcomGrd) (nodeInds: int[]) (cellInds: int[]) (nv: int[,] option)=
let nNodes = nodeInds.Length
let nCells = cellInds.Length
let noNv = Array2D.zeroCreate 3 nCells
let nvCheck = defaultArg nv noNv
if IO.File.Exists name then
IO.File.Delete name
let saveFileNC = new NetCDFDataSet (name)
saveFileNC.CreateDimension ("siglay", grd.nSiglays)
saveFileNC.CreateDimension ("siglev", grd.nSiglevs)
saveFileNC.CreateDimension ("time", 0)
saveFileNC.CreateDimension ("three", 3)
saveFileNC.CreateDimension ("node", nNodes)
saveFileNC.CreateDimension ("nele", nCells)
@@ -65,7 +115,8 @@ let createOutputFile (name: string) (grd: FvcomGrd) (nodeInds: int[]) (cellInds:
//|> ignore
saveFileNC.AddVariable<float> ("siglay", (Array2D.zeroCreate<float> grd.nSiglays nNodes), [| "siglay"; "node" |])
|> ignore
saveFileNC.AddVariable<float> ("siglev", (Array2D.zeroCreate<float> grd.nSiglevs nNodes), [| "siglev"; "node" |])
|> ignore
saveFileNC.AddVariable<float> ("xc", (Array.zeroCreate<float> nCells), [| "nele" |])
|> ignore
saveFileNC.AddVariable<float> ("yc", (Array.zeroCreate<float> nCells), [| "nele" |])
@@ -80,10 +131,21 @@ let createOutputFile (name: string) (grd: FvcomGrd) (nodeInds: int[]) (cellInds:
saveFileNC.AddVariable<float> ("latc", (Array.zeroCreate<float> nCells), [| "nele" |])
|> ignore
saveFileNC.AddVariable<int> ("nodeidx", (Array.zeroCreate<int> nNodes), [| "node" |])
|> ignore
saveFileNC.AddVariable<int> ("cellidx", (Array.zeroCreate<int> nCells), [| "nele" |])
|> ignore
//saveFileNC.AddVariable<int> ("nodeidx", (Array.zeroCreate<int> nNodes), [| "node" |])
//|> ignore
//saveFileNC.AddVariable<int> ("cellidx", (Array.zeroCreate<int> nCells), [| "nele" |])
//|> ignore
// check if nv was provided
if not (obj.ReferenceEquals(nvCheck, noNv)) then
saveFileNC.AddVariable<int32> (
"nv",
(Array2D.zeroCreate<int32> 3 nCells),
[| "three"; "nele" |]
)
|> ignore
else
printfn "nv was not requested --> nv will not be saved in the netCDF file"
saveFileNC.AddVariable<float32> (
"temp",
@@ -130,6 +192,7 @@ let createOutputFile (name: string) (grd: FvcomGrd) (nodeInds: int[]) (cellInds:
saveFileNC["h"].PutData ((readArray grd.h nodeInds))
//saveFileNC["hc"].PutData ((readArray grd.hc cellInds))
saveFileNC["siglay"].PutData ((readArray2D grd.siglay nodeInds))
saveFileNC["siglev"].PutData ((readArray2D grd.siglev nodeInds))
saveFileNC.Commit ()
saveFileNC.Dispose ()
@@ -205,9 +268,6 @@ let readNcVar (fname: string) (var: string) (nSiglays: int) (nodeInds: int[]) (c
data
let readNcVarTimeStep
(fname: string)
(var: string)
@@ -237,12 +297,10 @@ let readNcVarTimeStep
readArray3D d grdInd |> Array3D.map float32
let exportToNc (outfile: string) (fl: fvcomFile[]) (nodeInd: int[]) (cellInd: int[]) (nSiglays: int) =
let exportToNc (outfile: string) (fl: fvcomFile[]) (nodeInd: int[]) (cellInd: int[]) (nSiglays: int) (nv: int[,] option)=
let vars = [| "temp"; "salinity"; "u"; "v"; "ww" |]
let ncout = NetCDFDataSet.Open outfile
let noNv = Array2D.zeroCreate 3 cellInd.Length
fl
|> Array.iter (fun f ->
let t, It, It2 = readNcFileTime f.Path (f.startIndex, f.stopIndex)
@@ -255,8 +313,12 @@ let exportToNc (outfile: string) (fl: fvcomFile[]) (nodeInd: int[]) (cellInd: in
ncout["Itime"].PutData ([| f.startArrayInd |], It)
ncout["Itime2"].PutData ([| f.startArrayInd |], It2)
ncout["zeta"].PutData ([| f.startArrayInd; 0 |], zeta)
ncout["nodeidx"].PutData ([| 0 |], nodeInd)
ncout["cellidx"].PutData ([| 0 |], cellInd)
let nvCheck = defaultArg nv noNv
if not (obj.ReferenceEquals(nvCheck, noNv)) then
printfn "nv was requested --> saving it"
printfn "%s nv" f.Path
ncout["nv"].PutData ([| 0; 0 |], nvCheck)
ncout.Commit ())
//ncout.Commit()