zboundarytocenter_inexact added

src/CommonToAll.jl

@@ -16,7 +16,7 @@ export trimxft, assembleTat1, gettargtemps, checkns, getchi2forall, nicenup, plo
         plotprofile, gridpoints, splitsoundingsbyline, getsoundingsperline, docontinue, linestartend,
         compatidxwarn, dfn2hdr, getgdfprefix, readlargetextmatrix, pairinteractionplot, flipline, 
         summaryconductivity, plotsummarygrids1, getVE, writevtkfromsounding, trapezoidrule,
-        readcols, colstovtk, findclosestidxincolfile, zcentertoboundary, zboundarytocenter, 
+        readcols, colstovtk, findclosestidxincolfile, zcentertoboundary, zboundarytocenter, zboundarytocenter_inexact,
         writeijkfromsounding, nanmean, infmean, nanstd, infstd, infnanmean, infnanstd, 
         kde_sj, plotmanygrids, readwell, getlidarheight, plotblockedwellonimages, getdeterministicoutputs, writeasegdfnfromonerow,
         writeasegdat, getprobabilisticoutputs, readfzipped, readxyzrhoϕ, writevtkxmlforcurtain, getRandgridr, getallxyinr, getXYlast,
@@ -1279,15 +1279,15 @@ function zcentertoboundary(zall)
     zb
 end
 
-# function zboundarytocenter(zb; fudgelast=false)
+function zboundarytocenter_inexact(zb; fudgelast=false)
 # no more fudging, this is superseded
-#     thickness = diff(zb)
-#     zall = zb[1:end-1] + thickness/2
-#     if fudgelast
-#         zall = [zall; zb[end]+thickness[end]/2]
-#     end
-#     zall    
-# end    
+    thickness = diff(zb)
+    zall = zb[1:end-1] + thickness/2
+    if fudgelast
+        zall = [zall; zb[end]+thickness[end]/2]
+    end
+    zall    
+end    
 
 function zboundarytocenter(zb)
     # first get extendfrac r

zz_portalcurtains/comparelei.jl

@@ -12,7 +12,7 @@ line, X, Y, Z, thick, σ = transD_GP.readcols([5,7,8,9,[53,82],[23,52]],
     startfrom=1, decfactor=5)
 idx = line .== linewanted
 Xd, Yd, Zd, thick, σd = map(x->x[idx,:],(X, Y, Z, thick, σ))
-zalld = transD_GP.zboundarytocenter(cumsum(thick[1,:]))
+zalld = transD_GP.zboundarytocenter_inexact(cumsum(thick[1,:]))
 ## Get probabilistic
 Xp, Yp, Zp, zallp, ρlow, ρmid, ρhigh, = transD_GP.readxyzrhoϕ(linewanted, 52, 
      pathname="/Users/anray/Documents/work/projects/largeaem/final_01/summaries_all/AusAEM_01/2017_QLD/summary")

zz_portalcurtains/comparelei_2.jl

@@ -0,0 +1,31 @@
+using HiQGA, PyPlot
+linewanted = 2030001
+dr = 100
+xl, yl = [], [-250, 190]
+vmin, vmax = -2.5, 0.5
+linesmoothδ² = 1e-3
+donn = false
+delbin = 600
+## Get LEI
+line, X, Y, Z, thick, σ = transD_GP.readcols([5,7,8,11,[56,85],[26,55]], 
+    "/Users/anray/Documents/work/projects/largeaem/GA_LEI/vtk_WRC_AusAEM_parallel/South/AusAEM_Western_Resources_Corridor_GA_vsum_inversion_South.dat",
+    startfrom=1, decfactor=5)
+idx = line .== linewanted
+Xd, Yd, Zd, thick, σd = map(x->x[idx,:],(X, Y, Z, thick, σ))
+zalld = transD_GP.zboundarytocenter_inexact(cumsum(thick[1,:]), fudgelast=true)
+## Get probabilistic
+Xp, Yp, Zp, zallp, ρlow, ρmid, ρhigh, = transD_GP.readxyzrhoϕ(linewanted, 52, 
+     pathname="/Users/anray/Documents/work/projects/largeaem/final_01/summaries_all/AusAEM_03_WRC/2022_WRC_South/summary/")
+## put them together
+X = [Xd, Matrix(Matrix(Xp')'), Matrix(Matrix(Xp')'), Matrix(Matrix(Xp')')]
+Y = [Yd, Matrix(Matrix(Yp')'), Matrix(Matrix(Yp')'), Matrix(Matrix(Yp')')]
+Z = [Zd, Matrix(Matrix(Zp')'), Matrix(Matrix(Zp')'), Matrix(Matrix(Zp')') ]
+zall = [zalld, zallp, zallp, zallp]
+σ = ([Matrix(log10.(σd')), -ρhigh, -ρmid, -ρlow])
+## plot
+XYprofiles = nothing#[Xp[520];Yp[520]]
+titles = ["Deterministic conductivity", "10th Percentile conductivity", "50th Percentile conductivity", "90th Percentile conductivity"]
+xrd, yrd, axd, zatXY, satXY = transD_GP.plotmanygrids(deepcopy(σ), deepcopy(X), deepcopy(Y), deepcopy(Z), deepcopy(zall); 
+    dr, vmin, vmax, donn, xl, yl, δ²=linesmoothδ², figsize=(20,10), titles, XYprofiles,
+    preferEright=true, plotbinning=false, fontsize=10, delbin, hspace=0.22, spacefactor=0.1)
+savefig("Line_$linewanted.png", dpi=500)