devel: update tilt settings

fix: use consistent lonlat projection on all cached grids
Set cwd for atlantis/just run-client
2025-12-17 17:32:17 +01:00 · 2025-12-17 17:30:12 +01:00 · 2025-12-17 15:33:43 +01:00
16 changed files with 437 additions and 428 deletions
--- a/src/Atlantis/justfile
+++ b/src/Atlantis/justfile
@@ -56,8 +56,9 @@ run-server:
    dotnet watch run {{server_path}}

 # Run client only in watch mode
+[working-directory: 'src/Client']
 run-client:
-    fable watch --cwd {{client_path}} -e .jsx -o build --run {{vite}}
+    fable watch -e .jsx -o build --run {{vite}}

 # Format code with Fantomas
 format:
--- a/src/Atlantis/src/Client/Lib/Search.fs
+++ b/src/Atlantis/src/Client/Lib/Search.fs
@@ -9,8 +9,9 @@ let tryGetNearestNodeAndElement aid (x, y) : Async<Option<NodeIdx * ElemIdx>> =
    async {
        let dataSvc = Remoting.getDataUrl ()
        let api = Remoting.proximityApi dataSvc
-        let! nodes = api.GetNearestNodes (aid, [| x, y |])
-        let! elems = api.GetNearestElements (aid, [| x, y |])
+        let lonlat = (x, y) |> Utils.toWgs84'
+        let! nodes = api.GetNearestNodes (aid, [| lonlat |])
+        let! elems = api.GetNearestElements (aid, [| lonlat |])
        let firstNode = nodes |> Array.tryHead |> Option.flatten
        let firstElem = elems |> Array.tryHead |> Option.flatten
        let result =
@@ -24,8 +25,9 @@ let tryGetNearestNodeAndElement aid (x, y) : Async<Option<NodeIdx * ElemIdx>> =
 let tryGetNearestNodeAndElementRes dataSvc aid (x, y) : Async<Result<NodeIdx * ElemIdx, exn>> =
    asyncResult {
        let api = Remoting.proximityApi dataSvc
-        let! nodes = api.GetNearestNodes (aid, [| x, y |]) |> Remoting.tryCatch
-        let! elems = api.GetNearestElements (aid, [| x, y |]) |> Remoting.tryCatch
+        let lonlat = (x, y) |> Utils.toWgs84'
+        let! nodes = api.GetNearestNodes (aid, [| lonlat |])
+        let! elems = api.GetNearestElements (aid, [| lonlat |])
        let firstNode = nodes |> Array.tryHead |> Option.flatten
        let firstElem = elems |> Array.tryHead |> Option.flatten
        let result =
--- a/src/Atlantis/src/Client/Mapster/Mapster.fs
+++ b/src/Atlantis/src/Client/Mapster/Mapster.fs
@@ -191,7 +191,7 @@ let initializeArchive model =
                plainGrid.Grid |> async.Return
            | None ->
                console.debug $"Grid not found in db, fetching from server {model.archive.sha}"
-                dataApi.Archive.GetPlainGrid aid
+                dataApi.Archive.GetPlainGridEpsg3857 aid

        if plainGridFromDB.IsNone then
            console.log $"Storing grid in db: {model.archive.sha}"
--- a/src/Atlantis/src/Client/Mapster/Plots/ContourPlots.fs
+++ b/src/Atlantis/src/Client/Mapster/Plots/ContourPlots.fs
@@ -189,11 +189,10 @@ let contourPlots archive (frame: int) (probePoint: Probing) (pickLine: (float *
                console.debug("[ContourPlots] Created sampling line %o", pts)

                pts
+                |> Array.map Utils.toWgs84'
                |> fetchContourData archive.id frame probePoint.prop
                |> Async.StartAsPromise
-                |> Promise.iter (fun d ->
-                    createTraces d
-                )
+                |> Promise.iter createTraces
    )

    Hook.useEffectOnChange (
@@ -209,11 +208,10 @@ let contourPlots archive (frame: int) (probePoint: Probing) (pickLine: (float *
                console.debug("[ContourPlots] Created sampling line %o", pts)

                pts
+                |> Array.map Utils.toWgs84'
                |> fetchContourData archive.id frame newProp
                |> Async.StartAsPromise
-                |> Promise.iter (fun d ->
-                    createTraces d
-                )
+                |> Promise.iter createTraces
    )

    let spinner () =
--- a/src/Atlantis/src/Client/Mapster/Postdrift.fs
+++ b/src/Atlantis/src/Client/Mapster/Postdrift.fs
@@ -211,7 +211,7 @@ let fetchDriftersGrid m =
            |> Option.map _.Archive.archiveId
            |> Option.defaultValue Guid.Empty
        let fvcomApi = FvcomApi m.dataSvc
-        let! g = fvcomApi.Archive.GetPlainGrid aid
+        let! g = fvcomApi.Archive.GetPlainGridEpsg3857 aid

        let data: MapData = {
            Grid = g
--- a/src/Atlantis/src/Server/packages.lock.json
+++ b/src/Atlantis/src/Server/packages.lock.json
@@ -1545,7 +1545,7 @@
          "Hipster.Api": "[1.0.1, )",
          "Oceanbox.DataAgent": "[7.3.0, )",
          "Oceanbox.DataAgent.Api": "[7.2.1, )",
-          "Oceanbox.ServerPack": "[1.35.2, )",
+          "Oceanbox.ServerPack": "[1.36.0, )",
          "Petimeter.Api": "[1.0.0, )"
        }
      },
--- a/src/Interfaces/Sorcerer/Api.fs
+++ b/src/Interfaces/Sorcerer/Api.fs
@@ -15,7 +15,7 @@ module Api =
        type Archive = {
            GetTime: Guid -> FrameIdx -> Async<DateTime option>
            GetNLayers: Guid -> Async<int>
-            GetPlainGrid: Guid -> Async<PlainGrid>
+            GetPlainGridEpsg3857: Guid -> Async<PlainGrid>
            GetNumFrames: Guid -> Async<int>
            GetTimeStep: Guid -> Async<int>
            EquipArchiveCache: Guid -> Async<string option>
@@ -136,9 +136,10 @@ module Api =
        let routeBuilder = routeBuilder "prox"

        type Prox = {
-            GetNearestElements: Guid * (float * float)[] -> Async<ElemIdx option[]>
-            GetNearestNodes: Guid * (float * float)[] -> Async<NodeIdx option[]>
-            GetNodesInCircle: Guid * float * float * float -> Async<NodeIdx []>
+            GetNearestElements: Guid * (float * float) array -> Async<ElemIdx option array>
+            GetNearestNodes: Guid * (float * float) array -> Async<NodeIdx option array>
+            GetNodesInCircle: Guid * float * float * float -> Async<NodeIdx array>
+            GetNodesInPolygon: Guid * (float * float) array -> Async<NodeIdx array>
            GetGridDistanceCircle: Guid * float * float * float -> Async<(float*float) array>
            GetShortestNodePath: Guid * (float * float) * (float * float) -> Async<GridPath option>
            GetShortestElementPath: Guid * (float * float) * (float * float) -> Async<GridPath option>
--- a/src/Sorcerer/Tiltfile
+++ b/src/Sorcerer/Tiltfile
@@ -66,7 +66,7 @@ local_resource(

 local_resource(
    'build-server',
-    cmd='just bundle-debug-server',
+    cmd='just bundle-debug',
    deps=[
        './src/Server',
        './src/Shared'
--- a/src/Sorcerer/src/Server/Api.fs
+++ b/src/Sorcerer/src/Server/Api.fs
@@ -286,13 +286,13 @@ module Fvcom =
                    use _ = observer.trace ("GetTime", "{@log_data}", logData)
                    return Fvcom.archiveTime aid n
                }
-            GetPlainGrid =
+            GetPlainGridEpsg3857 =
                fun aid ->
                    async {
                        let logData = {| aid = aid |}
-                        use _ = observer.trace ("GetPlainGrid", "{@log_data}", logData)
+                        use _ = observer.trace ("GetPlainGridEpsg3857", "{@log_data}", logData)
                        return
-                            Fvcom.retrievePlainGrid aid
+                            Fvcom.retrievePlainGridEpsg3857 aid
                            |> Result.defaultValue PlainGrid.empty
                    }
            EquipArchiveCache = fun aid ->
@@ -579,8 +579,9 @@ module Proximity =
        let uid = ctx.User.Identity.Name
        let observer = ctx.GetService<ObserverFactory>().Create(uid, tag="Proximity.prox")
        {
-            GetNearestElements = fun (aid, coords) -> async { return Fvcom.elems aid coords }
-            GetNearestNodes = fun (aid, coords) -> async { return Fvcom.nodes aid coords }
+            GetNearestElements = fun (aid, coords) -> async { return Fvcom.nearestElements aid coords }
+            GetNearestNodes = fun (aid, coords) -> async { return Fvcom.nearestNodes aid coords }
+            GetNodesInPolygon = fun (aid, poly) -> async { return Fvcom.Proximity.nodesInPolygon aid poly }
            GetNodesInCircle = fun (aid, x, y, r) -> async { return Fvcom.Proximity.nodesInCircle aid x y r }
            GetGridDistanceCircle = fun (aid, x, y, r) -> async { return Fvcom.Proximity.gridDistanceCircle aid x y r }
            GetShortestNodePath = fun (aid, p1, p2) -> async { return Fvcom.Proximity.shortestNodePath aid p1 p2 }
@@ -732,14 +733,14 @@ module Ocean =
            |> Ocean.Property.Speed
        | _ -> Ocean.Property.Void "request failed"

-    let getGridIndex (prop: OceanProp) (aid: Guid) ((x, y): float * float) =
-        let nIdx = Fvcom.nodeLonLat aid x y
+    let getGridIndex (prop: OceanProp) (aid: Guid) (coord: float * float) =
+        let nIdx = Fvcom.nearestNode aid coord
        let idx =
            match prop with
            | OceanProp.Temperature
            | OceanProp.Salinity
            | OceanProp.Density -> nIdx
-            | _ -> Fvcom.elemLonLat aid x y
+            | _ -> Fvcom.nearestElement aid coord
        match nIdx, idx with
        | Some n, Some i -> Some (n, i)
        | _ -> None
--- a/src/Sorcerer/src/Server/Crop.fs
+++ b/src/Sorcerer/src/Server/Crop.fs
@@ -43,9 +43,10 @@ let calcVerticalInterpolationCoefs (depths: single[]) (z: single) =
 // coord0 is in WebMercator
 // dx, dy is in actual (not latitude scaled) meters
 let cropGrid (aid: Guid) (coord0: float * float) (dx: float, dy: float) (nx, ny) =
+    failwith "needs fix: projection"
    let contribs = Array2D.init<int option> nx ny (fun _ _ -> None)
    let x0, y0 = coord0
-    match Fvcom.retrieveExtendedGrid true aid with
+    match Fvcom.retrieveExtendedGrid aid with
    | Ok grid ->
        for j = 0 to ny - 1 do
            for i = 0 to nx - 1 do
@@ -80,9 +81,10 @@ let cropf (f: int -> Result<single[], string>) (aid: Guid) (coord0: float * floa
    result

 let cropf' (f: int -> Result<single[], string>) (aid: Guid) (coord0: float * float) (dx: float, dy: float) (nx, ny) (z: float) =
+    failwith "needs fix: projection"
    let result = Array.create (nx * ny) -99999f
    let contribs = cropGrid aid coord0 (dx, dy) (nx, ny)
-    match Fvcom.retrieveExtendedGrid true aid with
+    match Fvcom.retrieveExtendedGrid aid with
    | Ok grid ->
        let distinctMap =
            contribs
@@ -112,9 +114,10 @@ let cropf' (f: int -> Result<single[], string>) (aid: Guid) (coord0: float * flo
    result

 let cropf'' (f: int -> Result<(single * single)[], string>) (aid: Guid) (coord0: float * float) (dx: float, dy: float) (nx, ny) (z: float) =
+    failwith "needs fix: projection"
    let result = Array.zeroCreate (nx * ny)
    let contribs = cropGrid aid coord0 (dx, dy) (nx, ny)
-    match Fvcom.retrieveExtendedGrid true aid with
+    match Fvcom.retrieveExtendedGrid aid with
    | Ok grid ->
        for j = 0 to ny - 1 do
            for i = 0 to nx - 1 do
@@ -263,12 +266,13 @@ module Avg =
    // coord0 is in WebMercator
    // dx, dy is in actual (not latitude scaled) meters
    let cropGrid (aid: Guid) (coord0: float * float) (dx: float, dy: float) (nx, ny) =
+        failwith "needsfix: projection"
        let contribs = Array2D.init nx ny (fun _ _ -> [||])
        let x0, y0 = coord0
        // let _, lat = Utils.webMercatorToWGS84.project coord0
        // let metric = Utils.mercatorScaleFactor lat
        // let dx, dy = dx * metric, dy * metric
-        match Fvcom.retrieveExtendedGrid true aid with
+        match Fvcom.retrieveExtendedGrid aid with
        | Ok grid ->
            for j = 0 to ny - 1 do
                for i = 0 to nx - 1 do
@@ -309,9 +313,10 @@ module Avg =
        result

    let cropf' (f: Guid -> int -> int -> Result<single[], string>) (aid: Guid) t (coord0: float * float) (dx: float, dy: float) (nx, ny) (z: float) =
+        failwith "needsfix: projection"
        let result = Array.zeroCreate (nx * ny)
        let contribs = cropGrid aid coord0 (dx, dy) (nx, ny)
-        match Fvcom.retrieveExtendedGrid true aid with
+        match Fvcom.retrieveExtendedGrid aid with
        | Ok grid ->
            for j = 0 to ny - 1 do
                for i = 0 to nx - 1 do
@@ -336,9 +341,10 @@ module Avg =
        result

    let cropf'' (f: Guid -> int -> int -> Result<(single * single)[], string>) (aid: Guid) t (coord0: float * float) (dx: float, dy: float) (nx, ny) (z: float) =
+        failwith "needsfix: projection"
        let result = Array.zeroCreate (nx * ny)
        let contribs = cropGrid aid coord0 (dx, dy) (nx, ny)
-        match Fvcom.retrieveExtendedGrid true aid with
+        match Fvcom.retrieveExtendedGrid aid with
        | Ok grid ->
            for j = 0 to ny - 1 do
                for i = 0 to nx - 1 do
--- a/src/Sorcerer/src/Server/Downloads.fs
+++ b/src/Sorcerer/src/Server/Downloads.fs
@@ -33,10 +33,9 @@ let allPeriods = [| Year; for n = 0 to 11 do Month n |]
 let allMetrics = [| Mean; Std; Var; Q05; Q25; Q50; Q75; Q95; Q99 |]

 let tryGetNearestNodeAndElement aid (coord: float * float) =
-    let x, y = Utils.lngLatToWebMercator.project coord
    monad {
-        let! node = Fvcom.node aid x y
-        let! elem = Fvcom.elem aid x y
+        let! node = Fvcom.nearestNode aid coord
+        let! elem = Fvcom.nearestElement aid coord
        return node, elem
    }

@@ -384,8 +383,8 @@ module Contour =
        let ixs =
            match prop with
            | Prop.Temp
-            | Prop.Salt -> Fvcom.nodes aid coords
-            | Prop.Speed -> Fvcom.elems aid coords
+            | Prop.Salt -> Fvcom.nearestNodes aid coords
+            | Prop.Speed -> Fvcom.nearestElements aid coords
            | _ -> [||]
        let idx = // remove duplicates (same node/element)
            ixs
--- a/src/Sorcerer/src/Server/Drifters.fs
+++ b/src/Sorcerer/src/Server/Drifters.fs
@@ -3,6 +3,7 @@ module Drifters

 open System
 open FSharpPlus.Control
+open ProjNet.FSharp
 open Microsoft.Research.Science.Data.NetCDF4
 open Microsoft.Research.Science.Data
 open Serilog
@@ -576,37 +577,27 @@ let getReleaseSites aid =
 // CONTOURS

 let distanceContour aid npoints distance =
-    let mercatorScaleFactor (lat: float) = 1.0 / (lat * Math.PI / 180.0 |> Math.Cos)
-
-    let toWebMercator (x, y, r) =
-        let x', y' = (float x, float y) |> Utils.lngLatToWebMercator.project
-        let r' = (float r + distance) * mercatorScaleFactor (float y)
-        x', y', r'
-
-    let sampleCircle np (x, y, r)  =
-        [| 0..np |]
-        |> Array.map (fun i ->
-            let theta = 2.0 * Math.PI * float i / float np
-            let x' = x + r * cos theta
-            let y' = y + r * sin theta
-            (x', y'))
+    let geodesicCircle (x, y, r) =
+        [| 0..npoints |]
+        |> Array.map (fun i -> 2.0 * Math.PI * float i / float npoints)
+        |> Array.map (Utils.destinationPointSphere (x, y) r)

    aid
    |> getReleaseSites
    |> Result.map (fun sites ->
        sites
        |> Array.distinctBy (fun (x, y, _) -> $"%.6f{x} %.6f{y}")
-        |> Array.Parallel.map (toWebMercator >> sampleCircle npoints)
+        |> Array.map (fun (x, y, r) -> float x, float y, float r + distance)
+        |> Array.Parallel.map geodesicCircle
        |> Array.concat
        |> Contours.convexHull
-        |> Array.Parallel.map Utils.webMercatorToWGS84.project
        |> fun x -> [| x |])

 let sedimentationContour aid t field sediment ckind value =
    let name = $"{field.ToString()}_{sediment.ToString()}"

    let grid =
-        match retrieveExtendedGrid true aid with
+        match retrieveExtendedGrid aid with
        | Ok g -> g.ToGrid ()
        | _ -> failwith "error reading grid"

@@ -633,7 +624,6 @@ let sedimentationContour aid t field sediment ckind value =
            Node2D.readField ds t name
            |> Array.map float
            |> calcContour
-            |> Array.map (Array.Parallel.map Utils.webMercatorToWGS84.project)

    dataAgent.eval (fCont, aid, t)

@@ -645,7 +635,7 @@ let contour2D aid t filter ckind value =
        | _ -> failwith "error reading field names"

    let grid =
-        match retrieveExtendedGrid true aid with
+        match retrieveExtendedGrid aid with
        | Ok g -> g.ToGrid ()
        | _ -> failwith "error reading grid"

@@ -672,7 +662,6 @@ let contour2D aid t filter ckind value =
            Node2D.accumulateFields ds t names (+) 0.f
            |> Array.map float
            |> calcContour
-            |> Array.map (Array.Parallel.map Utils.webMercatorToWGS84.project)

    dataAgent.eval (fCont, aid, t)

@@ -684,7 +673,7 @@ let contour3D aid t filter ckind value =
        | _ -> failwith "error reading field names"

    let grid =
-        match retrieveExtendedGrid true aid with
+        match retrieveExtendedGrid aid with
        | Ok g -> g.ToGrid ()
        | _ -> failwith "error reading grid"

@@ -711,7 +700,6 @@ let contour3D aid t filter ckind value =
            Node3D.accumulateFields ds t names filter.layers (+) 0.f
            |> Array.map float
            |> calcContour
-            |> Array.map (Array.Parallel.map Utils.webMercatorToWGS84.project)

    dataAgent.eval (fCont, aid, t)

@@ -736,7 +724,7 @@ let calcConnectionMatrix2D aid t filter =
            | _ -> failwith "error reading release sites"

    let grid =
-        match retrieveExtendedGrid false aid with
+        match retrieveExtendedGrid aid with
        | Ok g -> g
        | _ -> failwith "error reading grid"

--- a/src/Sorcerer/src/Server/Fvcom.fs
+++ b/src/Sorcerer/src/Server/Fvcom.fs
@@ -9,14 +9,12 @@ open Oceanbox.FvcomKit.Grid
 open Serilog
 open MessagePack
 open Sorcerer.Types
-open ProjNet.FSharp
 open Settings
 open Oceanbox.FvcomKit
 open Oceanbox.DataAgent

-open Utils
-
 let private cache = CacheAgent.CacheAgent<ExtendedGrid>(60 * 10)
+let private cacheExtGrid aid (eg: ExtendedGrid) = cache.add(aid, eg); eg

 module Auxiliary =
    let toPlainGrid (grid: IGrid) : PlainGrid = {
@@ -34,12 +32,102 @@ module Auxiliary =
        uv
        |> Array.Parallel.map calcSpeed

+    /// Compute density of seawater from salinity, temperature, and pressure
+    let calcDensity (temp: float) (salt: float) (pres: float) =
+        // From Bjørn Aadlandsvik https://github.com/bjornaa/seawater/blob/master/seawater/density.py
+
+        /// Density of seawater at zero pressure
+        let dens0 (salt: float) (temp: float) =
+            // --- Define constants ---
+            let a0 = 999.842594
+            let a1 = 6.793952e-2
+            let a2 = -9.095290e-3
+            let a3 = 1.001685e-4
+            let a4 = -1.120083e-6
+            let a5 = 6.536332e-9
+
+            let b0 = 8.24493e-1
+            let b1 = -4.0899e-3
+            let b2 = 7.6438e-5
+            let b3 = -8.2467e-7
+            let b4 = 5.3875e-9
+
+            let c0 = -5.72466e-3
+            let c1 = 1.0227e-4
+            let c2 = -1.6546e-6
+
+            let d0 = 4.8314e-4
+
+            // --- Computations ---
+            // Density of pure water
+            let SMOW = a0 + (a1 + (a2 + (a3 + (a4 + a5 * temp) * temp) * temp) * temp) * temp
+
+            // More temperature polynomials
+            let RB = b0 + (b1 + (b2 + (b3 + b4 * temp) * temp) * temp) * temp
+            let RC = c0 + (c1 + c2 * temp) * temp
+
+            SMOW + RB * salt + RC * (salt ** 1.5) + d0 * salt * salt
+
+        /// Secant bulk modulus
+        let seck salt temp pres =
+            // --- Pure water terms ---
+            let h0 = 3.239908
+            let h1 = 1.43713e-3
+            let h2 = 1.16092e-4
+            let h3 = -5.77905e-7
+            let AW = h0 + (h1 + (h2 + h3 * temp) * temp) * temp
+
+            let k0 = 8.50935e-5
+            let k1 = -6.12293e-6
+            let k2 = 5.2787e-8
+            let BW = k0 + (k1 + k2 * temp) * temp
+
+            let e0 = 19652.21
+            let e1 = 148.4206
+            let e2 = -2.327105
+            let e3 = 1.360477e-2
+            let e4 = -5.155288e-5
+            let KW = e0 + (e1 + (e2 + (e3 + e4 * temp) * temp) * temp) * temp
+
+            // --- Seawater, P = 0 ---
+            let SR = salt ** 0.5
+
+            let i0 = 2.2838e-3
+            let i1 = -1.0981e-5
+            let i2 = -1.6078e-6
+            let j0 = 1.91075e-4
+            let A = AW + (i0 + (i1 + i2 * temp) * temp + j0 * SR) * salt
+
+            let f0 = 54.6746
+            let f1 = -0.603459
+            let f2 = 1.09987e-2
+            let f3 = -6.1670e-5
+            let g0 = 7.944e-2
+            let g1 = 1.6483e-2
+            let g2 = -5.3009e-4
+            let K0 =
+                KW
+                + (f0 + (f1 + (f2 + f3 * temp) * temp) * temp + (g0 + (g1 + g2 * temp) * temp) * SR)
+                  * salt
+
+            // --- General expression ---
+            let m0 = -9.9348e-7
+            let m1 = 2.0816e-8
+            let m2 = 9.1697e-10
+            let B = BW + (m0 + (m1 + m2 * temp) * temp) * salt
+
+            K0 + (A + B * pres) * pres
+
+        // Convert to bar
+        let pres = 0.1 * pres
+        (dens0 salt temp) / (1.0 - pres / seck salt temp pres)
+
    let toBBox (bbox: Types.BBox) : BBox = {
        minX = single bbox.minX
        minY = single bbox.minY
        maxX = single bbox.maxX
        maxY = single bbox.maxY
-        center = bimap single bbox.center
+        center = Utils.bimap single bbox.center
    }

    // TODO(SimenLK): Decide if this is Sorcerer's responsibility
@@ -74,18 +162,47 @@ module Auxiliary =
            Log.Error $"{e}"
            Error e.Message

-    let readUVRange aid t l e0 en =
-        let f ds t = Fvcom.readUVRange ds t l e0 en
-        dataAgent.eval (f, aid, t) |> Result.get
+    // Grid files assumed in LonLat
+    let retrieveGrid aid : Result<Grid, string>  =
+        let s = Diagnostics.Stopwatch()
+        s.Start()
+        let dto =
+            try
+                archiveAgent.tryGetArchive aid
+                |> Option.defaultWith (fun () -> failwith $"retrieveGrid: could not load archive {aid}")
+            with exn ->
+                failwith $"%A{exn}"
+        Log.Debug $"retrieveGid: dto: {s.ElapsedMilliseconds}ms"

-    let readUV aid t l e =
-        let f ds t = Fvcom.Siglay.readUv ds e t
-        dataAgent.eval (f, aid, t)
-        |> Result.get
-        |> Array.item l
-
-    let getGrid ds = Fvcom.getGrid(ds).ToGrid()
-    let getGridLonLat ds = Fvcom.getGridLonLat(ds).ToGrid()
+        s.Restart ()
+        let toPath f = IO.Path.Join [| dto.files.basePath; f |]
+        let f0 = toPath "grid.bin" // must be in lon, lat
+        let f1 = toPath "grid.dat" // must be in lon, lat
+        if IO.File.Exists f0 then
+            Log.Debug $"Reading grid from {f0}"
+            let bytes = IO.File.ReadAllBytes f0
+            let g = MessagePackSerializer.Deserialize<BinGrid> bytes
+            g.toGrid() |> Ok
+        elif IO.File.Exists f1 then
+            Log.Debug $"Reading grid from {f1}"
+            match readGrdFile f1 with
+            | Some g -> Ok g
+            | None -> Error $"Read grid.dat failed: {aid}"
+        else
+            match dto.props.archiveType with
+            | ArchiveType.Fvcom _ ->
+                Log.Debug $"Reading lonlat grid from NetCDF archive: {aid}"
+                let f ds _ = Fvcom.getGridLonLat(ds).ToGrid()
+                dataAgent.eval (f, aid)
+            | ArchiveType.Drifters _ ->
+                if List.contains dto.props.projection [ "WGS84"; "LonLat" ] then
+                    Log.Debug $"Reading grid from NetCDF archive: {aid}"
+                    let f ds _ = Fvcom.getGrid(ds).ToGrid()
+                    dataAgent.eval (f, aid)
+                else
+                    Error $"Drifter grid in wrong projection : {dto.props.projection}"
+            | _ ->
+                Error $"Not able to fetch grid in LonLat projection : {aid}"

    let reviewNix sha (eg: ExtendedGrid) =
        let nix = $"{appsettings.cacheDir}/{sha}-nix.dat"
@@ -137,72 +254,37 @@ module Auxiliary =

 open Auxiliary

-// Returns grid in WGS84 (LonLat) projection
-let private retrieveGrid aid : Result<ExtendedGrid, string>  =
+
+// Grid in LonLat -> PlainGrid in EPSG:3857
+// For use in Atlantis map, more convenient to do projection to WebMercator here
+let retrievePlainGridEpsg3857 (aid: Guid) =
    let s = Diagnostics.Stopwatch()
    s.Start()
-    let dto =
-        try
-            archiveAgent.tryGetArchive aid
-            |> Option.defaultWith (fun () -> failwith $"retrieveGrid: could not load archive {aid}")
-        with exn ->
-            failwith $"%A{exn}"
-    Log.Debug $"retrieveGid: dto: {s.ElapsedMilliseconds}ms"
-    let toLonLat (g: Grid) =
-        Transformations.stringToTransformation dto.props.projection
-        |> Option.defaultWith (fun () ->
-            failwith $"retrieveGrid: failed to get projecton coordsys for {dto.props.projection}: {aid}")
-        |> fun coordsys -> toLonLat coordsys g
+    let grid =
+        match cache.tryGet aid with
+        | Some(eg: ExtendedGrid) ->
+            Log.Debug $"plainGrid: in cache {aid}"
+            eg |> Ok
+        | None ->
+            Log.Debug $"plainGrid: fetch {aid}"
+            retrieveGrid aid
+            |> Result.map (ExtendedGrid >> cacheExtGrid aid)
+    Log.Debug $"retrievePlainGrid(): {s.ElapsedMilliseconds}ms"

-    s.Restart ()
-    let toPath f = IO.Path.Join [| dto.files.basePath; f |]
-    let f0 = toPath "grid.bin" // must be in lon, lat
-    let f1 = toPath "grid.dat" // must be in lon, lat
-    if IO.File.Exists f0 then
-        Log.Debug $"Reading grid from {f0}"
-        let bytes = IO.File.ReadAllBytes f0
-        let g = MessagePackSerializer.Deserialize<BinGrid> bytes
-        let eg = g.toGrid () |> ExtendedGrid
-        eg.initHash g.hash
-        Ok eg
-    elif IO.File.Exists f1 then
-        Log.Debug $"Reading grid from {f1}"
-        match readGrdFile f1 with
-        | Some g -> Ok (ExtendedGrid g)
-        | None -> Error $"Read grid.dat failed: {aid}"
-    else
-        Log.Debug $"Reading grid from NetCDF archive: {aid}"
-        let f ds _ =
-            if dto.props.projection = "WGS84" then
-                getGridLonLat ds
-            else
-                getGrid ds |> toLonLat
-        dataAgent.eval (f, aid)
-        |> Result.map ExtendedGrid
+    grid |> Result.map (fun g -> (g.Grid :?> Grid) |> toWebMercator CoordSys.WGS84 |> toPlainGrid)

-let private cacheExtGrid aid (grid: Result<ExtendedGrid, _>) =
-    grid
-    |> Result.map (fun x ->
-        Log.Debug("Caching archive {ArchiveId}")
-        cache.add (aid, x)
-        x)
-
-let private retrieveAndCacheExtGrid toWebMercator aid =
-    retrieveGrid aid
-    |> Result.map (fun grid ->
-        Log.Debug $"grid: {grid}"
-        try
-            if toWebMercator then
-                Log.Debug "projecting grid to WebMercator"
-                let g = grid.ToGrid ()
-                Grid.toWebMercator CoordSys.WGS84 g |> ExtendedGrid
-            else
-                grid
-        with exn ->
-            Log.Error $"{exn}"
-            failwith exn.Message
-        )
-    |> cacheExtGrid aid
+// Grid in LonLat
+let retrieveExtendedGrid (aid: Guid) =
+    let s = Diagnostics.Stopwatch()
+    s.Start()
+    let eg =
+        match cache.tryGet aid with
+        | Some (g: ExtendedGrid) -> g |> Ok
+        | None ->
+            retrieveGrid aid
+            |> Result.map (ExtendedGrid >> cacheExtGrid aid)
+    Log.Debug $"retrieveExtendedGrid(): {s.ElapsedMilliseconds}ms"
+    eg

 let private queryGridIndexCache (aid: Guid) (grid: ExtendedGrid) =
    let sha = grid.getHashString()
@@ -214,124 +296,10 @@ let private queryGridIndexCache (aid: Guid) (grid: ExtendedGrid) =
    }
    |> Async.Start

-/// Compute density of seawater from salinity, temperature, and pressure
-let private dens (temp: float) (salt: float) (pres: float) =
-    // From Bjørn Aadlandsvik https://github.com/bjornaa/seawater/blob/master/seawater/density.py
-
-    /// Density of seawater at zero pressure
-    let dens0 (salt: float) (temp: float) =
-        // --- Define constants ---
-        let a0 = 999.842594
-        let a1 = 6.793952e-2
-        let a2 = -9.095290e-3
-        let a3 = 1.001685e-4
-        let a4 = -1.120083e-6
-        let a5 = 6.536332e-9
-
-        let b0 = 8.24493e-1
-        let b1 = -4.0899e-3
-        let b2 = 7.6438e-5
-        let b3 = -8.2467e-7
-        let b4 = 5.3875e-9
-
-        let c0 = -5.72466e-3
-        let c1 = 1.0227e-4
-        let c2 = -1.6546e-6
-
-        let d0 = 4.8314e-4
-
-        // --- Computations ---
-        // Density of pure water
-        let SMOW = a0 + (a1 + (a2 + (a3 + (a4 + a5 * temp) * temp) * temp) * temp) * temp
-
-        // More temperature polynomials
-        let RB = b0 + (b1 + (b2 + (b3 + b4 * temp) * temp) * temp) * temp
-        let RC = c0 + (c1 + c2 * temp) * temp
-
-        SMOW + RB * salt + RC * (salt ** 1.5) + d0 * salt * salt
-
-    /// Secant bulk modulus
-    let seck salt temp pres =
-        // --- Pure water terms ---
-        let h0 = 3.239908
-        let h1 = 1.43713e-3
-        let h2 = 1.16092e-4
-        let h3 = -5.77905e-7
-        let AW = h0 + (h1 + (h2 + h3 * temp) * temp) * temp
-
-        let k0 = 8.50935e-5
-        let k1 = -6.12293e-6
-        let k2 = 5.2787e-8
-        let BW = k0 + (k1 + k2 * temp) * temp
-
-        let e0 = 19652.21
-        let e1 = 148.4206
-        let e2 = -2.327105
-        let e3 = 1.360477e-2
-        let e4 = -5.155288e-5
-        let KW = e0 + (e1 + (e2 + (e3 + e4 * temp) * temp) * temp) * temp
-
-        // --- Seawater, P = 0 ---
-        let SR = salt ** 0.5
-
-        let i0 = 2.2838e-3
-        let i1 = -1.0981e-5
-        let i2 = -1.6078e-6
-        let j0 = 1.91075e-4
-        let A = AW + (i0 + (i1 + i2 * temp) * temp + j0 * SR) * salt
-
-        let f0 = 54.6746
-        let f1 = -0.603459
-        let f2 = 1.09987e-2
-        let f3 = -6.1670e-5
-        let g0 = 7.944e-2
-        let g1 = 1.6483e-2
-        let g2 = -5.3009e-4
-        let K0 =
-            KW
-            + (f0 + (f1 + (f2 + f3 * temp) * temp) * temp + (g0 + (g1 + g2 * temp) * temp) * SR)
-              * salt
-
-        // --- General expression ---
-        let m0 = -9.9348e-7
-        let m1 = 2.0816e-8
-        let m2 = 9.1697e-10
-        let B = BW + (m0 + (m1 + m2 * temp) * temp) * salt
-
-        K0 + (A + B * pres) * pres
-
-    // Convert to bar
-    let pres = 0.1 * pres
-    (dens0 salt temp) / (1.0 - pres / seck salt temp pres)
-
-let nLayers (aid: Guid) =
-    let f ds _ = Fvcom.getNumSiglay ds
-    dataAgent.eval (f, aid)
-
-let retrievePlainGrid (aid: Guid) =
-    match cache.tryGet aid with
-    | Some(g: ExtendedGrid) ->
-        Log.Debug $"plainGrid: in cache {aid}"
-        g.Grid |> toPlainGrid |> Ok
-    | None ->
-        Log.Debug $"plainGrid: fetch {aid}"
-        retrieveAndCacheExtGrid true aid
-        |> Result.map (fun eg -> eg.Grid |> toPlainGrid)
-
-let retrieveExtendedGrid toWebMercator (aid: Guid) =
-    let s = Diagnostics.Stopwatch()
-    s.Start()
-    let r =
-        match cache.tryGet aid with
-        | Some (g: ExtendedGrid) -> g |> Ok
-        | None -> retrieveAndCacheExtGrid toWebMercator aid
-    Log.Debug $"obtainExtendedGrid(): {s.ElapsedMilliseconds}ms"
-    r
-
 let equipArchiveCache (aid: Guid) =
    Log.Debug("Trying to equip archive {ArchiveId}", aid)
    try
-        retrieveExtendedGrid true aid
+        retrieveExtendedGrid aid
        |> function
            | Ok g ->
                let sha = g.getHashString ()
@@ -346,6 +314,43 @@ let equipArchiveCache (aid: Guid) =
        Log.Error $"sorcerer: equipArchviveCache: {exn}"
        None

+// Grid in LonLat
+let private tryObtainExtGrid aid () =
+    retrieveExtendedGrid aid
+    |> Result.toOption
+
+// Coord in LonLat
+let nearestElement aid coord =
+    cache.tryGet aid
+    |> Option.orElseWith (tryObtainExtGrid aid)
+    |> Option.bind (fun cg ->
+        let eIdx = cg.tryGetElement coord
+        Log.Debug $"found eidx: %A{eIdx}"
+        eIdx)
+
+// Coord in LonLat
+let nearestNode aid coord =
+    cache.tryGet aid
+    |> Option.orElseWith (tryObtainExtGrid aid)
+    |> Option.bind (fun cg ->
+        let nIdx = cg.nearestNode coord
+        Log.Debug $"found nidx: %A{nIdx}"
+        nIdx)
+
+// Coord in LonLat
+let nearestElements aid coord =
+    cache.tryGet aid
+    |> Option.orElseWith (tryObtainExtGrid aid)
+    |> Option.map (fun eg -> coord |> Array.map eg.tryGetElement)
+    |> Option.defaultValue Array.empty
+
+// Coord in LonLat
+let nearestNodes aid coord =
+    cache.tryGet aid
+    |> Option.orElseWith (tryObtainExtGrid aid)
+    |> Option.map (fun eg -> coord |> Array.map eg.tryGetNode)
+    |> Option.defaultValue Array.empty
+
 let numFrames aid =
    archiveAgent.tryGetArchive aid
    |> Option.map (fun (a: ArchiveDto) -> a.props.frames)
@@ -365,7 +370,7 @@ let speedAtElement aid t e =

 let speedAtNodes aid t l =
    let f ds t = Fvcom.readUV ds t l
-    let grid = retrieveExtendedGrid true aid |> Result.get
+    let grid = retrieveExtendedGrid aid |> Result.get

    let nix =
        grid.NeighborIndex
@@ -375,7 +380,7 @@ let speedAtNodes aid t l =

    dataAgent.eval (f, aid, t)
    |> Result.map (fun uvs ->
-        let uvs' = uvs |> Array.map (bimap float)
+        let uvs' = uvs |> Array.map (Utils.bimap float)
        Util.Element.speedToNodal nix uvs'
        |> Array.map single)

@@ -383,7 +388,7 @@ let uvAtNodes aid t l =
    let f ds t = Fvcom.readUV ds t l

    let nix =
-        let grid = retrieveExtendedGrid true aid |> Result.get
+        let grid = retrieveExtendedGrid aid |> Result.get
        grid.NeighborIndex
        |> Option.defaultWith (fun () ->
            reviewNix aid grid
@@ -391,17 +396,17 @@ let uvAtNodes aid t l =

    dataAgent.eval (f, aid, t)
    |> Result.map (fun uvs ->
-        let us', vs' = uvs |> Array.map (bimap float) |> Array.unzip
+        let us', vs' = uvs |> Array.map (Utils.bimap float) |> Array.unzip
        let us'' = Util.Element.propToNodal nix us'
        let vs'' = Util.Element.propToNodal nix vs'
        Array.zip us'' vs''
-        |> Array.map (bimap single))
+        |> Array.map (Utils.bimap single))

 let uAtNodes aid t l =
    let f ds t = Fvcom.readU ds t l

    let nix =
-        let grid = retrieveExtendedGrid true aid |> Result.get
+        let grid = retrieveExtendedGrid aid |> Result.get
        grid.NeighborIndex
        |> Option.defaultWith (fun () ->
            reviewNix aid grid
@@ -417,7 +422,7 @@ let vAtNodes aid t l =
    let f ds t = Fvcom.readV ds t l

    let nix =
-        let grid = retrieveExtendedGrid true aid |> Result.get
+        let grid = retrieveExtendedGrid aid |> Result.get
        grid.NeighborIndex
        |> Option.defaultWith (fun () ->
            reviewNix aid grid
@@ -497,66 +502,11 @@ let archiveBBox aid =

 let archiveTime aid n =
    let f ds _ = Fvcom.getTime ds n
-
    dataAgent.eval (f, aid) |> Result.get

-let private obtainExtGrid aid () =
-    retrieveAndCacheExtGrid true aid
-    |> Result.toOption
-
-let elem aid x y =
-    cache.tryGet aid
-    |> Option.orElseWith (obtainExtGrid aid)
-    |> Option.bind (fun cg ->
-        Log.Debug $"x: {x}, y: {y}, as singles x: {single x}, y: {single y}"
-        let eidx = cg.tryGetElement (x, y)
-        Log.Debug $"found eidx: %A{eidx}"
-        eidx)
-
-let elemLonLat aid x y =
-    let toWebMercator = makeTransform CoordSys.WGS84 CoordSys.EPSG3857
-    let coord: Coordf = x, y
-    let x, y = toWebMercator.project coord
-    cache.tryGet aid
-    |> Option.orElseWith (obtainExtGrid aid)
-    |> Option.bind (fun cg ->
-        Log.Debug $"x: {x}, y: {y}, as singles x: {single x}, y: {single y}"
-        let eidx = cg.tryGetElement (x, y)
-        Log.Debug $"found eidx: %A{eidx}"
-        eidx)
-
-let nodeLonLat aid x y =
-    let toWebMercator = makeTransform CoordSys.WGS84 CoordSys.EPSG3857
-    let coord: Coordf = x, y
-    let x, y = toWebMercator.project coord
-    cache.tryGet aid
-    |> Option.orElseWith (obtainExtGrid aid)
-    |> Option.bind (fun cg ->
-        let nodeIdx = cg.nearestNode (x, y)
-        Log.Debug $"found node idx: %A{nodeIdx}"
-        nodeIdx)
-
-let elems aid coords : int option array =
-    cache.tryGet aid
-    |> Option.orElseWith (obtainExtGrid aid)
-    |> function
-        | Some cg -> coords |> Array.map cg.tryGetElement
-        | None -> Array.empty
-
-let node aid x y =
-    cache.tryGet aid
-    |> Option.orElseWith (obtainExtGrid aid)
-    |> Option.bind (fun cg ->
-        let nodeIdx = cg.nearestNode (x, y)
-        Log.Debug $"found node idx: %A{nodeIdx}"
-        nodeIdx
-    )
-
-let nodes aid coords =
-    cache.tryGet aid
-    |> Option.orElseWith (obtainExtGrid aid)
-    |> Option.map (fun cg -> coords |> Array.map cg.tryGetNode)
-    |> Option.defaultValue [||]
+let nLayers (aid: Guid) =
+    let f ds _ = Fvcom.getNumSiglay ds
+    dataAgent.eval (f, aid)

 let zetaAtNode aid t n =
    let f ds t = Fvcom.Singular.readZeta ds n t
@@ -566,20 +516,17 @@ let densityAtNode aid t n =
    monad {
        let! temp = tempAtNode aid t n
        let! salt = salinityAtNode aid t n
-        let! h = bathymetryAtNode aid n
-        let! zeta = zetaAtNode aid t n
-        let! siglay = siglay aid n

        return
-            Array.map3 (fun t_val s_val sigma ->
-                let depth = abs (sigma * h + zeta) |> float
+            Array.map2 (fun t_val s_val ->
                let t_val = float t_val
                let s_val = float s_val
+
                // NOTE: Pressure set to surface pressure which is
                // approx 0 using σ_t(T, S, 0)
-                let pres = 0
-                (dens t_val s_val pres - 1000.) |> single
-            ) temp salt siglay
+                let pres = 0.
+                (calcDensity t_val s_val pres - 1000.) |> single
+            ) temp salt
    }

 let getDepthsAtNode aid node =
@@ -596,31 +543,90 @@ let getDepthsAtElement aid elem =
        return siglay |> Array.map (fun z -> -z * h)
    } |> Result.get

-module Proximity =
-    let nodesInCircle aid lon lat radius : NodeIdx array =
-        let gridSys =
-            archiveAgent.tryGetArchive aid
-            |> Option.bind (_.props.projection >> Transformations.stringToTransformation)
-            |> Option.defaultWith (fun () -> failwith $"retrieveGrid: failed to get projecton coordsys for : {aid}")

-        let otmSys = CoordSys.OTMn (int lon)
-        let lonLatToGrid = Coordf >> (makeTransform CoordSys.WGS84 gridSys).project
-        let gridToOTM = Coordf >> (makeTransform gridSys otmSys).project
+module Proximity =
+    /// AI generated:
+    /// Returns true if a point (x, y) is inside the polygon.
+    /// - polygon: array of vertices in order (either clockwise or counter‑clockwise).
+    /// - point: the point to test.
+    /// Includes points on the boundary as "inside".
+    let private pointInPolygon (polygon: (float * float) array) (px, py: float) : bool =
+        let n = polygon.Length
+        if n < 3 then
+            false
+        else
+            let mutable inside = false
+
+            // Iterate over each edge (xi, yi) -> (xj, yj)
+            let mutable i = 0
+            let mutable j = n - 1
+
+            while i < n do
+                let (xi, yi) = polygon[i]
+                let (xj, yj) = polygon[j]
+
+                // Check if the edge crosses the horizontal ray from the point
+                let intersect =
+                    ((yi > py) <> (yj > py)) &&
+                    (px < (xj - xi) * (py - yi) / (yj - yi + 0.0) + xi)
+
+                if intersect then
+                    inside <- not inside
+
+                j <- i
+                i <- i + 1
+
+            inside
+
+    let nodesInPolygon aid (poly: (float*float) array) : NodeIdx array =
+        let lon = poly |> Array.averageBy fst
+        let lat = poly |> Array.averageBy snd

        let grid : ExtendedGrid =
-            match retrieveGrid aid with
+            match retrieveExtendedGrid aid with
            | Ok g -> g
            | _ -> failwith "error reading grid"

-        let r0 = (lon, lat) |> gridToOTM
        let getVertex = (grid :> IGrid).getVertex
-        let calcDistance (x1, y1) (x0, y0) = ((x1 - x0)**2.0 + (y1 - y0)**2.0) |> sqrt
-        let calcNodeDistance = getVertex >> gridToOTM >> (calcDistance r0)
+        let nodeInPolygon = getVertex >> pointInPolygon poly
+        let n0 =
+            (lon, lat)
+            |> grid.nearestNode
+            |> Option.defaultWith (fun () -> failwith "nearestNode: failed to find nearest node")
+
+        let rec getSurroundingNodes (allNodes: NodeIdx array) (newNodes: NodeIdx array) : NodeIdx array =
+            if Array.isEmpty newNodes then
+                allNodes
+            else
+                let newNodes' =
+                    newNodes
+                    |> Array.map grid.getNodesSurroundingNode
+                    |> Array.concat
+                    |> Array.distinct
+                    |> Array.filter nodeInPolygon
+                    |> Array.filter (fun i -> not (Array.contains i allNodes))
+
+                let allNodes' =
+                    newNodes'
+                    |> Array.append allNodes
+                    |> Array.distinct
+
+                getSurroundingNodes allNodes' newNodes'
+
+        getSurroundingNodes [| n0 |] [| n0 |]
+
+    let nodesInCircle aid lon lat radius : NodeIdx array =
+        let grid : ExtendedGrid =
+            match retrieveExtendedGrid aid with
+            | Ok g -> g
+            | _ -> failwith "error reading grid"
+
+        let getVertex = (grid :> IGrid).getVertex
+        let calcNodeDistance = getVertex >> Utils.haversineDistance (lon, lat)
        let filterNodeDistance = calcNodeDistance >> (>) radius

        let n0 =
            (lon, lat)
-            |> lonLatToGrid
            |> grid.nearestNode
            |> Option.defaultWith (fun () -> failwith "nearestNode: failed to find nearest node")

@@ -660,10 +666,6 @@ module Proximity =
        (startElem: ElemIdx, endElem: ElemIdx)
        : (float * ElemIdx array) option =

-        /// Inlined calculation of squared distance
-        let inline calcDistanceSquared (x1, y1) (x0, y0) =
-            ((x1 - x0) ** 2.0 + (y1 - y0) ** 2.0) |> sqrt
-
        let distances = Dictionary<ElemIdx, float> ()
        let previous = Dictionary<ElemIdx, ElemIdx> ()

@@ -695,8 +697,8 @@ module Proximity =
                    | Some neighs ->
                        for neighbor in neighs do
                            if not (visited.Contains neighbor) then
-                                let distSq = calcDistanceSquared centroids[elem] centroids[neighbor]
-                                let newDistance = distanceTraveled + distSq
+                                let distance = Utils.haversineDistance centroids[elem] centroids[neighbor]
+                                let newDistance = distanceTraveled + distance

                                match distances.TryGetValue neighbor with
                                | true, oldDistance when newDistance >= oldDistance -> ()
@@ -704,8 +706,8 @@ module Proximity =
                                    distances[neighbor] <- newDistance
                                    previous[neighbor] <- elem

-                                    let distSq = calcDistanceSquared centroids[neighbor] centroids[endElem]
-                                    let priority = newDistance + distSq
+                                    let distance = Utils.haversineDistance centroids[neighbor] centroids[endElem]
+                                    let priority = newDistance + distance
                                    queue.Enqueue(neighbor, priority)

        if found then
@@ -738,10 +740,6 @@ module Proximity =
        (startNode: NodeIdx, endNode: NodeIdx)
        : (float * NodeIdx array) option =

-        /// Inlined calculation of squared distance
-        let inline calcDistanceSquared (x1, y1) (x0, y0) =
-            ((x1 - x0) ** 2.0 + (y1 - y0) ** 2.0) |> sqrt
-
        let distances = Dictionary<NodeIdx, float> ()
        let previous = Dictionary<NodeIdx, NodeIdx> ()

@@ -773,8 +771,8 @@ module Proximity =
                    | Some neighs ->
                        for neighbor in neighs do
                            if not (visited.Contains neighbor) then
-                                let distSq = calcDistanceSquared vertices[node] vertices[neighbor]
-                                let newDistance = distanceTraveled + distSq
+                                let distance = Utils.haversineDistance vertices[node] vertices[neighbor]
+                                let newDistance = distanceTraveled + distance

                                match distances.TryGetValue neighbor with
                                | true, oldDistance when newDistance >= oldDistance -> ()
@@ -782,8 +780,8 @@ module Proximity =
                                    distances[neighbor] <- newDistance
                                    previous[neighbor] <- node

-                                    let distSq = calcDistanceSquared vertices[neighbor] vertices[endNode]
-                                    let priority = newDistance + distSq
+                                    let distance = Utils.haversineDistance vertices[neighbor] vertices[endNode]
+                                    let priority = newDistance + distance
                                    queue.Enqueue(neighbor, priority)

        if found then
@@ -803,45 +801,39 @@ module Proximity =


    let shortestElementPath aid r0 r1 : GridPath option =
-        let gridToLonLat = Coordf >> (makeTransform CoordSys.EPSG3857 CoordSys.WGS84).project
-        let lonLatToGrid = Coordf >> (makeTransform CoordSys.WGS84 CoordSys.EPSG3857).project
+        match retrieveExtendedGrid aid with
+        | Ok grid ->
+            let eOpt0 = r0 |> grid.tryGetElement
+            let eOpt1 = r1 |> grid.tryGetElement

-        // Fetches in WebMercator (EPSG3857)
-        let gridRes = retrieveExtendedGrid true aid
-        let eOpt0 = [| r0 |> lonLatToGrid |] |> elems aid |> Array.tryHead |> Option.flatten
-        let eOpt1 = [| r1 |> lonLatToGrid |] |> elems aid |> Array.tryHead |> Option.flatten
+            match eOpt0, eOpt1 with
+            | Some e0, Some e1 ->
+                let centroids = grid.getCentroids()
+                let getCentroid i = centroids[i]

-        match eOpt0, eOpt1, gridRes with
-        | Some e0, Some e1, Ok grid ->
-            let centroids = grid.getCentroids()
-            let getCentroid i = centroids[i]
+                let nix =
+                    grid.NeighborIndex
+                    |> Option.defaultWith (fun () ->
+                        reviewNix aid grid
+                        grid.NeighborIndex.Value
+                    )

-            let nix =
-                grid.NeighborIndex
-                |> Option.defaultWith (fun () ->
-                    reviewNix aid grid
-                    grid.NeighborIndex.Value
+                (e0, e1)
+                |> shortestPathBetweenElements nix centroids
+                |> Option.map (fun (dist, elems) ->
+                    {
+                        length = dist
+                        lineSegments = elems |> Array.map getCentroid
+                    }
                )
-
-            (e0, e1)
-            |> shortestPathBetweenElements nix centroids
-            |> Option.map (fun (dist, elems) ->
-                {
-                    length = dist
-                    lineSegments = elems |> Array.map (getCentroid >> gridToLonLat)
-                }
-            )
+            | _ -> None
        | _ -> None

    let shortestNodePath aid r0 r1 : GridPath option =
-        let gridToLonLat = Coordf >> (makeTransform CoordSys.EPSG3857 CoordSys.WGS84).project
-        let lonLatToGrid = Coordf >> (makeTransform CoordSys.WGS84 CoordSys.EPSG3857).project
-
-        // Fetches in WebMercator (EPSG3857)
-        match retrieveExtendedGrid true aid with
+        match retrieveExtendedGrid aid with
        | Ok grid ->
-            let nOpt0 = r0 |> lonLatToGrid |> grid.nearestNode
-            let nOpt1 = r1 |> lonLatToGrid |> grid.nearestNode
+            let nOpt0 = r0 |> grid.nearestNode
+            let nOpt1 = r1 |> grid.nearestNode

            match nOpt0, nOpt1 with
            | Some n0, Some n1 ->
@@ -860,7 +852,7 @@ module Proximity =
                |> Option.map (fun (dist, nodes) ->
                    {
                        length = dist
-                        lineSegments = nodes |> Array.map (getVertex >> gridToLonLat)
+                        lineSegments = nodes |> Array.map getVertex
                    }
                )
            | _ -> None
@@ -880,10 +872,6 @@ module Proximity =
        (maxDistance: float)
        : (ElemIdx * float) array =

-        /// Inlined calculation of squared distance
-        let inline calcDistanceSquared (x1, y1) (x0, y0) =
-            (x1 - x0) ** 2.0 + (y1 - y0) ** 2.0 |> sqrt
-
        let distances = Dictionary<ElemIdx, float> ()
        let queue = PriorityQueue<ElemIdx, float> ()

@@ -907,13 +895,13 @@ module Proximity =
                    | None -> ()
                    | Some neighs ->
                        for neighbor in neighs do
-                            let distSq = calcDistanceSquared centroids[elem] centroids[neighbor]
-                            let newDistance = distanceTraveled + distSq
+                            let distance = Utils.haversineDistance centroids[elem] centroids[neighbor]
+                            let newDistance = distanceTraveled + distance
                            let _, oldDistance = distances.TryGetValue neighbor

                            if newDistance <= maxDistance && (not (distances.ContainsKey neighbor) || newDistance < oldDistance) then
                                distances[neighbor] <- newDistance
-                                let priority = newDistance + distSq
+                                let priority = newDistance + distance
                                queue.Enqueue(neighbor, priority)

        // Return all triangles in the distance map which are within the maxDistance
@@ -934,19 +922,10 @@ module Proximity =
                if not (boundary.Remove(u, v)) then
                    boundary.Add(u, v) |> ignore

-
-    let private mercatorScaleFactor (lat: float) =
-        1.0 / (lat * Math.PI / 180.0 |> Math.Cos)
-
    let gridDistanceCircle aid lon lat radius : (float * float) array =
-        let gridToLonLat = Coordf >> (makeTransform CoordSys.EPSG3857 CoordSys.WGS84).project
-        let lonLatToGrid = Coordf >> (makeTransform CoordSys.WGS84 CoordSys.EPSG3857).project
-        let radius = radius * mercatorScaleFactor lat
-
-        // Fetches in WebMercator (EPSG3857)
-        match retrieveExtendedGrid true aid with
+        match retrieveExtendedGrid aid with
        | Ok grid ->
-            let eOpt = [| (lon, lat) |> lonLatToGrid |] |> elems aid |> Array.tryHead |> Option.flatten
+            let eOpt = (lon, lat) |> grid.tryGetElement

            match eOpt with
            | Some eIdx ->
@@ -978,16 +957,12 @@ module Proximity =
                    adjacency[j].Add(i)

                let centerPoint = centroids[eIdx]
-                let distSq (vIdx: NodeIdx) =
-                    let vx, vy = getVertex vIdx
-                    let dx = vx - fst centerPoint
-                    let dy = vy - snd centerPoint
-                    dx * dx + dy * dy
+                let calcHaversineDistance = getVertex >> Utils.haversineDistance centerPoint

                // Make sure not to start on an island boundary
                let startNode =
                    adjacency.Keys
-                    |> Seq.maxBy distSq
+                    |> Seq.maxBy calcHaversineDistance

                // Traverse the boundary to create an ordered outline
                let orderedOutline = ResizeArray<NodeIdx>()
@@ -1007,7 +982,7 @@ module Proximity =
                                    if Seq.isEmpty x then
                                        None
                                    else
-                                        x |> Seq.maxBy distSq |> Some
+                                        x |> Seq.maxBy calcHaversineDistance |> Some

                            match next with
                            | Some n ->
@@ -1023,7 +998,7 @@ module Proximity =
                        | _ -> continue' <- false

                orderedOutline
-                |> Seq.map (getVertex >> gridToLonLat)
+                |> Seq.map getVertex
                |> Array.ofSeq
            | _ -> [||]
        | _ -> [||]
@@ -1189,16 +1164,15 @@ module TimeSeries =
                    monad {
                        let! temp = tempAtNode aid frame idx vidx
                        let! salt = salinityAtNode aid frame idx vidx
-                        let! h = bathymetryAtNode aid idx
-                        let! zeta = zetaAtNode aid frame idx

-                        let depth = 0.
                        let t_val = float temp
                        let s_val = float salt
-                        let pres = depth * 10.0

-                        let dens' = (dens t_val s_val pres - 1000.) |> single
-                        return dens'
+                        // NOTE: Pressure set to surface pressure which is
+                        // approx 0 using σ_t(T, S, 0)
+                        let pres = 0.
+
+                        (calcDensity t_val s_val pres - 1000.) |> single
                    }
                )
                |> Result.defaultValue Array.empty<single>
--- a/src/Sorcerer/src/Server/Stats.fs
+++ b/src/Sorcerer/src/Server/Stats.fs
@@ -136,7 +136,7 @@ let getSalinityByLayer (aid: Guid) (dt: Period) (m: StatMetric) (siglay: int) =
    getStatByLayer aid StatProp.Salinity dt m siglay

 let speedAtNodes aid (speed: single[])=
-    let grid = Fvcom.retrieveExtendedGrid true aid |> Result.get
+    let grid = Fvcom.retrieveExtendedGrid aid |> Result.get
    let nix =
        grid.NeighborIndex
        |> Option.defaultWith (fun () ->
--- a/src/Sorcerer/src/Server/Utils.fs
+++ b/src/Sorcerer/src/Server/Utils.fs
@@ -2,11 +2,18 @@ module Utils

 open ProjNet.FSharp
 open System
-open FSharpPlus

 let bimap f (a, b) = f a, f b

-let radToDeg x = (180. / Math.PI) * x
+let earthRadiusM = 6.371e6 // Earth's radius in meters
+let inline degToRad degrees = degrees * Math.PI / 180.0
+let inline radToDeg radians = radians * 180.0 / Math.PI
+
+let wrapLonRad lonRad =
+    let twoPi = 2.0 * Math.PI
+    let x = (lonRad + Math.PI) % twoPi
+    let x = if x < 0.0 then x + twoPi else x
+    x - Math.PI

 let tileToLatLng x y z =
    let n = 2. ** z
@@ -26,3 +33,35 @@ let lngLatToLCC =

 let mercatorScaleFactor (lat: float) =
   1.0 / (lat * Math.PI / 180.0 |> Math.Cos)
+
+let destinationPointSphere (lon, lat) distance theta =
+    let phi1 = degToRad lat
+    let lambda1 = degToRad lon
+    let delta = distance / earthRadiusM
+
+    let sinPhi1, cosPhi1 = Math.Sin phi1, Math.Cos phi1
+    let sinDelta, cosDelta = Math.Sin delta, Math.Cos delta
+    let sinTheta, cosTheta = Math.Sin theta, Math.Cos theta
+
+    let phi2 = Math.Asin (sinPhi1 * cosDelta + cosPhi1 * sinDelta * cosTheta)
+    let sinPhi2 = Math.Sin phi2
+
+    let y = sinTheta * sinDelta * cosPhi1
+    let x = cosDelta - sinPhi1 * sinPhi2
+    let lambda2 = wrapLonRad (lambda1 + Math.Atan2(y, x))
+
+    radToDeg lambda2, radToDeg phi2
+
+let haversineDistance (lon1: float, lat1: float) (lon2: float, lat2: float) =
+    let dLat = degToRad(lat2 - lat1)
+    let dLon = degToRad(lon2 - lon1)
+    let lat1 = degToRad(lat1)
+    let lat2 = degToRad(lat2)
+
+    let a = Math.Sin(dLat/2.0) * Math.Sin(dLat/2.0) +
+            Math.Sin(dLon/2.0) * Math.Sin(dLon/2.0) *
+            Math.Cos(lat1) * Math.Cos(lat2)
+
+    let c = 2.0 * Math.Asin(Math.Sqrt(a))
+    earthRadiusM * c // Returns distance in meters
+
--- a/src/Sorcerer/tilt/tilt/appsettings.json
+++ b/src/Sorcerer/tilt/tilt/appsettings.json
@@ -42,10 +42,10 @@
    },
    "plainAuthUsers": [],
    "fga": {
-      "apiUrl": "https://openfga.srv.oceanbox.io",
-      "apiKey": "",
-      "storeId": "01JKTZXMP7ANN4GG2P5W8Y56M6",
-      "modelId": "01JKTZYMCZZBVSBG66W27XMW0A"
+        "apiUrl": "https://openfga.srv.oceanbox.io",
+        "apiKey": "",
+        "storeId": "01JKTZXMP7ANN4GG2P5W8Y56M6",
+        "modelId": "01JKTZYMCZZBVSBG66W27XMW0A"
    },
    "sentryUrl": "https://5e6e3584098dc006de18038cf85d2cbe@o4509530141622272.ingest.de.sentry.io/4509547350065232",
    "redis": "localhost:6379,user=default,password=secret",
@@ -66,6 +66,6 @@
    "otelCollector": "http://10.255.241.12:4317",
    "archiveSvc": "https://<x>-atlantis.dev.oceanbox.io",
    "dataDir": "/data/archives",
-    "cacheDir": "/data/archives/cache",
+    "cacheDir": "/data/archives/cache/review",
    "authDomain": "staging"
 }
Author	SHA1	Message	Date
Stig Rune Jensen	a1e1a1bc1e	devel: update tilt settings	2025-12-17 17:32:17 +01:00
Stig Rune Jensen	f019b6839b	fix: use consistent lonlat projection on all cached grids	2025-12-17 17:30:12 +01:00
Simen Kirkvik	d6dc637476	Set cwd for atlantis/just run-client	2025-12-17 15:33:43 +01:00