Merge pull request #176 from JuliaGeodynamics/mt_profiledata

Added a screenshot field to the ProfileData structure

Author: mthielma

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
 """

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: