Merge pull request #70 from GeoscienceAustralia/curtainstuff

Curtainstuff

zz_portalcurtains/RDP.jl

@@ -2,6 +2,7 @@ module RDP
 using LinearAlgebra, Dates, ArchGDAL, Printf, 
     PyPlot, Images, FileIO, HiQGA, Interpolations
 import GeoFormatTypes as GFT
+import GeoDataFrames as GDF
 const mpl = PyPlot.matplotlib
 const tilesize = 512
 const epsg_GDA94 = 4283
@@ -311,6 +312,25 @@ function writevtkfromxyzrhodir(nlayers::Int; src_dir="", dst_dir="", src_epsg=0,
     end
 end
 
+function writeshpfromxyzrhodir(nlayers::Int; prefix="", src_dir="", dst_dir="", src_epsg=0, writeesri=false)
+    lines = transD_GP.getprobabilisticlinesfromdirectory(src_dir)
+    isdir(dst_dir) || mkpath(dst_dir)
+    table = map(lines) do ln
+        X, Y, _ = transD_GP.readxyzrhoϕ(ln, nlayers; pathname=src_dir)
+        plist = makeplist(X, Y)
+        longlat = reverse.(reprojecttoGDA94(plist, src_epsg))
+        R = transD_GP.CommonToAll.cumulativelinedist(X,Y)
+        (;geom=ArchGDAL.createlinestring(longlat), Line=ln, Length=R[end], soundings_inverted=length(R))
+    end
+    if writeesri
+        fn = joinpath(dst_dir, prefix*".shp")
+        GDF.write(fn, table; geom_columns=(:geom,), crs=GFT.EPSG(epsg_GDA94))
+    else    
+        fn = joinpath(dst_dir, prefix*".gpkg")
+        GDF.write(fn, table; geom_columns=(:geom,), crs=GFT.EPSG(epsg_GDA94))
+    end    
+end
+
 function writeaseggdffromxyzrho(nlayers::Int; src_dir="", dst_dir="", 
          fname="", src_epsg=0)
     isdir(dst_dir) || mkpath(dst_dir)

zz_portalcurtains/makeshp_from_all_xyzrho.jl

@@ -0,0 +1,14 @@
+using HiQGA, PyPlot
+includet("RDP.jl")
+nlayers = 52
+rootdir = "/Users/anray/Documents/work/projects/largeaem/final_01/shapefiles/"
+## multiple surveys
+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 (sdir, epsg)
+    parts = split(sdir,"/")
+    prefix = parts[end-2]*"_"*parts[end-1]
+    RDP.writeshpfromxyzrhodir(nlayers; prefix, src_dir=sdir, dst_dir=rootdir, src_epsg=epsg)
+end