diff --git a/src/CommonToAll.jl b/src/CommonToAll.jl index 9823878..952abf9 100644 --- a/src/CommonToAll.jl +++ b/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 diff --git a/zz_portalcurtains/RDP.jl b/zz_portalcurtains/RDP.jl index e3e39d6..c7794ab 100644 --- a/zz_portalcurtains/RDP.jl +++ b/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)