colstvtk with longlat possible

GeoscienceAustralia · Jul 3, 2024 · 5a680c5 · 5a680c5
1 parent 83c5036
commit 5a680c5
Show file tree

Hide file tree

Showing 2 changed files with 16 additions and 2 deletions.
diff --git a/src/CommonToAll.jl b/src/CommonToAll.jl
@@ -1225,6 +1225,10 @@ end
 function colstovtk(cols::Dict, fname::String; decfactor=1, hasthick=true, islog10=false, prefix="")
     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, lines, fname::String; hasthick=true, islog10=false, prefix="")
     linenos  = unique(Int.(lines))
     thick = thick[1,:] # assumes all thicknesses are same
     zall = thicktodepth(thick; hasthick)
@@ -1236,7 +1240,7 @@ function colstovtk(cols::Dict, fname::String; decfactor=1, hasthick=true, islog1
         fpath = joinpath(dstring, "LEI_Line_"*prefix*fstring)
         colstovtk(X[idx], Y[idx], Z[idx], σ[idx,:], zall, fpath, islog10)
     end
-end    
+end
 
 function thicktodepth(thick; hasthick=true)
     if hasthick

diff --git a/zz_portalcurtains/RDP.jl b/zz_portalcurtains/RDP.jl
@@ -3,7 +3,7 @@ using LinearAlgebra, Dates, ArchGDAL, Printf, DataFrames, CSV,
     PyPlot, Images, FileIO, HiQGA, Interpolations, NearestNeighbors, DelimitedFiles
 import GeoFormatTypes as GFT
 import GeoDataFrames as GDF
-import HiQGA.transD_GP.LineRegression.getsmoothline
+import HiQGA.transD_GP.LineRegression.getsmoothline, HiQGA.transD_GP.CommonToAll.colstovtk
 
 const mpl = PyPlot.matplotlib
 const tilesize = 512
@@ -80,6 +80,16 @@ reprojecttoGDA94(plist, fromepsg) = ArchGDAL.reproject(plist, GFT.EPSG(fromepsg)
 
 makexyfromlatlonglist(latlonglist) = [[l[i] for l in latlonglist] for i in 2:-1:1]
 
+function colstovtk(cols::Dict, fname::String, src_epsg::Int; decfactor=1, hasthick=true, islog10=false, prefix="")
+    # outputs in long, lat
+    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)
+    plist = makeplist(X, Y)
+    latlonglist = reprojecttoGDA94(plist, src_epsg)
+    long, lat = makexyfromlatlonglist(latlonglist)
+    colstovtk(long, lat, Z, σ, thick, lines, fname; hasthick, islog10, prefix)
+end
+
 function getimgsize(line, dst_dir, suffix)
     img = load(joinpath(dst_dir, "jpeg", "$line"*suffix*".jpg"))
     h, w = size(img)