segy for multiple lines from one ASEG-GDF

src/CommonToAll.jl

@@ -20,7 +20,7 @@ export trimxft, assembleTat1, gettargtemps, checkns, getchi2forall, nicenup, plo
         writeijkfromsounding, nanmean, infmean, nanstd, infstd, infnanmean, infnanstd, 
         kde_sj, plotmanygrids, readwell, getlidarheight, plotblockedwellonimages, getdeterministicoutputs, writeasegdfnfromonerow,
         writeasegdat, getprobabilisticoutputs, readfzipped, readxyzrhoϕ, writevtkxmlforcurtain, getRandgridr, getallxyinr, getXYlast,
-        getprobabilisticlinesfromdirectory, readxyzrhoϕnu, plotgausshist, findMAE, getclosestidx, getpostatXY, getpostatwell
+        getprobabilisticlinesfromdirectory, readxyzrhoϕnu, plotgausshist, findMAE, getclosestidx, getpostatXY, getpostatwell, thicktodepth
 
 # Kernel Density stuff
 abstract type KDEtype end
@@ -1206,7 +1206,8 @@ function writevtkxmlforcurtain(;vtkfname="",
 end
 
 function colstovtk(cols::Vector, fname::String; decfactor=1, hasthick=true, islog10=false, prefix="")
-    # for one file in a direcoty
+    # for one file with one line, not used much I don't think 
+    # TBD: delete function
     X, Y, Z, σ, thick = readcols(cols, fname; decfactor)
     thick = thick[1,:]
     zall = thicktodepth(thick; hasthick)
@@ -1217,6 +1218,7 @@ function colstovtk(cols::Vector, fname::String; decfactor=1, hasthick=true, islo
 end
 
 function colstovtk(X, Y, Z, σ, zall, fpath::String, islog10::Bool)
+    # lowest level call
     Ni = length(X)
     Nk = length(zall)
     x = [X[i] for i = 1:Ni, j = 1:1, k = 1:Nk]
@@ -1231,12 +1233,14 @@ function colstovtk(X, Y, Z, σ, zall, fpath::String, islog10::Bool)
 end
 
 function colstovtk(cols::Dict, fname::String; decfactor=1, hasthick=true, islog10=false, prefix="")
+    # dictionary call for a file with multiple lines
     Xc, Yc, Zc, σc, thickc, linec = map(x->get(cols, x, 0), ["X", "Y", "Z", "cond", "thick", "line"])
     X, Y, Z, σ, thick, lines = readcols([Xc, Yc, Zc, σc, thickc, linec], fname; decfactor)
     colstovtk(X, Y, Z, σ, thick, lines, fname; hasthick, islog10, prefix)
 end    
 
 function colstovtk(X, Y, Z, σ, thick_in, lines, fname::String; hasthick=true, islog10=false, prefix="")
+    # no dictionary version next level call
     linenos  = unique(Int.(lines))
     if ndims(thick_in) == 2 
         thick = thick_in[1,:] # assumes all thicknesses are same

zz_portalcurtains/RDP.jl

@@ -70,9 +70,9 @@ function XYZ_zmid_gridtoSEGY(σ, X, Y, Z; dr=nothing, zall=nothing, dz=nothing,
             dr, zall, dz)
     img[isnan.(img)] .= nanval       
     segypath = dst_dir
-    xymid, = transD_GP.getallxyinr(X, Y, dr)
+    xymid, rm = transD_GP.getallxyinr(X, Y, dr)
     xm, ym, = map(i->xymid[i,:], 1:2)
-    rm = transD_GP.CommonToAll.cumulativelinedist(xm, ym)
+    # rm = transD_GP.CommonToAll.cumulativelinedist(xm, ym)
     topom = (interpolate((gridr,), topofine, Gridded(Linear())))(rm)
     block = SeisBlock(Float32.(img))
     set_header!(block, :dt, dz*1000) # ms
@@ -107,6 +107,40 @@ function writesegyfromxyzrhodir(nlayers::Int; src_dir="", src_epsg=0, dst_dir=""
     end
 end
 
+function colstosegy(cols::Dict, fname::String; dr=nothing, dz=nothing, decfactor=1, hasthick=true, islog10=false, suffix="", dst_dir="", src_epsg=0)
+    # dictionary call for a file with multiple lines
+    @assert !isnothing(dr)
+    @assert !isnothing(dz)
+    Xc, Yc, Zc, σc, thickc, linec = map(x->get(cols, x, 0), ["X", "Y", "Z", "cond", "thick", "line"])
+    X, Y, Z, σ, thick, lines = transD_GP.readcols([Xc, Yc, Zc, σc, thickc, linec], fname; decfactor)
+    colstosegy(X, Y, Z, σ, thick, lines; dr, dz, hasthick, islog10, suffix, dst_dir, src_epsg)
+end    
+
+function colstosegy(X, Y, Z, σ, thick_in, lines; dr=nothing, dz=nothing, hasthick=true, islog10=false, suffix="", dst_dir="", src_epsg=0)
+    # no dictionary version next level call with decimation already done, if any
+    @assert !isnothing(dr)
+    @assert !isnothing(dz)
+    linenos  = unique(Int.(lines))
+    if ndims(thick_in) == 2 
+        thick = thick_in[1,:] # assumes all thicknesses are same
+    else
+        thick = thick_in
+    end
+    zall = transD_GP.thicktodepth(thick; hasthick)
+    isdir(dst_dir) || mkpath(dst_dir)
+    ioproj = open(joinpath(dst_dir, "0000_projection.txt"), "w")
+    write(ioproj, "EPSG: $src_epsg")
+    close(ioproj)
+    for l in linenos
+        idx = lines .== l
+        fstring = "LEI_Line_$l"
+        @info "doing Line "*fstring
+        T = islog10 ?  x->x : x->log10(x) # if log 10 leave alone
+        # [@info size(x) for x in (T.(σ[idx,:]), X[idx], Y[idx], Z[idx], zall)] 
+        XYZ_zmid_gridtoSEGY(T.(σ[idx,:])', X[idx], Y[idx], Z[idx]; dr, zall, dz, dst_dir, fname=fstring, suffix)
+    end
+end
+
 function writepathextentfiles(line, nrows, ncols, gridr, topo, gridz, X, Y; suffix="", dst_dir="", src_epsg=0)
     geompath = joinpath(dst_dir,"geometry")
     isdir(geompath) || mkpath(geompath)