From 31125e98b77e0ad66b0fdf515d1bb417fdec95ee Mon Sep 17 00:00:00 2001 From: Anandaroop Ray Date: Wed, 7 Feb 2024 18:02:44 +1100 Subject: [PATCH 1/6] read xyzrho files and make vtk --- src/CommonToAll.jl | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/CommonToAll.jl b/src/CommonToAll.jl index 5432a240..9227d825 100644 --- a/src/CommonToAll.jl +++ b/src/CommonToAll.jl @@ -1054,6 +1054,10 @@ function writevtkfromsounding(lineofsoundings::Array{S, 1}, zall) where S<:Sound @info("opening summary: Line $(lnum)") rholow, rhomid, rhohigh = map(x->readdlm(x*"_line_$(lnum)_"*"summary.txt"), ["rho_low", "rho_mid", "rho_hi"]) + writevtkfromxyzrho(rholow, rhomid, rhohigh, X, Y, Z, zall) +end + +function writevtkfromxyzrho(rholow, rhomid, rhohigh, X, Y, Z, zall) Ni, Nj = map(x->length(x), (X, Y)) Nk = length(zall) x = [X[i] for i = 1:Ni, j = 1:1, k = 1:Nk] @@ -1069,7 +1073,7 @@ function writevtkfromsounding(lineofsoundings::Array{S, 1}, zall) where S<:Sound vtk["cond_high"] = σhigh end nothing -end +end function writevtkfromsounding(s::Vector{Array{S, 1}}, zall) where S<:Sounding pmap(s) do x From 97b7c9af227a8130937f335d8a2941e0dc53a1e2 Mon Sep 17 00:00:00 2001 From: Anandaroop Ray Date: Thu, 8 Feb 2024 15:49:14 +1100 Subject: [PATCH 2/6] changes for vtk viz --- Project.toml | 2 ++ src/CommonToAll.jl | 15 ++++++++++----- 2 files changed, 12 insertions(+), 5 deletions(-) diff --git a/Project.toml b/Project.toml index 032a03c0..0e586224 100644 --- a/Project.toml +++ b/Project.toml @@ -16,6 +16,7 @@ Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838" FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549" Formatting = "59287772-0a20-5a39-b81b-1366585eb4c0" +Glob = "c27321d9-0574-5035-807b-f59d2c89b15c" ImageDistances = "51556ac3-7006-55f5-8cb3-34580c88182d" ImageIO = "82e4d734-157c-48bb-816b-45c225c6df19" Images = "916415d5-f1e6-5110-898d-aaa5f9f070e0" @@ -55,6 +56,7 @@ Distributions = "0.25.58" FastGaussQuadrature = "0.4.9" FileIO = "1.10" Formatting = "0.4.2" +Glob = "1.3.1" ImageDistances = "0.2.15" ImageIO = "0.5.9" Images = "0.25.2" diff --git a/src/CommonToAll.jl b/src/CommonToAll.jl index 9227d825..0d79eae3 100644 --- a/src/CommonToAll.jl +++ b/src/CommonToAll.jl @@ -1,7 +1,7 @@ module CommonToAll using PyPlot, StatsBase, Statistics, Distances, LinearAlgebra, DelimitedFiles, ..AbstractOperator, NearestNeighbors, Printf, ReadVTK, - KernelDensitySJ, KernelDensity, Interpolations, CSV, WriteVTK, Distributed + KernelDensitySJ, KernelDensity, Interpolations, CSV, WriteVTK, Distributed, Glob import ..Options, ..OptionsStat, ..OptionsNonstat, ..OptionsNuisance, ..history, ..GP.κ, ..calcfstar!, ..AbstractOperator.Sounding, @@ -19,7 +19,7 @@ export trimxft, assembleTat1, gettargtemps, checkns, getchi2forall, nicenup, plo readcols, colstovtk, findclosestidxincolfile, zcentertoboundary, zboundarytocenter, writeijkfromsounding, nanmean, infmean, nanstd, infstd, kde_sj, plotmanygrids, readwell, getlidarheight, plotblockedwellonimages, getdeterministicoutputs, getprobabilisticoutputs, - readfzipped, readxyzrhoϕ, writevtkxmlforcurtain + readfzipped, readxyzrhoϕ, writevtkxmlforcurtain, getprobabilisticlinesfromdirectory, writevtkfromxyzrho # Kernel Density stuff abstract type KDEtype end @@ -1054,10 +1054,15 @@ function writevtkfromsounding(lineofsoundings::Array{S, 1}, zall) where S<:Sound @info("opening summary: Line $(lnum)") rholow, rhomid, rhohigh = map(x->readdlm(x*"_line_$(lnum)_"*"summary.txt"), ["rho_low", "rho_mid", "rho_hi"]) - writevtkfromxyzrho(rholow, rhomid, rhohigh, X, Y, Z, zall) + writevtkfromxyzrho(rholow, rhomid, rhohigh, X, Y, Z, zall, lnum) end -function writevtkfromxyzrho(rholow, rhomid, rhohigh, X, Y, Z, zall) +function getprobabilisticlinesfromdirectory(src_dir) + fn = readdir(glob"rho_low*_summary_xyzrho.txt", src_dir) + [parse(Int, split(split(f, "_summary")[1], "line_")[2]) for f in fn] +end + +function writevtkfromxyzrho(rholow, rhomid, rhohigh, X, Y, Z, zall, lnum; dst_dir="") Ni, Nj = map(x->length(x), (X, Y)) Nk = length(zall) x = [X[i] for i = 1:Ni, j = 1:1, k = 1:Nk] @@ -1067,7 +1072,7 @@ function writevtkfromxyzrho(rholow, rhomid, rhohigh, X, Y, Z, zall) # switch from rho to sigma so low, hi interchanged [-rho[k, i] for i = 1:Ni, j = 1:1, k = 1:Nk] end - vtk_grid("Line_$(lnum)", x, y, z) do vtk + vtk_grid(joinpath(dst_dir, "Line_$(lnum)"), x, y, z) do vtk vtk["cond_low"] = σlow vtk["cond_mid"] = σmid vtk["cond_high"] = σhigh From f1bca9ad649377a8161fd9f7c4830ec6d57eed40 Mon Sep 17 00:00:00 2001 From: Anandaroop Ray Date: Sun, 11 Feb 2024 21:29:23 +1100 Subject: [PATCH 3/6] writing phid to vtk --- src/CommonToAll.jl | 18 +++++++++++++++++- 1 file changed, 17 insertions(+), 1 deletion(-) diff --git a/src/CommonToAll.jl b/src/CommonToAll.jl index 0d79eae3..df0f8275 100644 --- a/src/CommonToAll.jl +++ b/src/CommonToAll.jl @@ -19,7 +19,8 @@ export trimxft, assembleTat1, gettargtemps, checkns, getchi2forall, nicenup, plo readcols, colstovtk, findclosestidxincolfile, zcentertoboundary, zboundarytocenter, writeijkfromsounding, nanmean, infmean, nanstd, infstd, kde_sj, plotmanygrids, readwell, getlidarheight, plotblockedwellonimages, getdeterministicoutputs, getprobabilisticoutputs, - readfzipped, readxyzrhoϕ, writevtkxmlforcurtain, getprobabilisticlinesfromdirectory, writevtkfromxyzrho + readfzipped, readxyzrhoϕ, writevtkxmlforcurtain, getprobabilisticlinesfromdirectory, + writevtkfromxyzrho, writevtkphifromsummary # Kernel Density stuff abstract type KDEtype end @@ -1080,6 +1081,21 @@ function writevtkfromxyzrho(rholow, rhomid, rhohigh, X, Y, Z, zall, lnum; dst_di nothing end +function writevtkphifromsummary(phid, sdev_phid, X, Y, Z, lnum; dst_dir="") + Ni, Nj = map(x->length(x), (X, Y)) + x = [X[i] for i = 1:Ni, j = 1:1, k = 1:1] + y = [Y[i] for i = 1:Ni, j = 1:1, k = 1:1] + z = [Z[i] for i = 1:Ni, j = 1:1, k = 1:1] + ϕmean, ϕsdev = map((phid, sdev_phid)) do p + [p[i] for i = 1:Ni, j = 1:1, k = 1:1] + end + vtk_grid(joinpath(dst_dir, "phid_Line_$(lnum)"), x, y, z) do vtk + vtk["phid_mean"] = ϕmean + vtk["phid_sdev"] = ϕsdev + end + nothing +end + function writevtkfromsounding(s::Vector{Array{S, 1}}, zall) where S<:Sounding pmap(s) do x writevtkfromsounding(x, zall) From 5de740cb780d6ac2640c4e5f98830ef2f51e3359 Mon Sep 17 00:00:00 2001 From: Anand Ray Date: Fri, 23 Feb 2024 12:31:44 +1100 Subject: [PATCH 4/6] changes to VTK XML for GA portal --- src/CommonToAll.jl | 22 +++++++++++++++------- 1 file changed, 15 insertions(+), 7 deletions(-) diff --git a/src/CommonToAll.jl b/src/CommonToAll.jl index df0f8275..549174b5 100644 --- a/src/CommonToAll.jl +++ b/src/CommonToAll.jl @@ -1129,25 +1129,33 @@ function getvtkwholeextent(fname) whole_extent end -function writevtkxmlforcurtain(vtkfname::String; src_epsg=0, dst_epsg=0, suffix="") +function writevtkxmlforcurtain(vtkfname::String; src_epsg=0, dst_epsg=0, suffix="", vmin=0, vmax=0) whole_extent = getvtkwholeextent(vtkfname) nrows = whole_extent[2] ncols = whole_extent[6] - writevtkxmlforcurtain(;vtkfname, nrows, ncols, src_epsg, dst_epsg, suffix) + writevtkxmlforcurtain(;vtkfname, nrows, ncols, src_epsg, dst_epsg, suffix, vmin, vmax) end -function writevtkxmlforcurtain(;vtkfname="", nrows=0, ncols=0, src_epsg=0, dst_epsg=0, suffix="") +function writevtkxmlforcurtain(;vtkfname="", + nrows=0, ncols=0, src_epsg=0, dst_epsg=0, suffix="", vmin=0, vmax=0) # nrows are usually ndepths # ncols are usually nsoundings # 0 based indexing so -1 than actual out = """ - - + 0 $(nrows) 0 0 0 $(ncols) - $(vtkfname*suffix) - EPSG:$src_epsgEPSG:$dst_epsgvtkXmlReader + $(basename(vtkfname)*suffix) + EPSG:$src_epsg + EPSG:$dst_epsg + vtkXmlReader + Turbo + true + true + $vmin + $vmax + """ outname = split(vtkfname,".vts")[1]*".xml" f = open(outname, "w") From 4c6b5d0a4fddaa66ad58857c4a17f37a2be866e8 Mon Sep 17 00:00:00 2001 From: Anand Ray Date: Fri, 23 Feb 2024 12:44:38 +1100 Subject: [PATCH 5/6] portal curtain stuff --- zz_portalcurtains/RDP.jl | 275 ++++++++++++++++++ .../makecurtains_probabilistic_earthsci.jl | 23 ++ zz_portalcurtains/makevtk_from_all_xyzrho.jl | 24 ++ zz_portalcurtains/makevtk_from_xyzrho.jl | 7 + 4 files changed, 329 insertions(+) create mode 100644 zz_portalcurtains/RDP.jl create mode 100644 zz_portalcurtains/makecurtains_probabilistic_earthsci.jl create mode 100644 zz_portalcurtains/makevtk_from_all_xyzrho.jl create mode 100644 zz_portalcurtains/makevtk_from_xyzrho.jl diff --git a/zz_portalcurtains/RDP.jl b/zz_portalcurtains/RDP.jl new file mode 100644 index 00000000..2bf3c7d7 --- /dev/null +++ b/zz_portalcurtains/RDP.jl @@ -0,0 +1,275 @@ +module RDP +using LinearAlgebra, Dates, ArchGDAL, Printf, + PyPlot, Images, FileIO, HiQGA +import GeoFormatTypes as GFT +const mpl = PyPlot.matplotlib +const tilesize = 512 +function perpdist(a, b, p) + # a and b are two points on the line in space + # p is the point to compute distance of from line + n̂ = normalize(b-a) + dperp = p - a - dot(p-a, n̂)n̂ + norm(dperp) +end + +function rdpreduce(pointlist, ϵ=1.0) + dists = perpdist.(Ref(pointlist[1]), Ref(pointlist[end]), pointlist) + d, idx = findmax(dists) + if d > ϵ + # split line at the farthest point, doing RDP on left and right segments + left = rdpreduce(pointlist[1:idx], ϵ) + right = rdpreduce(pointlist[idx:end], ϵ) + # concatenate results of RDP on segments + keep = vcat(left[1:end-1], right) + else + # don't keep the most distant point from the line, only the end points of the line + keep = [pointlist[1], pointlist[end]] + end + keep +end + +function worldcoordinates(;gridr=nothing, gridz=nothing) + # from https://support.esri.com/en-us/knowledge-base/faq-what-is-the-format-of-the-world-file-used-for-geore-000002860 + A = step(gridr) + D = 0 + B = 0 + E = step(gridz) + C = A/2 + F = gridz[1] + E/2 + [A, D, B, E, C, F] +end + +makeplist(X,Y) = [[x,y] for (x,y) in zip(X,Y)] + +reprojecttoGDA94(plist, fromepsg) = ArchGDAL.reproject(plist, GFT.EPSG(fromepsg), GFT.EPSG(4283) ) + +makexyfromlatlonglist(latlonglist) = [[l[i] for l in latlonglist] for i in 2:-1:1] + +function writeworldXML(latlonglist; dst_dir= "", line="", suffix="", fnamebar="colorbar.jpg", gridz=nothing) + img = load(joinpath(dst_dir, "jpeg", "$line"*suffix*".jpg")) + h, w = size(img) + ntilelevels = gettilelevels(w,h,tilesize) + firstpart = + """ + + + $(string(line)*suffix) + $fnamebar + + ./ + + + LocalRequester + TransparentColorTransformer(255,255,255,0.2) + ResizeTransformer($tilesize,$tilesize) + + $(now()) + $(string(line)*suffix) + AWS_CACHE_NAME/$(string(line)*suffix) + image/jpg + .jpg + + image/jpg + + + + + + + + + + """ + middlepart = writelatlonglist(latlonglist) + lastpart = + """ + + $(gridz[1] + step(gridz)/2) + $(gridz[end] - step(gridz)/2) + false + 2 + true + true + image/dds + true + 0.5 + + """ + xmlpath = joinpath(dst_dir,"xml_and_tiles") + isdir(xmlpath) || mkpath(xmlpath) + f = open(joinpath(xmlpath, "$line"*suffix*".xml"), "w") + map([firstpart, middlepart, lastpart]) do part + write(f, part) + end + close(f) +end + +function writelatlonglist(latlonglist) + s = "" + for p in latlonglist + if p != last(latlonglist) + s*=@sprintf("\t\t\n", p[1], p[2]) + else + s*=@sprintf("\t\t", p[1], p[2]) + end + end + s +end + +function writedosbatchfile(jpegsavename::String; isfirst=false, islast=false) + lname = basename(jpegsavename) + dst_dir = dirname(dirname(jpegsavename)) + fname = joinpath(dst_dir, "tileall.bat") + if isfirst + f = open(fname, "w") + write(f, "@echo off\n") + else + f = open(fname, "a") + end + write(f, "call ribbon.bat -tilesize $tilesize -noLayerDef -source jpeg\\$lname -output xml_and_tiles\\\n") + islast && write(f, "pause\n") + close(f) + nothing +end + +function pltpoint(ax, pts) + pt = reduce(hcat, pts)' + ax.plot(pt[:,1], pt[:,2], "*-") +end + +function writeimageandcolorbar(img::Array, gridr, gridz, line::Int; cmap="turbo", dpi=300, dst_dir="", + vmin=-Inf, vmax=Inf, fnamebar="colorbar.jpg", barfigsize=(1,6), suffix="", isfirst=false, islast=false, + writecolorbar=false, VE=20, shrink=25_000) + if (isinf(vmin) || isinf(vmax)) + vmin, vmax = extrema(filter(!isnan, img)) + end + jpegpath = joinpath(dst_dir,"jpeg") + if writecolorbar + fig, ax = plt.subplots(1, 2, gridspec_kw=Dict("width_ratios" => [0.3,1]), figsize=barfigsize) + norm = mpl.colors.Normalize(vmin, vmax) + cbar = fig.colorbar(mpl.cm.ScalarMappable(;norm, cmap), cax=ax[1], orientation="vertical") + cbar.ax.tick_params(labelsize=4, pad=0) + cbar.set_label("Log₁₀ S/m", fontsize=4, labelpad=0) + ax[2].axis("off") + fig.tight_layout() + isdir(jpegpath) || mkpath(jpegpath) + xmlpath = joinpath(dst_dir,"xml_and_tiles") + savefig(joinpath(jpegpath, fnamebar); dpi) + isdir(xmlpath) || mkpath(xmlpath) + savefig(joinpath(xmlpath, fnamebar); dpi) + close(fig) + end + figsize = gridr[end]/shrink, abs(diff([extrema(gridz)...])[1])*VE/shrink + fig, ax = plt.subplots(;figsize) + ax.imshow(img; cmap, extent=[gridr[1], gridr[end], gridz[end], gridz[1]]) + ax.set_aspect(VE) + ax.set_aspect("auto") + ax.axis("off") + jpegsavename = joinpath(jpegpath, "$line"*suffix*".jpg") + savefig(jpegsavename; dpi, bbox_inches = "tight", pad_inches = 0) + writedosbatchfile(jpegsavename; isfirst, islast) + close(fig) + nothing +end + +function gettilelevels(w,h,tilesize) + xc = w/tilesize + yc = h/tilesize + levels = 0 + + while (4 * xc * yc >= 1) + levels+=1; + xc /= 2. + yc /= 2. + end + levels +end + +function doonecurtaintriad(line::Int; nlayers=0, pathname="", dr=0, dz=0, dst_dir="", + donn=false, ϵfrac=0, src_epsg=0, barfigsize=(0.2, 1.2), isfirst=false, islast=false, + dpi=400, cmap="turbo", fnamebar="colorbar.jpg", vmin=-Inf, vmax=Inf, + writecolorbar=true, VE=20, shrink = 25_000) + X, Y, Z, zall, ρlow, ρmid, ρhigh, ρavg, ϕmean, ϕsdev = + transD_GP.readxyzrhoϕ(line, nlayers; pathname) + plist = makeplist(X,Y) + R = transD_GP.CommonToAll.cumulativelinedist(X,Y) + pgood = rdpreduce(plist, ϵfrac*R[end]) + latlonglist = RDP.reprojecttoGDA94(pgood, src_epsg) + map(zip([-ρhigh, -ρmid, -ρlow], ["_low", "_mid", "_high"],(1:3))) do (σ, suffix, i) + img, gridr, gridz, topofine, R = transD_GP.makegrid(σ, X, Y, Z; donn, + dr, zall, dz) + writecb = false + if writecolorbar == true + i == 1 && (writecb = true) + end + isl = false + if islast == true + i == 3 && (isl = true) + end + writeimageandcolorbar(img, gridr, gridz, line; cmap, fnamebar, suffix, dst_dir, isfirst=writecb, islast=isl, + vmin, vmax, dpi, writecolorbar=writecb, + barfigsize, VE, shrink) + writeworldXML(latlonglist;line, gridz, suffix, dst_dir) + end +end + +function writebunchxmlfile(lines::Vector{Int}, dst_dir::String; prefix="", suffix="") + xmlpath = joinpath(dst_dir,"xml_and_tiles") + bunchfname = prefix*suffix + f = open(joinpath(xmlpath, bunchfname*".xml"), "w") + out = + """ + + + """ + write(f, out) + for line in lines + write(f,"\t\n") + end + out = + """ + + + """ + write(f, out) + close(f) + nothing +end + +function doallcurtaintriads(;src_dir="", dst_dir="curtains", prefix="", + nlayers=0, dr=0, dz=0, + donn=false, ϵfrac=0, src_epsg=0, barfigsize=(0.2, 1.2), + dpi=400, cmap="turbo", fnamebar="colorbar.jpg", vmin=-Inf, vmax=Inf, + VE=20, shrink = 25_000) + + isdir(dst_dir) || mkpath(dst_dir) + lines = transD_GP.getprobabilisticlinesfromdirectory(src_dir) + map(lines) do line + @info "Doing line $line" + writecb = line == first(lines) ? true : false + isfirst = writecb + islast = line == last(lines) ? true : false + doonecurtaintriad(line; nlayers, pathname=src_dir, dr, dz, ϵfrac, src_epsg, isfirst, islast, + barfigsize, dpi, cmap, fnamebar, writecolorbar=writecb, dst_dir, shrink, VE, + vmin, vmax) + end + for suffix in ("low", "mid", "high") + writebunchxmlfile(lines, dst_dir; prefix, suffix) + end + nothing +end + +function writevtkfromxyzrhodir(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 = reprojecttoGDA94(plist, src_epsg) + X, Y = makexyfromlatlonglist(latlonglist) + transD_GP.writevtkfromxyzrho(ρlow, ρmid, ρhigh, X, Y, Z, zall, ln; dst_dir) + transD_GP.writevtkphifromsummary(ϕmean, ϕsdev, X, Y, Z, ln; dst_dir) + end +end + +end # module \ No newline at end of file diff --git a/zz_portalcurtains/makecurtains_probabilistic_earthsci.jl b/zz_portalcurtains/makecurtains_probabilistic_earthsci.jl new file mode 100644 index 00000000..d1c93175 --- /dev/null +++ b/zz_portalcurtains/makecurtains_probabilistic_earthsci.jl @@ -0,0 +1,23 @@ +using HiQGA, PyPlot +includet("RDP.jl") +## +nlayers = 52 +dr = 480 +dz = 2 +dpi = 400 +fnamebar = "colorbar.jpg" +cmap="turbo" +vmin, vmax = -0.5, 2.5 +src_epsg = 28354 +ϵfrac = 0.001 +barfigsize = (0.4, 1.2) +shrink = 8000 +VE = 30 +prefix = "ERC_01" +## multiple lines +# destdir for multiple lines +dst_dir = "/Users/anray/Documents/work/projects/curtainstuff/ERC_01_curtains_for_mgn" +src_dir = "/Users/anray/Documents/work/projects/curtainstuff/summaries" +RDP.doallcurtaintriads.(;src_dir, dst_dir, nlayers, dr, dz, ϵfrac, src_epsg, + barfigsize, dpi, cmap, fnamebar, shrink, VE, prefix, + vmin, vmax) \ No newline at end of file diff --git a/zz_portalcurtains/makevtk_from_all_xyzrho.jl b/zz_portalcurtains/makevtk_from_all_xyzrho.jl new file mode 100644 index 00000000..e97cb320 --- /dev/null +++ b/zz_portalcurtains/makevtk_from_all_xyzrho.jl @@ -0,0 +1,24 @@ +using HiQGA, PyPlot +includet("RDP.jl") +nlayers = 52 +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] +map(zip(src_dir, src_epsg)) do (dir, epsg) + dst_dir = joinpath(dirname(dir),"GDA94_vtk") + RDP.writevtkfromxyzrhodir(nlayers; src_dir=dir, dst_dir, src_epsg=epsg) +end +## XML files +using HiQGA, Glob +idx_summary = [basename(it[1]) == "GDA94_vtk" for it in items] +src_dir = [it[1] for it in items[idx_summary]] +vmin, vmax = -2.5, 0.5 +map(src_dir) do (dir) + fn = readdir(glob"Line*.vts", dir) + # @info fn[1] + # GDA94 to WGS84 + map(fn) do f + transD_GP.writevtkxmlforcurtain(f; src_epsg=4283, dst_epsg=4326, suffix="", vmin, vmax) + end +end \ No newline at end of file diff --git a/zz_portalcurtains/makevtk_from_xyzrho.jl b/zz_portalcurtains/makevtk_from_xyzrho.jl new file mode 100644 index 00000000..b0f01036 --- /dev/null +++ b/zz_portalcurtains/makevtk_from_xyzrho.jl @@ -0,0 +1,7 @@ +using HiQGA, PyPlot +includet("RDP.jl") +src_epsg = 28354 +nlayers = 52 +src_dir = "/Users/anray/Documents/work/projects/curtainstuff/summaries" +dst_dir = "/Users/anray/Documents/work/projects/curtainstuff/ERC_01_vtk_GDA94" +RDP.writevtkfromxyzrhodir(nlayers; src_dir, dst_dir, src_epsg) \ No newline at end of file From 1bd4c1dca97974e08c9f9fd2f816938719c46cf3 Mon Sep 17 00:00:00 2001 From: Anandaroop Ray Date: Fri, 23 Feb 2024 14:10:41 +1100 Subject: [PATCH 6/6] whole dir xml --- src/CommonToAll.jl | 23 ++++++++++---------- zz_portalcurtains/RDP.jl | 9 ++++++-- zz_portalcurtains/makevtk_from_all_xyzrho.jl | 18 +++------------ 3 files changed, 22 insertions(+), 28 deletions(-) diff --git a/src/CommonToAll.jl b/src/CommonToAll.jl index 549174b5..52b65732 100644 --- a/src/CommonToAll.jl +++ b/src/CommonToAll.jl @@ -1144,17 +1144,18 @@ function writevtkxmlforcurtain(;vtkfname="", out = """ - - 0 $(nrows) 0 0 0 $(ncols) - $(basename(vtkfname)*suffix) - EPSG:$src_epsg - EPSG:$dst_epsg - vtkXmlReader - Turbo - true - true - $vmin - $vmax + + 0 $(nrows) 0 0 0 $(ncols) + $(basename(vtkfname)*suffix) + EPSG:$src_epsg + EPSG:$dst_epsg + vtkXmlReader + Turbo + true + true + $vmin + $vmax + Log10 S/m """ outname = split(vtkfname,".vts")[1]*".xml" diff --git a/zz_portalcurtains/RDP.jl b/zz_portalcurtains/RDP.jl index 2bf3c7d7..dcdcea1d 100644 --- a/zz_portalcurtains/RDP.jl +++ b/zz_portalcurtains/RDP.jl @@ -4,6 +4,9 @@ using LinearAlgebra, Dates, ArchGDAL, Printf, import GeoFormatTypes as GFT const mpl = PyPlot.matplotlib const tilesize = 512 +const epsg_GDA94 = 4283 +const epsg_WGS84 = 4326 + function perpdist(a, b, p) # a and b are two points on the line in space # p is the point to compute distance of from line @@ -41,7 +44,7 @@ end makeplist(X,Y) = [[x,y] for (x,y) in zip(X,Y)] -reprojecttoGDA94(plist, fromepsg) = ArchGDAL.reproject(plist, GFT.EPSG(fromepsg), GFT.EPSG(4283) ) +reprojecttoGDA94(plist, fromepsg) = ArchGDAL.reproject(plist, GFT.EPSG(fromepsg), GFT.EPSG(epsg_GDA94) ) makexyfromlatlonglist(latlonglist) = [[l[i] for l in latlonglist] for i in 2:-1:1] @@ -259,7 +262,7 @@ function doallcurtaintriads(;src_dir="", dst_dir="curtains", prefix="", nothing end -function writevtkfromxyzrhodir(nlayers::Int; src_dir="", dst_dir="", src_epsg=0) +function writevtkfromxyzrhodir(nlayers::Int; src_dir="", dst_dir="", src_epsg=0, vmin=0, vmax=0) lines = transD_GP.getprobabilisticlinesfromdirectory(src_dir) isdir(dst_dir) || mkpath(dst_dir) map(lines) do ln @@ -269,6 +272,8 @@ function writevtkfromxyzrhodir(nlayers::Int; src_dir="", dst_dir="", src_epsg=0) X, Y = makexyfromlatlonglist(latlonglist) transD_GP.writevtkfromxyzrho(ρlow, ρmid, ρhigh, X, Y, Z, zall, ln; dst_dir) transD_GP.writevtkphifromsummary(ϕmean, ϕsdev, X, Y, Z, ln; dst_dir) + fn = joinpath(dst_dir, "Line_$(ln).vts") + transD_GP.writevtkxmlforcurtain(fn; src_epsg=epsg_GDA94, dst_epsg=epsg_WGS84, suffix="", vmin, vmax) end end diff --git a/zz_portalcurtains/makevtk_from_all_xyzrho.jl b/zz_portalcurtains/makevtk_from_all_xyzrho.jl index e97cb320..9573cddf 100644 --- a/zz_portalcurtains/makevtk_from_all_xyzrho.jl +++ b/zz_portalcurtains/makevtk_from_all_xyzrho.jl @@ -1,24 +1,12 @@ using HiQGA, PyPlot 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] map(zip(src_dir, src_epsg)) do (dir, epsg) dst_dir = joinpath(dirname(dir),"GDA94_vtk") - RDP.writevtkfromxyzrhodir(nlayers; src_dir=dir, dst_dir, src_epsg=epsg) -end -## XML files -using HiQGA, Glob -idx_summary = [basename(it[1]) == "GDA94_vtk" for it in items] -src_dir = [it[1] for it in items[idx_summary]] -vmin, vmax = -2.5, 0.5 -map(src_dir) do (dir) - fn = readdir(glob"Line*.vts", dir) - # @info fn[1] - # GDA94 to WGS84 - map(fn) do f - transD_GP.writevtkxmlforcurtain(f; src_epsg=4283, dst_epsg=4326, suffix="", vmin, vmax) - end -end \ No newline at end of file + RDP.writevtkfromxyzrhodir(nlayers; src_dir=dir, dst_dir, src_epsg=epsg, vmin, vmax) +end \ No newline at end of file