From aa59ae2e88f8c02a74468a090e04ca5c53c9feb9 Mon Sep 17 00:00:00 2001 From: Boris Kaus Date: Sat, 2 Mar 2024 19:42:00 +0100 Subject: [PATCH] create surfaces function; addition and subtraction included; rename DrapeOnTopo -> drape_on_topo; add tests --- docs/src/man/Tutorial_Jura.md | 2 +- docs/src/man/tools.md | 2 +- src/GeophysicalModelGenerator.jl | 1 + src/data_types.jl | 26 ++++- src/surface_functions.jl | 160 ++++++++++++++++++++++++++++ src/utils.jl | 174 +------------------------------ test/runtests.jl | 3 + tutorials/Tutorial_Jura.jl | 2 +- 8 files changed, 192 insertions(+), 178 deletions(-) create mode 100644 src/surface_functions.jl diff --git a/docs/src/man/Tutorial_Jura.md b/docs/src/man/Tutorial_Jura.md index 43beca0a..2fb34182 100644 --- a/docs/src/man/Tutorial_Jura.md +++ b/docs/src/man/Tutorial_Jura.md @@ -47,7 +47,7 @@ Geology = Screenshot_To_GeoData("SchoriM_Encl_01_Jura-map_A1.png", lowerleft, u You can "drape" this image on the topographic map with ```julia -TopoGeology = DrapeOnTopo(Topo, Geology) +TopoGeology = drape_on_topo(Topo, Geology) ``` ```julia diff --git a/docs/src/man/tools.md b/docs/src/man/tools.md index b6791937..8ab34382 100644 --- a/docs/src/man/tools.md +++ b/docs/src/man/tools.md @@ -18,7 +18,7 @@ GeophysicalModelGenerator.PointData2NearestGrid GeophysicalModelGenerator.Convert2UTMzone GeophysicalModelGenerator.Convert2CartData GeophysicalModelGenerator.ProjectCartData -GeophysicalModelGenerator.DrapeOnTopo +GeophysicalModelGenerator.drape_on_topo GeophysicalModelGenerator.LithostaticPressure! GeophysicalModelGenerator.CountMap ``` diff --git a/src/GeophysicalModelGenerator.jl b/src/GeophysicalModelGenerator.jl index 96e5cf6e..4f4a7bce 100644 --- a/src/GeophysicalModelGenerator.jl +++ b/src/GeophysicalModelGenerator.jl @@ -43,6 +43,7 @@ include("stl.jl") include("ProfileProcessing.jl") include("IO.jl") include("event_counts.jl") +include("surface_functions.jl") # Add optional routines (only activated when the packages are loaded) diff --git a/src/data_types.jl b/src/data_types.jl index f783f291..768e9e13 100644 --- a/src/data_types.jl +++ b/src/data_types.jl @@ -1,7 +1,7 @@ # This is data_types.jl # contains type definitions to be used in GeophysicalModelGenerator -import Base: show +import Base: show, size export GeoData, ParaviewData, UTMData, CartData, LonLatDepthGrid, XYZGrid, Velocity_SphericalToCartesian!, @@ -216,6 +216,7 @@ struct GeoData <: AbstractGeneralGrid end end +size(d::GeoData) = size(d.lon.val) # Print an overview of the Geodata struct: function Base.show(io::IO, d::GeoData) @@ -247,6 +248,23 @@ function Base.show(io::IO, d::GeoData) end +""" + GeoData(lld::Tuple{Array,Array,Array}) + +This creates a `GeoData` struct if you have a Tuple with 3D coordinates as input. +# Example +```julia +julia> data = GeoData(LonLatDepthGrid(-10:10,-5:5,0)) +GeoData + size : (21, 11, 1) + lon ϵ [ -10.0 : 10.0] + lat ϵ [ -5.0 : 5.0] + depth ϵ [ 0.0 : 0.0] + fields : (:Z,) +``` +""" +GeoData(lld::Tuple) = GeoData(lld[1],lld[2],lld[3],(Z=lld[3],)) + """ ParaviewData(x::GeoUnit, y::GeoUnit, z::GeoUnit, values::NamedTuple) @@ -265,6 +283,7 @@ mutable struct ParaviewData <: AbstractGeneralGrid z :: GeoUnit fields :: NamedTuple end +size(d::ParaviewData) = size(d.x.val) # Print an overview of the ParaviewData struct: function Base.show(io::IO, d::ParaviewData) @@ -464,6 +483,7 @@ struct UTMData <: AbstractGeneralGrid end end +size(d::UTMData) = size(d.EW.val) # Print an overview of the UTMData struct: function Base.show(io::IO, d::UTMData) @@ -748,6 +768,7 @@ struct CartData <: AbstractGeneralGrid end end +size(d::CartData) = size(d.x.val) # Print an overview of the UTMData struct: function Base.show(io::IO, d::CartData) @@ -1032,6 +1053,7 @@ struct CartGrid{FT, D} <: AbstractGeneralGrid coord1D :: NTuple{D,StepRangeLen{FT}} # Tuple with 1D vectors in all directions coord1D_cen :: NTuple{D,StepRangeLen{FT}} # Tuple with 1D vectors of center points in all directions end +size(d::CartGrid) = d.N """ @@ -1149,8 +1171,6 @@ function CreateCartGrid(; end - - # view grid object function show(io::IO, g::CartGrid{FT, DIM}) where {FT, DIM} diff --git a/src/surface_functions.jl b/src/surface_functions.jl new file mode 100644 index 00000000..2c228f97 --- /dev/null +++ b/src/surface_functions.jl @@ -0,0 +1,160 @@ +# This contains a number of routines that deal with surfaces +using NearestNeighbors + +export remove_NaN_Surface!, drape_on_topo, is_surface +import Base: +,- + +""" + issurf = is_surface(surf::AbstractGeneralGrid) + +Returns true if `surf` is a horizontal 3D surface. +""" +function is_surface(surf::AbstractGeneralGrid) + if size(surf)[3] == 1 + issurf = true + else + issurf = false + end + return issurf +end + + +function +(a::_T, b::_T) where _T<:AbstractGeneralGrid + @assert size(a) == size(b) + return _addSurfaces(a,b) +end + +function -(a::_T, b::_T) where _T<:AbstractGeneralGrid + @assert size(a) == size(b) + return _subtractSurfaces(a,b) +end + +# Internal routines +_addSurfaces(a::_T, b::_T) where _T<:GeoData = GeoData(a.lon.val, a.lat.val, a.depth.val + b.depth.val, merge(a.fields,b.fields)) +_addSurfaces(a::_T, b::_T) where _T<:UTMData = UTMData(a.EW.val, a.NS.val, a.depth.val + b.depth.val, merge(a.fields,b.fields)) +_addSurfaces(a::_T, b::_T) where _T<:CartData = CartData(a.x.val, a.y.val, a.z.val + b.z.val, merge(a.fields,b.fields)) +_addSurfaces(a::_T, b::_T) where _T<:ParaviewData = ParaviewData(a.x.val, a.y.val, a.z.val + b.z.val, merge(a.fields,b.fields)) + +_subtractSurfaces(a::_T, b::_T) where _T<:GeoData = GeoData(a.lon.val, a.lat.val, a.depth.val - b.depth.val, merge(a.fields,b.fields)) +_subtractSurfaces(a::_T, b::_T) where _T<:UTMData = UTMData(a.EW.val, a.NS.val, a.depth.val - b.depth.val, merge(a.fields,b.fields)) +_subtractSurfaces(a::_T, b::_T) where _T<:CartData = CartData(a.x.val, a.y.val, a.z.val - b.z.val, merge(a.fields,b.fields)) +_subtractSurfaces(a::_T, b::_T) where _T<:ParaviewData = ParaviewData(a.x.val, a.y.val, a.z.val - b.z.val, merge(a.fields,b.fields)) + +""" + remove_NaN_Surface!(Z::Array,X::Array,Y::Array) + +Removes NaN's from a grid `Z` by taking the closest points as specified by `X` and `Y`. +""" +function remove_NaN_Surface!(Z,X,Y) + @assert size(Z) == size(X) == size(Y) + + # use nearest neighbour to interpolate data + id = findall(isnan.(Z) .== false) + id_NaN = findall(isnan.(Z)) + + coord = [X[id]'; Y[id]']; + kdtree = KDTree(coord; leafsize = 10); + + points = [X[id_NaN]'; Y[id_NaN]']; + idx,dist = nn(kdtree, points); + + Z[id_NaN] = Z[id[idx]] + + return nothing +end + + +""" + Topo = drape_on_topo(Topo::GeoData, Data::GeoData) + +This drapes fields of a data set `Data` on the topography `Topo` + +""" +function drape_on_topo(Topo::GeoData, Data::GeoData) + + Lon,Lat,_ = LonLatDepthGrid( Topo.lon.val[:,1,1], Topo.lat.val[1,:,1],Topo.depth.val[1,1,:]); + + # use nearest neighbour to interpolate data + coord = [vec(Data.lon.val)'; vec(Data.lat.val)']; + kdtree = KDTree(coord; leafsize = 10); + points = [vec(Lon)';vec(Lat)']; + idx,_ = nn(kdtree, points); + + idx_out = findall( (Lon .< minimum(Data.lon.val)) .| (Lon .> maximum(Data.lon.val)) .| + (Lat .< minimum(Data.lat.val)) .| (Lat .> maximum(Data.lat.val)) ) + + fields_new = Topo.fields; + field_names = keys(Data.fields); + + for i = 1:length(Data.fields) + + if typeof(Data.fields[i]) <: Tuple + + # vector or anything that contains more than 1 field + data_tuple = Data.fields[i] # we have a tuple (likely a vector field), so we have to loop + + data_array = zeros(typeof(data_tuple[1][1]),size(Topo.lon.val,1),size(Topo.lon.val,2),size(Topo.lon.val,3),length(Data.fields[i])); # create a 3D array that holds the 2D interpolated values + unit_array = zeros(size(data_array)); + + for j=1:length(data_tuple) + data_field = data_tuple[j]; + tmp = data_array[:,:,:,1]; + tmp = data_field[idx] + data_array[:,:,:,j] = tmp + end + + data_new = tuple([data_array[:,:,:,c] for c in 1:size(data_array,4)]...) # transform 4D matrix to tuple + + # remove points outside domain + for j=1:length(data_tuple) + tmp = data_new[j]; + tmp[idx_out] .= NaN + data_array[:,:,:,j] = tmp + end + data_new = tuple([data_array[:,:,:,c] for c in 1:size(data_array,4)]...) # transform 4D matrix to tuple + + else + + # scalar field + data_new = zeros(typeof(Data.fields[i][1]), size(Topo.lon.val,1),size(Topo.lon.val,2),size(Topo.lon.val,3)); + data_new = Data.fields[i][idx] # interpolate data field + + end + + # replace the one + new_field = NamedTuple{(field_names[i],)}((data_new,)) # Create a tuple with same name + fields_new = merge(fields_new, new_field); # replace the field in fields_new + + end + + Topo_new = GeoData(Topo.lon.val,Topo.lat.val,Topo.depth.val, fields_new) + + return Topo_new +end + + +""" + drape_on_topo(Topo::CartData, Data::CartData) + +Drapes Cartesian Data on topography +""" +function drape_on_topo(Topo::CartData, Data::CartData) + Topo_lonlat = GeoData(ustrip.(Topo.x.val),ustrip.(Topo.y.val), ustrip.(Topo.z.val), Topo.fields ) + Data_lonlat = GeoData(ustrip.(Data.x.val),ustrip.(Data.y.val), ustrip.(Data.z.val), Data.fields ) + + Topo_new_lonlat = drape_on_topo(Topo_lonlat, Data_lonlat) + + Topo_new = CartData(Topo_new_lonlat.lon.val, Topo_new_lonlat.lat.val, Topo_new_lonlat.depth.val, Topo_new_lonlat.fields) + + return Topo_new +end + + +""" + +This projects 3D points described by the vectors `lon_pt`,`lat_pt`, `depth_pt` to a GeoData surface `surf` + +""" +function project_points_to_surface(surf::GeoData, lon_pt::Vector, lat_pt::Vector, depth_pt::Vector) + +end diff --git a/src/utils.jl b/src/utils.jl index 5561c1a5..3c2b6543 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -4,10 +4,9 @@ export meshgrid, CrossSection, CrossSectionVolume, CrossSectionSurface, CrossSec export ParseColumns_CSV_File, AboveSurface, BelowSurface, VoteMap, CountMap export InterpolateDataOnSurface, InterpolateDataFields2D, InterpolateDataFields, InterpolateTopographyOnPlane export RotateTranslateScale -export DrapeOnTopo, LithostaticPressure! +export LithostaticPressure! export FlattenCrossSection export AddField, RemoveField -export SubtractSurfaces!, AddSurfaces!, remove_NaN_Surface! using NearestNeighbors @@ -1841,96 +1840,6 @@ function RotateTranslateScale(Data::Union{ParaviewData, CartData}; Rotate::Numbe end -""" - Topo = DrapeOnTopo(Topo::GeoData, Data::GeoData) - -This drapes fields of a data set `Data` on the topography `Topo` - - -""" -function DrapeOnTopo(Topo::GeoData, Data::GeoData) - - - Lon,Lat,Depth = LonLatDepthGrid( Topo.lon.val[:,1,1], Topo.lat.val[1,:,1],Topo.depth.val[1,1,:]); - - # use nearest neighbour to interpolate data - coord = [vec(Data.lon.val)'; vec(Data.lat.val)']; - kdtree = KDTree(coord; leafsize = 10); - points = [vec(Lon)';vec(Lat)']; - idx,dist = nn(kdtree, points); - - - idx_out = findall( (Lon .< minimum(Data.lon.val)) .| (Lon .> maximum(Data.lon.val)) .| - (Lat .< minimum(Data.lat.val)) .| (Lat .> maximum(Data.lat.val)) ) - - fields_new = Topo.fields; - field_names = keys(Data.fields); - - for i = 1:length(Data.fields) - - if typeof(Data.fields[i]) <: Tuple - - # vector or anything that contains more than 1 field - data_tuple = Data.fields[i] # we have a tuple (likely a vector field), so we have to loop - - data_array = zeros(typeof(data_tuple[1][1]),size(Topo.lon.val,1),size(Topo.lon.val,2),size(Topo.lon.val,3),length(Data.fields[i])); # create a 3D array that holds the 2D interpolated values - unit_array = zeros(size(data_array)); - - for j=1:length(data_tuple) - data_field = data_tuple[j]; - tmp = data_array[:,:,:,1]; - tmp = data_field[idx] - data_array[:,:,:,j] = tmp - end - - data_new = tuple([data_array[:,:,:,c] for c in 1:size(data_array,4)]...) # transform 4D matrix to tuple - - # remove points outside domain - for j=1:length(data_tuple) - tmp = data_new[j]; - tmp[idx_out] .= NaN - data_array[:,:,:,j] = tmp - end - data_new = tuple([data_array[:,:,:,c] for c in 1:size(data_array,4)]...) # transform 4D matrix to tuple - - else - - # scalar field - data_new = zeros(typeof(Data.fields[i][1]), size(Topo.lon.val,1),size(Topo.lon.val,2),size(Topo.lon.val,3)); - data_new = Data.fields[i][idx] # interpolate data field - - end - - # replace the one - new_field = NamedTuple{(field_names[i],)}((data_new,)) # Create a tuple with same name - fields_new = merge(fields_new, new_field); # replace the field in fields_new - - end - - - Topo_new = GeoData(Topo.lon.val,Topo.lat.val,Topo.depth.val, fields_new) - - return Topo_new - -end - - -""" - DrapeOnTopo(Topo::CartData, Data::CartData) - -Drapes Cartesian Data on topography -""" -function DrapeOnTopo(Topo::CartData, Data::CartData) - Topo_lonlat = GeoData(ustrip.(Topo.x.val),ustrip.(Topo.y.val), ustrip.(Topo.z.val), Topo.fields ) - Data_lonlat = GeoData(ustrip.(Data.x.val),ustrip.(Data.y.val), ustrip.(Data.z.val), Data.fields ) - - Topo_new_lonlat = DrapeOnTopo(Topo_lonlat, Data_lonlat) - - Topo_new = CartData(Topo_new_lonlat.lon.val, Topo_new_lonlat.lat.val, Topo_new_lonlat.depth.val, Topo_new_lonlat.fields) - - return Topo_new -end - """ LithostaticPressure!(Plithos::Array, Density::Array, dz::Number; g=9.81) @@ -1946,83 +1855,4 @@ function LithostaticPressure!(Plithos::Array{T,N}, Density::Array{T,N}, dz::Numb Plithos[:] = reverse!(cumsum(reverse!(Plithos),dims=N)) return nothing -end - - -""" - AddSurfaces!(Surface1::Union{GeoData,UTMData}, Surface2::Union{GeoData,UTMData}) - -Adds `Surface2` to `Surface1`. The addition happens on the `Surface1.depth`; the fields remain unchanged - -""" -function AddSurfaces!(Surface1::Union{GeoData,UTMData}, Surface2::Union{GeoData,UTMData}) - - Surface1.depth.val .= Surface1.depth.val + Surface2.depth.val - - return nothing -end - -""" - AddSurfaces!(Surface1::Union{CartData,ParaviewData}, Surface2::Union{CartData,ParaviewData}) - -Adds `Surface2` to `Surface1`. The addition happens on the `Surface1.z`; the fields remain unchanged - -""" -function AddSurfaces!(Surface1::Union{CartData,ParaviewData}, Surface2::Union{CartData,ParaviewData}) - - Surface1.z.val .= Surface1.z.val + Surface2.z.val - - return nothing -end - - -""" - SubtractSurfaces!(Surface1::Union{GeoData,UTMData}, Surface2::Union{GeoData,UTMData}) - -Subtracts `Surface2` to `Surface1`. The addition happens on the `Surface1.depth`; the fields remain unchanged - -""" -function SubtractSurfaces!(Surface1::Union{GeoData,UTMData}, Surface2::Union{GeoData,UTMData}) - - Surface1.depth.val .= Surface1.depth.val - Surface2.depth.val - - return nothing -end - - -""" - SubtractSurfaces!(Surface1::Union{CartData,ParaviewData}, Surface2::Union{CartData,ParaviewData}) - -Subtracts `Surface2` to `Surface1`. The addition happens on the `Surface1.z`; the fields remain unchanged - -""" -function SubtractSurfaces!(Surface1::Union{CartData,ParaviewData}, Surface2::Union{CartData,ParaviewData}) - - Surface1.z.val .= Surface1.z.val - Surface2.z.val - - return nothing -end - - - - -""" - remove_NaN_Surface!(Z,X,Y) - -Removes NaN's from a grid `Z` by taking the closest points as specified by `X` and `Y`. -""" -function remove_NaN_Surface!(Z,X,Y) - # use nearest neighbour to interpolate data - id = findall(isnan.(Z) .== false) - id_NaN = findall(isnan.(Z)) - - coord = [X[id]'; Y[id]']; - kdtree = KDTree(coord; leafsize = 10); - - points = [X[id_NaN]'; Y[id_NaN]']; - idx,dist = nn(kdtree, points); - - Z[id_NaN] = Z[id[idx]] - - return nothing -end +end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 90ec22bc..29ecf489 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -19,6 +19,9 @@ end @testset "Transformations" begin include("test_transformation.jl") end +@testset "Surfaces" begin + include("test_surfaces.jl") +end @testset "LaMEM" begin include("test_lamem.jl") diff --git a/tutorials/Tutorial_Jura.jl b/tutorials/Tutorial_Jura.jl index a6ba2623..67875625 100644 --- a/tutorials/Tutorial_Jura.jl +++ b/tutorials/Tutorial_Jura.jl @@ -27,7 +27,7 @@ upperright = [8.948117154811715, 47.781282316442606, 0.0] Geology = Screenshot_To_GeoData("SchoriM_Encl_01_Jura-map_A1.png", lowerleft, upperright, fieldname=:geology_colors) # name should have "colors" in it # You can "drape" this image on the topographic map with -TopoGeology = DrapeOnTopo(Topo, Geology) +TopoGeology = drape_on_topo(Topo, Geology) # ```julia # GeoData # size : (3721, 2641, 1)