Merge pull request #72 from GeoscienceAustralia/MAE

Mae

Project.toml

@@ -1,7 +1,7 @@
 name = "HiQGA"
 uuid = "01574e87-63b6-4466-925a-334cf7e21b6b"
 authors = ["richardt94 <29562583+richardt94@users.noreply.github.com>"]
-version = "0.4.4"
+version = "0.4.5"
 
 [deps]
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"

src/CommonToAll.jl

@@ -18,9 +18,9 @@ export trimxft, assembleTat1, gettargtemps, checkns, getchi2forall, nicenup, plo
         summaryconductivity, plotsummarygrids1, getVE, writevtkfromsounding, 
         readcols, colstovtk, findclosestidxincolfile, zcentertoboundary, zboundarytocenter, 
         writeijkfromsounding, nanmean, infmean, nanstd, infstd, infnanmean, infnanstd, 
-        kde_sj, plotmanygrids, readwell, getlidarheight, plotblockedwellonimages, getdeterministicoutputs, getprobabilisticoutputs, 
-        readfzipped, readxyzrhoϕ, writevtkxmlforcurtain, getRandgridr, getallxyinr, getXYlast,
-        getprobabilisticlinesfromdirectory, readxyzrhoϕnu, plotgausshist
+        kde_sj, plotmanygrids, readwell, getlidarheight, plotblockedwellonimages, getdeterministicoutputs, 
+        getprobabilisticoutputs,readfzipped, readxyzrhoϕ, writevtkxmlforcurtain, getRandgridr, getallxyinr, getXYlast,
+        getprobabilisticlinesfromdirectory, readxyzrhoϕnu, plotgausshist, findMAE, getclosestidx, getpostatXY, getpostatwell
 
 # Kernel Density stuff
 abstract type KDEtype end
@@ -790,6 +790,39 @@ function plot_posterior(F::Operator1D,
     CI[:,1], CI[:,2], CI[:,3], meanimage, meandiffimage, sdslope
 end
 
+function getpostatXY(F::Operator, soundings::Vector{T}, opt_in::OptionsStat, Xwanted, Ywanted; 
+        burninfrac = 0.5,
+        nbins = 50, 
+        qp1 = nothing, 
+        qp2 = nothing,
+        usekde = true,
+        pdfnormalize = false,
+        ) where T<:Sounding
+    @assert !isnothing(qp1) | !isnothing(qp2)
+    opt = deepcopy(opt_in)
+    M = getclosestensemble(soundings, opt, Xwanted, Ywanted; burninfrac)
+    rhomin, rhomax = extrema(opt.fbounds)
+    himage, edges, CI, _ = gethimage(F, M, opt; temperaturenum=1,
+                nbins, qp1, qp2, rhomin, rhomax, usekde,
+                islscale=false, pdfnormalize)
+    himage, edges, CI
+end
+
+function getpostatwell(F::Operator, soundings::Vector{T}, opt_in::OptionsStat, Xwell, Ywell, zc_rho; 
+    burninfrac = 0.5,
+    nbins = 50, 
+    qp1 = nothing, 
+    qp2 = nothing,
+    usekde = true,
+    pdfnormalize = false,
+    ) where T<:Sounding
+    zcentre = opt_in.xall
+    Mwell = blocktomodel(zcentre, zc_rho)
+    himage, edges, CI = getpostatXY(F, soundings, opt_in, Xwell, Ywell; burninfrac,
+                                nbins, qp1, qp2, usekde, pdfnormalize)
+    himage, edges, CI, Mwell
+end
+
 function getbounds(CI, bounds)
     propmin = min(minimum(CI), minimum(bounds))
     propmax = max(maximum(CI), maximum(bounds))
@@ -2202,7 +2235,8 @@ end
 function plotmanygrids(σ, X, Y, Z, zcentre; yl=[], xl=[],
         cmapσ="turbo", vmin=-Inf, vmax=Inf, topowidth=1, fontsize=12, spacefactor=5,
         dr=nothing, dz=nothing, plotbinning=true, δ²=1e-3, regtype=:R1, donn=false, hspace=1,
-        figsize=(10,10), smallratio=0.1, preferEright=true, delbin=15., titles=fill("", length(σ)))
+        figsize=(10,10), smallratio=0.1, preferEright=true, delbin=15., 
+        titles=fill("", length(σ)), XYprofiles = nothing, drawprofiles=true)
     @assert !isnothing(dr) # pass as variable as it is used by other functions too       
     if isnothing(dz)
         mins = [zc[1] for zc in zcentre] # depth to first centre
@@ -2221,18 +2255,37 @@ function plotmanygrids(σ, X, Y, Z, zcentre; yl=[], xl=[],
     x, y, xr, yr = get_x_y(r, m, coord_mle, preferEright)
     plotbinning && plotbinningresults(X, Y, x, y, xr, yr)
     outmap = map(zip(σ, X, Y, Z, zcentre)) do (s, xx, yy, topo, zc)
-         id = getclosestidx([xx';yy'], xr', yr', showinfo=false)
-         makegrid(s[:,id], xr, yr, topo[id]; donn, dr, zall=zc, dz)
+        id = getclosestidx([xx';yy'], xr', yr', showinfo=false)
+        if !isnothing(XYprofiles)
+            @assert size(XYprofiles, 1) == 2
+            idxatXY = getclosestidx([xx';yy'], XYprofiles[1,:]', XYprofiles[2,:]', showinfo=false)
+            zatXY =  zc
+            satXY = s[:,idxatXY]
+        else
+            zatXY = nothing 
+            satXY = nothing
+        end
+        img, gridr, gridz, topofine, R = makegrid(s[:,id], xr, yr, topo[id]; donn, dr, zall=zc, dz)
+        img, gridr, gridz, topofine, R, zatXY, satXY
     end
-    img, gridr, gridz, topofine, R = [[out[i] for out in outmap] for i in 1:5]
+    img, gridr, gridz, topofine, R, zatXY, satXY = [[out[i] for out in outmap] for i in 1:7]
     if (isinf(vmin) || isinf(vmax))
         vmin, vmax = extrema(reduce(vcat, [[extrema(s)...] for s in σ]))
     end
+    if !isnothing(XYprofiles)
+        @assert size(XYprofiles, 1) == 2
+        idx_ = getclosestidx([xr'; yr'], XYprofiles[1,:]', XYprofiles[2,:]', showinfo=false)
+        Rtoplotat = [cumulativelinedist(xr, yr)[id_] for id_ in idx_]
+    end
     imhandle = map(zip(ax, img, gridr, gridz, topofine, titles)) do (
         ax_, img_, gridr_, gridz_, topofine_, ti) 
         imhandle_ = ax_.imshow(img_, extent=[gridr_[1], gridr_[end], gridz_[end], gridz_[1]]; 
             cmap=cmapσ, aspect="auto", vmin, vmax)
         ax_.plot(gridr_, topofine_, linewidth=topowidth, "-k")
+        if (!isnothing(XYprofiles) && drawprofiles)
+            idx, = nn(KDTree(gridr_'), [Rtoplotat...]')
+            plotprofile(ax_, idx, topofine_, gridr_)
+        end    
         isempty(ti) || ax_.set_title(ti)
         imhandle_
     end
@@ -2251,7 +2304,7 @@ function plotmanygrids(σ, X, Y, Z, zcentre; yl=[], xl=[],
     cb.set_label("Log₁₀ S/m", labelpad=0)
     nicenup(fig, fsize=fontsize, h_pad=0)
     fig.subplots_adjust(hspace = all(isempty.(titles)) ? 0 : hspace)
-    xr, yr, ax # return easting northing of grid and figure axes
+    xr, yr, ax, zatXY, satXY # return easting northing of grid and figure axes
 end
 
 function getbinby(X, Y, preferEright)
@@ -2361,6 +2414,31 @@ function readwell(fname, skipstart; lidarfile=nothing)
     name, X, Y, Z, zc_rho
 end
 
+function findMAE(soundings::Vector{T}, opt_in::OptionsStat, Xwell, Ywell, zc_rho; burninfrac=0.5) where T<:Sounding
+    opt = deepcopy(opt_in)
+    zcentre = opt.xall
+    Mwell = blocktomodel(zcentre, zc_rho)
+    m = getclosestensemble(soundings, opt, Xwell, Ywell; burninfrac)
+    findMAE(m, vec(Mwell))
+end
+
+function blocktomodel(zcentre, zc_rho)
+    zboundaries = zcentertoboundary(zcentre)
+    Mwell = block1Dvalues([zc_rho[:,2]], zc_rho[:,1], [zboundaries[1:end-1] zboundaries[2:end]], :median)[:]
+    Mwell = vec([Mwell;Mwell[end,:]']) # dummy last cell in depth
+end
+
+function getclosestensemble(soundings::Vector{T}, opt::OptionsStat, Xwell, Ywell; burninfrac=0.5) where T<:Sounding
+    idx = getclosestidx(Xwell, Ywell, soundings)
+    opt.fdataname = soundings[idx].sounding_string*"_"
+    m = assembleTat1(opt, :fstar; burninfrac, temperaturenum=1)
+end
+
+function findMAE(mm::Vector{T}, Mwell::Vector{S}) where T<:Array where S<:Real
+    absdevs = [abs.(vec(m) - Mwell) for m in mm]
+    nanmean(reduce(hcat, absdevs), 2) # mean of ndepths × namples in nsamples dirn
+end
+
 function makeblockedwellimage(readwellarray, zall, xr, yr; distblank=50, dr=nothing, donn=false, displaydistanceaway=true )
     # xr, yr are the line path along which to find closest well index
     @assert !isnothing(dr)
@@ -2446,6 +2524,15 @@ function plotblockedwellonimages(ax, wellarray, zall, xr, yr; donn=false,
         @assert !isnothing(vmax)
     img, gridr, gridz, hsegs, vsegs = makeblockedwellimage(wellarray, zall, xr, yr; distblank, dr, donn)
     plotwelloutline(ax, img, hsegs, vsegs, gridr, gridz, vmin, vmax; cmap, color, linewidth)
+    wnames = [String(w[1]) for w in wellarray]
+    annotatewells(ax[end], wnames, hsegs)
+
+end
+
+function annotatewells(ax, wnames, hsegs)
+    for (h, wn) in zip(hsegs[2:2:end], wnames)
+        ax.text(h[end,1], h[end,2], wn, fontsize=8)
+    end
 end
 
 end # module CommonToAll

zz_portalcurtains/RDP.jl

@@ -1,8 +1,10 @@
 module RDP
 using LinearAlgebra, Dates, ArchGDAL, Printf, 
-    PyPlot, Images, FileIO, HiQGA, Interpolations
+    PyPlot, Images, FileIO, HiQGA, Interpolations, NearestNeighbors, DelimitedFiles
 import GeoFormatTypes as GFT
 import GeoDataFrames as GDF
+import HiQGA.transD_GP.LineRegression.getsmoothline
+
 const mpl = PyPlot.matplotlib
 const tilesize = 512
 const epsg_GDA94 = 4283
@@ -312,6 +314,17 @@ function writevtkfromxyzrhodir(nlayers::Int; src_dir="", dst_dir="", src_epsg=0,
     end
 end
 
+function writegiantfilefromxyzrhodir(nlayers::Int; src_dir="", dst_dir="", src_epsg=0)
+    lines = transD_GP.getprobabilisticlinesfromdirectory(src_dir)
+    isdir(dst_dir) || mkpath(dst_dir)
+    map(lines) do ln
+        X, Y, Z, zall, ρlow, ρmid, ρhigh, ρavg, ϕmean, ϕsdev = transD_GP.readxyzrhoϕ(ln, nlayers; pathname=src_dir)
+        plist = makeplist(X, Y)
+        latlonglist = reduce(hcat, reprojecttoGDA94(plist, src_epsg))'
+        [latlonglist src_epsg*ones(size(X)) X Y Z transD_GP.zcentertoboundary(zall)'.*ones(size(X)) ρlow' ρmid' ρhigh' ϕmean ϕsdev]
+    end
+end
+
 function writeshpfromxyzrhodir(nlayers::Int; prefix="", src_dir="", dst_dir="", src_epsg=0)
     lines = transD_GP.getprobabilisticlinesfromdirectory(src_dir)
     isdir(dst_dir) || mkpath(dst_dir)
@@ -381,4 +394,57 @@ function writeaseggdffromxyzrho(nlayers::Int; src_dir="", dst_dir="",
     end
 end
 
+mutable struct XY
+    x
+    y
+end 
+
+function collectpoints(;npoints=1000)
+    xy = ginput(npoints, timeout=0)
+    x, y = [[pts[i] for pts in xy] for i =  1:2]
+    XY(x, y)
+end
+
+function smoothline(xy::XY; λ²=0.01, finefactor=100, regtype=:R1, fname=nothing)
+    xmin, xmax = extrema(xy.x)
+    Δx = (xmax - xmin)/finefactor
+    gridx = range(start=xmin, stop=xmax, step=Δx)
+    gridy = snaptogrid(gridx, xy.x, xy.y)
+    sd = 1 # identity weighting matrix for points
+    ysmooth = getsmoothline(gridy, sd; δ²=λ², regtype)
+    if !isnothing(fname)
+        @assert !isfile(fname)
+        io = open(fname,"w")
+        write(io, "$(λ²)\n")
+        write(io, "$finefactor\n")
+        write(io, "$regtype\n")
+        writedlm(io, [xy.x xy.y])
+        close(io)
+    end
+    gridx, ysmooth
+end    
+
+function readpoints(fname::String)
+    io = open(fname, "r")
+    λ² = parse(Float64, readline(io))
+    finefactor = parse(Float64, readline(io))
+    regtype = Symbol(readline(io))
+    xy = readdlm(io)
+    xy = XY(xy[:,1], xy[:,2])
+    smoothline(xy; λ², finefactor, regtype)
+end
+
+function readpoints(fnames::Vector{String})
+    xy = map(fnames) do fn
+        readpoints(fn)
+    end
+end
+
+function snaptogrid(gridx, x, y)
+    idx, _ = nn(KDTree(gridx'), x')
+    gridy = NaN .+ zeros(size(gridx))
+    gridy[idx] = y
+    gridy
+end
+
 end # module

zz_portalcurtains/comparelei.jl

@@ -26,5 +26,5 @@ zall = [zalld, zallp, zallp, zallp]
 titles = ["Deterministic conductivity", "10th Percentile conductivity", "50th Percentile conductivity", "90th Percentile conductivity"]
 xrd, yrd, axd = transD_GP.plotmanygrids(deepcopy(σ), deepcopy(X), deepcopy(Y), deepcopy(Z), deepcopy(zall); 
     dr, vmin, vmax, donn, xl, yl, δ²=linesmoothδ², figsize=(20,10), titles, 
-    preferEright=true, plotbinning=true, fontsize=10, delbin, hspace=0.22, spacefactor=0.1)
+    preferEright=true, plotbinning=false, fontsize=10, delbin, hspace=0.22, spacefactor=0.1)
 savefig("Line_$linewanted.png", dpi=500)

zz_portalcurtains/makegiantfile_from_all_xyzrho.jl

@@ -0,0 +1,19 @@
+using HiQGA, PyPlot, DelimitedFiles
+includet("RDP.jl")
+nlayers = 52
+vmin, vmax = -2.5, 0.5
+items = [item for item in walkdir("/Users/anray/Documents/work/projects/largeaem/final_01/summaries_all")]
+idx_summary = [basename(it[1]) == "summary" for it in items]
+src_dir = [it[1] for it in items[idx_summary]]
+src_epsg = [28353, 28354, 28351, 28354, 28354, 28354, 28354, 28352, 28352, 28352, 28352]
+out = map(zip(src_dir, src_epsg)) do (dir, epsg)
+    dst_dir = joinpath(dirname(dir),"GDA94_vtk")
+    RDP.writegiantfilefromxyzrhodir(nlayers; src_dir=dir, dst_dir, src_epsg=epsg)
+end
+## write to file
+fname = "/Users/anray/Documents/work/projects/largeaem/final_01/giantfile.txt"
+format = "lat long src_epsg X Y Z zboundaries log10ρlow log10ρmid log10ρhigh ϕmean ϕsdev"
+write(fname, format)
+io = open(fname, "a")
+writedlm(io, reduce(vcat, reduce(vcat, out)))
+close(io)

zz_portalcurtains/testdraw.jl

@@ -0,0 +1,30 @@
+using PyPlot
+includet("/home/547/ar0754/.julia/dev/HiQGA/zz_portalcurtains/RDP.jl")
+##
+XY = RDP.collectpoints()
+x, y = RDP.smoothline(XY, λ²=10000, finefactor=1000, regtype=:R2, fname="base_menindee.txt")
+##
+XY1 = RDP.collectpoints()
+x, y = RDP.smoothline(XY1, λ²=10000, finefactor=1000, regtype=:R2, fname="base_1_menindee.txt")
+##
+XY2a = RDP.collectpoints()
+x, y = RDP.smoothline(XY2a, λ²=10000, finefactor=1000, regtype=:R2, fname="base_2a_menindee.txt")
+##
+XY2b = RDP.collectpoints()
+x, y = RDP.smoothline(XY2b, λ²=10000, finefactor=1000, regtype=:R2, fname="base_2b_menindee.txt")
+## 
+XY3a = RDP.collectpoints()
+x, y = RDP.smoothline(XY1a, λ²=10000, finefactor=1000, regtype=:R2, fname="base_3a_menindee.txt")
+## 
+XY3b = RDP.collectpoints()
+x, y = RDP.smoothline(XY1b, λ²=10000, finefactor=1000, regtype=:R2, fname="base_3b_menindee.txt")
+## read them all
+x, y = [[xy[i] for xy in RDP.readpoints(["base_menindee.txt", "base_1_menindee.txt", "base_2a_menindee.txt", "base_2b_menindee.txt", 
+                    "base_3a_menindee.txt", "base_3b_menindee.txt"])] for i in 1:2]
+map(axp) do ax
+    for a in ax[1:6]
+        for (xx, yy) in zip(x, y)
+            a.plot(xx, yy, "--k", linewidth=0.5)
+        end
+    end
+end