diff --git a/src/ProfileProcessing.jl b/src/ProfileProcessing.jl index 8930da9a..27476b0a 100644 --- a/src/ProfileProcessing.jl +++ b/src/ProfileProcessing.jl @@ -17,6 +17,7 @@ Structure that holds profile data (interpolated/projected on the profile) VolData :: GeophysicalModelGenerator.GeoData SurfData :: Union{Nothing, NamedTuple} PointData :: Union{Nothing, NamedTuple} + ScreenshotData ::Union{Nothing, NamedTuple} end Structure to store cross section data @@ -29,6 +30,7 @@ mutable struct ProfileData VolData::Union{Nothing, GeophysicalModelGenerator.GeoData} SurfData::Union{Nothing, NamedTuple} PointData::Union{Nothing, NamedTuple} + ScreenshotData::Union{Nothing, NamedTuple} function ProfileData(; kwargs...) # this constructor allows to define only certain fields and leave the others blank K = new(true, nothing, nothing, nothing, nothing, nothing, nothing) @@ -51,7 +53,6 @@ mutable struct ProfileData end end - function show(io::IO, g::ProfileData) if g.vertical println(io, "Vertical ProfileData") @@ -69,7 +70,9 @@ function show(io::IO, g::ProfileData) if !isnothing(g.PointData) println(io, " PointData: $(keys(g.PointData)) ") end - + if !isnothing(g.ScreenshotData) + println(io, " ScreenshotData: $(keys(g.ScreenshotData)) ") + end return nothing end @@ -282,6 +285,31 @@ function create_profile_volume!(Profile::ProfileData, VolData::AbstractGeneralGr return nothing end +### internal function to process screenshot data - contrary to the volume data, we here have to save lon/lat/depth pairs for every screenshot data set, so we create a NamedTuple of GeoData data sets +function create_profile_screenshot!(Profile::ProfileData, DataSet::NamedTuple) + num_datasets = length(DataSet) + + tmp = NamedTuple() # initialize empty one + DataSetName = keys(DataSet) # Names of the datasets + + for idata in 1:num_datasets + # load data set --> each data set is a single GeoData structure, so we'll only have to get the respective key to load the correct type + data_tmp = DataSet[idata] + if Profile.vertical + x_profile = flatten_cross_section(data_tmp, Start = Profile.start_lonlat) # compute the distance along the profile + data_tmp = addfield(data_tmp, "x_profile", x_profile) + + # add the data set as a NamedTuple + data_NT = NamedTuple{(DataSetName[idata],)}((data_tmp,)) + tmp = merge(tmp, data_NT) + else + # we do not have this implemented + #error("horizontal profiles not yet implemented") + end + end + Profile.SurfData = tmp # assign to profile data structure + return +end ### internal function to process surface data - contrary to the volume data, we here have to save lon/lat/depth pairs for every surface data set, so we create a NamedTuple of GeoData data sets function create_profile_surface!(Profile::ProfileData, DataSet::NamedTuple; DimsSurfCross = (100,)) @@ -362,21 +390,38 @@ end """ - extract_ProfileData!(Profile::ProfileData,VolData::GeoData, SurfData::NamedTuple, PointData::NamedTuple; DimsVolCross=(100,100),Depth_extent=nothing,DimsSurfCross=(100,),section_width=50) + extract_ProfileData!(Profile::ProfileData,VolData::GeoData, SurfData::NamedTuple, PointData::NamedTuple, ScreenshotData::NamedTuple; DimsVolCross=(100,100),Depth_extent=nothing,DimsSurfCross=(100,),section_width=50) Extracts data along a vertical or horizontal profile """ -function extract_ProfileData!(Profile::ProfileData, VolData::Union{Nothing, GeoData} = nothing, SurfData::NamedTuple = NamedTuple(), PointData::NamedTuple = NamedTuple(); DimsVolCross = (100, 100), Depth_extent = nothing, DimsSurfCross = (100,), section_width = 50km) +function extract_ProfileData!(Profile::ProfileData, VolData::Union{Nothing, GeoData} = nothing, SurfData::NamedTuple = NamedTuple(), PointData::NamedTuple = NamedTuple(),ScreenshotData::NamedTuple = NamedTuple(); DimsVolCross = (100, 100), Depth_extent = nothing, DimsSurfCross = (100,), section_width = 50km) if !isnothing(VolData) create_profile_volume!(Profile, VolData; DimsVolCross = DimsVolCross, Depth_extent = Depth_extent) end create_profile_surface!(Profile, SurfData, DimsSurfCross = DimsSurfCross) create_profile_point!(Profile, PointData, section_width = section_width) + if !isempty(ScreenshotData) + create_profile_screenshot!(Profile, ScreenshotData) + end + return nothing +end +""" + extract_ProfileData!(Profile::ProfileData,VolData::GeoData, SurfData::NamedTuple, PointData::NamedTuple; DimsVolCross=(100,100),Depth_extent=nothing,DimsSurfCross=(100,),section_width=50) + +Extracts data along a vertical or horizontal profile. +""" +function extract_ProfileData!(Profile::ProfileData, VolData::Union{Nothing, GeoData} = nothing, SurfData::NamedTuple = NamedTuple(), PointData::NamedTuple = NamedTuple(); DimsVolCross = (100, 100), Depth_extent = nothing, DimsSurfCross = (100,), section_width = 50km) + + # call the actual function to extract the profile data + extract_ProfileData!(Profile, VolData, SurfData, PointData,NamedTuple(); DimsVolCross = DimsVolCross, Depth_extent = Depth_extent, DimsSurfCross = DimsSurfCross, section_width = section_width) return nothing end + + + """ This reads the picked profiles from disk and returns a vector of ProfileData """ diff --git a/test/test_ProfileProcessing.jl b/test/test_ProfileProcessing.jl index 605d887f..0a27c173 100644 --- a/test/test_ProfileProcessing.jl +++ b/test/test_ProfileProcessing.jl @@ -22,7 +22,6 @@ data_Vol2 = GMG_Dataset("Plomerova2022", "Volume", "https://seafile.rlp.net/f/ab #data_Vol2 = GMG_Dataset("Zhao2016","Volume","https://seafile.rlp.net/f/e81a6d075f6746609973/?dl=1", true) # Now load these datasets into NamedTuples - SurfData = load_GMG(data_Surf) PointData = load_GMG(data_EQ) ScreenshotData = load_GMG(data_SS) @@ -59,6 +58,7 @@ prof1 = ProfileData(start_lonlat = (5, 45), end_lonlat = (15, 49)) prof2 = ProfileData(depth = -100) prof3 = ProfileData(start_lonlat = (5, 45), end_lonlat = (5, 49)) prof4 = ProfileData(depth = -20) +prof5 = ProfileData(start_lonlat = (12, 49), end_lonlat = (12, 45)) # test internal routines to intersect profile with volumetric data: GeophysicalModelGenerator.create_profile_volume!(prof1, VolData_combined1) @@ -80,6 +80,10 @@ GeophysicalModelGenerator.create_profile_point!(prof4, Data.Point, section_width @test length(prof1.PointData[1].lon) == 13 @test length(prof4.PointData[1].lon) == 445 +# test screenshot data +GeophysicalModelGenerator.create_profile_screenshot!(prof5, Data.Screenshot) +@test prof5.SurfData[1].fields.x_profile[1,1,1] == 0 + # Test the main profile extraction routines: extract_ProfileData!(prof1, VolData_combined1, Data.Surface, Data.Point) @@ -96,6 +100,7 @@ extract_ProfileData!(prof1, VolData_combined3, Data.Surface, Data.Point) extract_ProfileData!(prof2, VolData_combined3, Data.Surface, Data.Point) extract_ProfileData!(prof3, VolData_combined3, Data.Surface, Data.Point) extract_ProfileData!(prof4, VolData_combined3, Data.Surface, Data.Point) +extract_ProfileData!(prof5, VolData_combined3, Data.Surface, Data.Point, Data.Screenshot) # Test that it works if only EQ's are provided: