Oceanbox.FvcomKit/xtest/Arome.fs

module Arome

open System

open Microsoft.Research.Science.Data

open Xunit
open FsUnit.Xunit
open FsUnit.CustomMatchers
open Serilog
open ProjNet.FSharp

open Oceanbox.FvcomKit

let logger =
    LoggerConfiguration()
        .MinimumLevel.Is(Events.LogEventLevel.Verbose)
        .WriteTo.Console()
        .CreateLogger()

do Log.Logger <- logger

[<Literal>]
let private path = "/data/archives/Arome/meps_det_sfc_20220102T00Z.nc"

let trans = makeTransform CoordSys.WGS84 (CoordSys.LCCMet ())

// Use Netcdf to open dataset
let private openDataset (path: string) : DataSet =
    let sw = Diagnostics.Stopwatch.StartNew()
    let uri = NetCDF4.NetCDFUri()
    do uri.FileName <- path
    do uri.OpenMode <- ResourceOpenMode.ReadOnly
    do uri.Deflate <- NetCDF4.DeflateLevel.Off
    let ds = NetCDF4.NetCDFDataSet.Open uri
    do Log.Debug $"openDataSet: {path} completed in {sw.ElapsedMilliseconds}ms"

    ds

let private aromeGrid = {
    Arome.SquareGrid.empty with
        dimensions = 949, 1069
        BBox = {
            minX = -18.12241427
            maxX = 54.24126163
            minY = 49.7653854
            maxY = 75.22869642
            center = 27.12063081, 37.61434821
        }
        squareSize = 2.444065489
}

[<Fact>]
let ``Create grid``() =
    use ds = openDataset path
    let res = Arome.getGrid ds

    //  "Should be able to create square grid from /data/archives/Arome/meps_det_sfc_20220102T00Z.nc"
    res |> should be (ofCase <@ Result<Arome.SquareGrid, string>.Ok @>)

[<Fact>]
let ``test grid regularity``() =
    let trans = makeTransform CoordSys.WGS84 (CoordSys.LCCMet ())
    use ds = openDataset path
    let dimensions = ds.Dimensions["x"].Length, ds.Dimensions["y"].Length
    let longs : float array2d = ds["longitude"].GetData () :?> float[,]
    let lats : float array2d = ds["latitude"].GetData () :?> float[,]
    // NOTE: The netcdf file dimensions are defined as (y, x)
    let width = Array2D.length2 longs
    let height = Array2D.length1 lats
    let points : Arome.Points array =
        let result = Array.create (width * height) Arome.Points.Zero
        for i in 0 .. height - 1 do
            for j in 0 .. width - 1 do
                let lat =  lats[i, j]
                let lon = longs[i, j]
                let p = lon, lat
                // NOTE(simkir): Convert to lambert and undo all I've done... :(
                let x, y = trans.project p
                result[i * width + j] <- Arome.Points.OfTuple (single x, single y)
        result

    let mutable acc = 0.0f
    for i in 0 .. height - 1 do
        for j in 0 .. width - 2 do
            let p0 = points[i * width + j]
            let p1 = points[i * width + j + 1]
            let dx = if p0.x < p1.x then p1.x - p0.x else p0.x - p1.x

            acc <- acc + dx

    let count = width * height
    let avgDiffX = acc / single count
    Log.Debug("Avg. X distance: {Avg}m", avgDiffX)

    acc <- 0.0f
    for i in 0 .. height - 2 do
        for j in 0 .. width - 1 do
            let p0 = points[i * width + j]
            let p1 = points[(i + 1) * width + j]
            let dy = if p0.y < p1.y then p1.y - p0.y else p0.y - p1.y

            acc <- acc + dy
    let avgDiffY = acc / single count
    Log.Debug("Avg. Y distance: {Avg}m", avgDiffY)

[<Fact>]
let ``point within first rect``() =
    use ds = openDataset path
    let res = Arome.getGrid ds

    match res with
    | Ok grid ->
        let width, _height = grid.dimensions
        let p0 = grid.points[0 * width + 0]
        let p1 = grid.points[1 * width + 1]
        let rect = Arome.Rect.OfPoints p0 p1
        // { x = -1064984.25f y = -1353489.375f },  { x = -1062474.625f y = -1350985.875f }
        let p2 = Arome.Points.OfTuple (-1063000.25f, -1353486.0f)
        let within = Arome.Rect.pointWithin rect p2

        within |> should equal true
    | Error _ ->
        failwith "Should create grid"

[<Fact>]
let ``Read uvs for index (1, 1)``() =
    use ds = openDataset path
    let res = Arome.getGrid ds
    match res with
    | Ok grid ->
        let idx =
            let lat = 50.35
            let lon =  0.304
            let p = trans.project((lon, lat))
            Arome.tryFindIndex grid p

        match idx with
        | Some (i, j) ->
            let u, v = Arome.readUV ds 0 i j
            (u, v) |> should equal (8.530239f, 6.659096f) // "Find uv (8.530239, 6.659096) at index (1, 1)"
        | None ->
            failwith "Should find idx"
    | Error _ ->
        failwith "Should find grid"

[<Fact>]
let ``Try search index (1, 1069) based on map coords``() =
    let lat = 72.800
    let lon = -17.850
    let p = trans.project((lon, lat))
    let res =
        use ds = openDataset path
        Arome.getGrid ds

    match res with
    | Ok grid ->
        let idx = Arome.tryFindIndex grid p

        idx |> should equal (Some (4, 1067)) // "Should find idx (1, 1069)"
    | Error _ ->
        failwith "Should find grid"

[<Fact>]
let ``Read uvs for pos (60.3562, 5.2178)``() =
    let lat = 60.3562
    let lon = 5.2178
    let p = trans.project((lon, lat))
    use ds = openDataset path
    let res = Arome.getGrid ds
    match res with
    | Ok grid ->
        let idx = Arome.tryFindIndex grid p
        match idx with
        | Some (i, j) ->
            let u, v = Arome.readUV ds 0 i j
            // "Find uv (2.367153168, -0.3955917358) at pos (60.3562, 5.2178) almost in Bergen"
            (u, v) |> should equal (-0.3750343323f, 7.744056702f)
        | None ->
            failwith "Should find idx"
    | Error _ ->
        failwith "Should find grid"

[<Fact>]
let ``Read uvs for t 23 x 43 y 144``() =
    use ds = openDataset path
    let u, v = Arome.readUV ds 23 43 144
    // "Find uv (2.367153168, -0.3955917358) at pos (60.3562, 5.2178) almost in Bergen"
    (u, v) |> should equal (11.77926254f, 5.764093399f)