read in Seb files

zz_portalcurtains/RDP.jl

@@ -1,5 +1,5 @@
 module RDP
-using LinearAlgebra, Dates, ArchGDAL, Printf, 
+using LinearAlgebra, Dates, ArchGDAL, Printf, DataFrames, CSV,
     PyPlot, Images, FileIO, HiQGA, Interpolations, NearestNeighbors, DelimitedFiles
 import GeoFormatTypes as GFT
 import GeoDataFrames as GDF
@@ -446,4 +446,19 @@ function snaptogrid(gridx, x, y)
     gridy
 end
 
+function readegspoints(fname::String, dict::Dict)
+    musthave = ("Xkey", "Ykey", "elevkey", "segkey")
+    map() do key
+        @assert haskey(dict, key)
+    end
+    headers  = map(k->get(dict, k, 0), musthave)
+    df = DataFrame(CSV.File(fname))
+    out = map(headers) do h # X, Y, elev, seg
+        df[:,h]
+    end
+    # split each segment
+    segs = unique(out[end]) # the last one is segment ID
+    X, Y, elev, seg = [[o[out[end] .== s] for s in segs] for o in out]
+end
+
 end # module

zz_portalcurtains/getsebpoints.jl

@@ -0,0 +1,28 @@
+using PyPlot
+includet("/home/547/ar0754/.julia/dev/HiQGA/zz_portalcurtains/RDP.jl")
+cd(@__DIR__)
+dict = Dict("Xkey"=>"X", "Ykey"=>"Y", "elevkey"=>"ELEVATION", "segkey"=>"SegmentID")
+X, Y, elev, seg = RDP.readegspoints("Menindee_Lakes_test_line_interp_MGA54_20240509.egs", dict)
+## get xy corresponding to fine grid
+using NearestNeighbors
+R, gridr = transD_GP.getRandgridr(xrd, yrd, dr)
+ncols = length(gridr)
+xyfine, gridrfine = transD_GP.getallxyinr(xrd, yrd, 0.; rangelenth=ncols)
+tree = KDTree(xyfine)
+## get axp, axd first
+# probabilistic
+map(axp) do ax
+    for a in ax[1:6]
+        for (xx, yy, mahd) in zip(X, Y, elev)
+            idxs, = nn(tree,[xx';yy'])
+            a.plot(gridrfine[idxs], mahd, "--k", linewidth=2)
+        end
+    end
+end
+# deterministic
+for a in axd[1:6]
+    for (xx, yy, mahd) in zip(X, Y, elev)
+        idxs, = nn(tree,[xx';yy'])
+        a.plot(gridrfine[idxs], mahd, "--k", linewidth=2)
+    end
+end