From 7f55db7f22141a08c850717e7f5700d47ad15ad3 Mon Sep 17 00:00:00 2001 From: Marcel Thielmann Date: Wed, 23 Sep 2026 16:31:32 +0200 Subject: [PATCH 01/25] Added a version of create_profile_volume! that allows to provide the datasets as named tuple. --- src/ProfileProcessing.jl | 64 +++++++++++++++++++++++++++++++++++++++- 1 file changed, 63 insertions(+), 1 deletion(-) diff --git a/src/ProfileProcessing.jl b/src/ProfileProcessing.jl index 6fe0aff95..294b10d5c 100644 --- a/src/ProfileProcessing.jl +++ b/src/ProfileProcessing.jl @@ -285,6 +285,69 @@ function create_profile_volume!(Profile::ProfileData, VolData::AbstractGeneralGr return nothing end +""" + create_profile_volume!(Profile::ProfileData, VolData::NamedTuple; DimsVolCross::NTuple=(100,100), Depth_extent=nothing) +Creates a cross-section through a volumetric 3D dataset `VolData` with the data supplied in `Profile`. `Depth_extent` can be the minimum & maximum depth for vertical profiles. This function allows to pass the data as NamedTuples instead of a GeoData object. +""" +function create_profile_volume!(Profile::ProfileData, VolData::NamedTuple; DimsVolCross::NTuple = (100, 100), Depth_extent = nothing) + + datasetnames = String.(keys(VolData)) # get the names of the datasets + + if Profile.vertical # take a vertical cross section + for ivol in eachindex(VolData) # loop over the different datasets and create a cross section through each of them + if ivol == 1 + cross_tmp = cross_section(VolData[1], dims = DimsVolCross, Start = Profile.start_lonlat, End = Profile.end_lonlat, Depth_extent = Depth_extent) # create the cross section + # flatten cross section and add this data to the structure + x_profile = flatten_cross_section(cross_tmp, Start = Profile.start_lonlat) # in the frist iteration, we create the x_profile field, which is the same for all datasets, so we only need to do this once + cross_tmp = addfield(cross_tmp, "x_profile", x_profile) + + # the issue is now that the fields do not contain any information about the originating dataset, so we add the name of the dataset to the field names + # we now do this by creating a new data structure named cross_add, which is then built up in the first iteration, and then merged with the next datasets in the following iterations + + cross_add = GeoData(cross_tmp.lon.val, cross_tmp.lat.val, cross_tmp.depth.val, (x_profile = cross_tmp.x_profile.val,)) # create a new GeoData structure with the x_profile field + + names_fields = String.(keys(cross_tmp.fields)) + for ifield in eachindex(names_fields) + name_new_field = datasetnames[ivol] * "_" * names_fields[ifield] # name of new field includes name of dataset + cross_add = addfield(cross_add, name_new_field, ustrip.(cross_tmp.fields[ifield])) # Note: we use ustrip here, and thereby remove the values, as the cross-section routine made problems + end + + else + cross_tmp = cross_section(VolData[ivol], dims = DimsVolCross, Start = Profile.start_lonlat, End = Profile.end_lonlat, Depth_extent = Depth_extent) # create the cross section + # add new fields to the cross_add structure, which already contains the data from the previous datasets + names_fields = String.(keys(cross_tmp.fields)) + for ifield in eachindex(names_fields) + name_new_field = datasetnames[ivol] * "_" * names_fields[ifield] # name of new field includes name of dataset + cross_add = addfield(cross_add, name_new_field, ustrip.(cross_tmp.fields[ifield])) # Note: we use ustrip here, and thereby remove the values, as the cross-section routine made problems + end + end + + end + + else # take a horizontal cross section - + if ivol == 1 + cross_tmp = cross_section(VolData, Depth_level = Profile.depth, Interpolate = true, dims = DimsVolCross) # create a horizontal cross section + cross_add = GeoData(cross_tmp.lon.val, cross_tmp.lat.val, cross_tmp.depth.val, (FlatCrossSection = cross_tmp.fields.FlatCrossSection,)) # create a basic cross section structure with the FlatCrossSection field + names_fields = String.(keys(cross_tmp.fields)) + for ifield in eachindex(names_fields) + name_new_field = datasetnames[ivol] * "_" * names_fields[ifield] # name of new field includes name of dataset + cross_add = addfield(cross_add, name_new_field, ustrip.(cross_tmp.fields[ifield])) # Note: we use ustrip here, and thereby remove the values, as the cross-section routine made problems + end + else + cross_tmp = cross_section(VolData[ivol], Depth_level = Profile.depth, Interpolate = true, dims = DimsVolCross) # create a horizontal cross section + names_fields = String.(keys(cross_tmp.fields)) + for ifield in eachindex(names_fields) + name_new_field = datasetnames[ivol] * "_" * names_fields[ifield] # name of new field includes name of dataset + cross_add = addfield(cross_add, name_new_field, ustrip.(cross_tmp.fields[ifield])) # Note: we use ustrip here, and thereby remove the values, as the cross-section routine made problems + end + end + end + Profile.VolData = cross_add # assign to Profile data structure + 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) @@ -413,7 +476,6 @@ function extract_ProfileData!(Profile::ProfileData, VolData::Union{Nothing, GeoD return nothing end - """ This reads the picked profiles from disk and returns a vector of ProfileData """ From c020e4e645a0bf0d1d0ffb9e62549e5207715de4 Mon Sep 17 00:00:00 2001 From: Marcel Thielmann Date: Wed, 23 Sep 2026 16:40:32 +0200 Subject: [PATCH 02/25] fixed typo --- src/ProfileProcessing.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ProfileProcessing.jl b/src/ProfileProcessing.jl index 294b10d5c..55970363f 100644 --- a/src/ProfileProcessing.jl +++ b/src/ProfileProcessing.jl @@ -298,7 +298,7 @@ function create_profile_volume!(Profile::ProfileData, VolData::NamedTuple; DimsV if ivol == 1 cross_tmp = cross_section(VolData[1], dims = DimsVolCross, Start = Profile.start_lonlat, End = Profile.end_lonlat, Depth_extent = Depth_extent) # create the cross section # flatten cross section and add this data to the structure - x_profile = flatten_cross_section(cross_tmp, Start = Profile.start_lonlat) # in the frist iteration, we create the x_profile field, which is the same for all datasets, so we only need to do this once + x_profile = flatten_cross_section(cross_tmp, Start = Profile.start_lonlat) # in the first iteration, we create the x_profile field, which is the same for all datasets, so we only need to do this once cross_tmp = addfield(cross_tmp, "x_profile", x_profile) # the issue is now that the fields do not contain any information about the originating dataset, so we add the name of the dataset to the field names From 5d9964cff9c4eaecab7c7da385770fdd8c46aff7 Mon Sep 17 00:00:00 2001 From: Marcel Thielmann Date: Wed, 23 Sep 2026 16:55:34 +0200 Subject: [PATCH 03/25] added variants of extract_ProfileData! that allow VolData to be passed as NamedTuple --- src/ProfileProcessing.jl | 34 +++++++++++++++++++++++++++++++++- 1 file changed, 33 insertions(+), 1 deletion(-) diff --git a/src/ProfileProcessing.jl b/src/ProfileProcessing.jl index 55970363f..853cf1bbf 100644 --- a/src/ProfileProcessing.jl +++ b/src/ProfileProcessing.jl @@ -264,7 +264,8 @@ end """ create_profile_volume!(Profile::ProfileData, VolData::AbstractGeneralGrid; DimsVolCross::NTuple=(100,100), Depth_extent=nothing) -Creates a cross-section through a volumetric 3D dataset `VolData` with the data supplied in `Profile`. `Depth_extent` can be the minimum & maximum depth for vertical profiles +Creates a cross-section through a volumetric 3D dataset `VolData` with the data supplied in `Profile`. `Depth_extent` can be the minimum & maximum depth for vertical profiles. + """ function create_profile_volume!(Profile::ProfileData, VolData::AbstractGeneralGrid; DimsVolCross::NTuple = (100, 100), Depth_extent = nothing) @@ -287,7 +288,9 @@ end """ create_profile_volume!(Profile::ProfileData, VolData::NamedTuple; DimsVolCross::NTuple=(100,100), Depth_extent=nothing) + Creates a cross-section through a volumetric 3D dataset `VolData` with the data supplied in `Profile`. `Depth_extent` can be the minimum & maximum depth for vertical profiles. This function allows to pass the data as NamedTuples instead of a GeoData object. + """ function create_profile_volume!(Profile::ProfileData, VolData::NamedTuple; DimsVolCross::NTuple = (100, 100), Depth_extent = nothing) @@ -452,6 +455,31 @@ function create_profile_point!(Profile::ProfileData, DataSet::NamedTuple; sectio end +""" + extract_ProfileData!(Profile::ProfileData,VolData::NamedTuple, SurfData::NamedTuple, PointData::NamedTuple; DimsVolCross=(100,100),Depth_extent=nothing,DimsSurfCross=(100,),section_width=50, ScreenshotData=nothing) + +Extracts data along a vertical or horizontal profile. Allows VolData to be passed as a NamedTuple. +""" +function extract_ProfileData!(Profile::ProfileData, VolData::NamedTuple = NamedTuple(), SurfData::NamedTuple = NamedTuple(), PointData::NamedTuple = NamedTuple(); DimsVolCross = (100, 100), Depth_extent = nothing, DimsSurfCross = (100,), section_width = 50km, ScreenshotData = nothing) + + return extract_ProfileData!(Profile, VolData, SurfData, PointData, ScreenshotData; DimsVolCross = DimsVolCross, Depth_extent = Depth_extent, DimsSurfCross = DimsSurfCross, section_width = section_width) +end + +# Internal method - called by the main method with ScreenshotData as positional argument, allows VolData as NamedTuple +function extract_ProfileData!(Profile::ProfileData, VolData::NamedTuple, SurfData::NamedTuple, PointData::NamedTuple, ScreenshotData::Union{Nothing, NamedTuple}; DimsVolCross = (100, 100), Depth_extent = nothing, DimsSurfCross = (100,), section_width = 50km) + + if !isempty(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 !isnothing(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, ScreenshotData=nothing) @@ -476,6 +504,10 @@ function extract_ProfileData!(Profile::ProfileData, VolData::Union{Nothing, GeoD return nothing end + + + + """ This reads the picked profiles from disk and returns a vector of ProfileData """ From dfddc05e10050299077be1f7a3b2b371e767f8b3 Mon Sep 17 00:00:00 2001 From: Marcel Thielmann Date: Wed, 23 Sep 2026 17:56:55 +0200 Subject: [PATCH 04/25] changes the input argumernt type from Union{Nothing, GeoData} to GeoData --- src/ProfileProcessing.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/ProfileProcessing.jl b/src/ProfileProcessing.jl index 853cf1bbf..14df38e8b 100644 --- a/src/ProfileProcessing.jl +++ b/src/ProfileProcessing.jl @@ -485,13 +485,13 @@ end 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, ScreenshotData = nothing) +function extract_ProfileData!(Profile::ProfileData, VolData::GeoData, SurfData::NamedTuple = NamedTuple(), PointData::NamedTuple = NamedTuple(); DimsVolCross = (100, 100), Depth_extent = nothing, DimsSurfCross = (100,), section_width = 50km, ScreenshotData = nothing) return extract_ProfileData!(Profile, VolData, SurfData, PointData, ScreenshotData; DimsVolCross = DimsVolCross, Depth_extent = Depth_extent, DimsSurfCross = DimsSurfCross, section_width = section_width) end # Internal method - called by the main method with ScreenshotData as positional argument -function extract_ProfileData!(Profile::ProfileData, VolData::Union{Nothing, GeoData}, SurfData::NamedTuple, PointData::NamedTuple, ScreenshotData::Union{Nothing, NamedTuple}; DimsVolCross = (100, 100), Depth_extent = nothing, DimsSurfCross = (100,), section_width = 50km) +function extract_ProfileData!(Profile::ProfileData, VolData::GeoData, SurfData::NamedTuple, PointData::NamedTuple, ScreenshotData::Union{Nothing, 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) From f318e52ed2ce28d9dc92dce81e9bb6f82fe8a44d Mon Sep 17 00:00:00 2001 From: Marcel Thielmann Date: Wed, 23 Sep 2026 22:39:42 +0200 Subject: [PATCH 05/25] added check for field additon --- src/ProfileProcessing.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/ProfileProcessing.jl b/src/ProfileProcessing.jl index 14df38e8b..1de37f269 100644 --- a/src/ProfileProcessing.jl +++ b/src/ProfileProcessing.jl @@ -311,7 +311,7 @@ function create_profile_volume!(Profile::ProfileData, VolData::NamedTuple; DimsV names_fields = String.(keys(cross_tmp.fields)) for ifield in eachindex(names_fields) - name_new_field = datasetnames[ivol] * "_" * names_fields[ifield] # name of new field includes name of dataset + name_new_field = datasetnames[1] * "_" * names_fields[ifield] # name of new field includes name of dataset cross_add = addfield(cross_add, name_new_field, ustrip.(cross_tmp.fields[ifield])) # Note: we use ustrip here, and thereby remove the values, as the cross-section routine made problems end @@ -321,6 +321,7 @@ function create_profile_volume!(Profile::ProfileData, VolData::NamedTuple; DimsV names_fields = String.(keys(cross_tmp.fields)) for ifield in eachindex(names_fields) name_new_field = datasetnames[ivol] * "_" * names_fields[ifield] # name of new field includes name of dataset + println("Adding field $name_new_field to cross_add") cross_add = addfield(cross_add, name_new_field, ustrip.(cross_tmp.fields[ifield])) # Note: we use ustrip here, and thereby remove the values, as the cross-section routine made problems end end From c044eff4214135075f95f50a64e6ca2767daa122 Mon Sep 17 00:00:00 2001 From: Marcel Thielmann Date: Wed, 23 Sep 2026 22:43:44 +0200 Subject: [PATCH 06/25] another test --- src/ProfileProcessing.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ProfileProcessing.jl b/src/ProfileProcessing.jl index 1de37f269..135a8f995 100644 --- a/src/ProfileProcessing.jl +++ b/src/ProfileProcessing.jl @@ -320,8 +320,8 @@ function create_profile_volume!(Profile::ProfileData, VolData::NamedTuple; DimsV # add new fields to the cross_add structure, which already contains the data from the previous datasets names_fields = String.(keys(cross_tmp.fields)) for ifield in eachindex(names_fields) + println("dataset: $(datasetnames[ivol]), field: $(names_fields[ifield])") name_new_field = datasetnames[ivol] * "_" * names_fields[ifield] # name of new field includes name of dataset - println("Adding field $name_new_field to cross_add") cross_add = addfield(cross_add, name_new_field, ustrip.(cross_tmp.fields[ifield])) # Note: we use ustrip here, and thereby remove the values, as the cross-section routine made problems end end From 7550f00a28fe3600fef490da6d1f469a94f0437d Mon Sep 17 00:00:00 2001 From: Marcel Thielmann Date: Wed, 23 Sep 2026 22:49:12 +0200 Subject: [PATCH 07/25] another try --- src/ProfileProcessing.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/ProfileProcessing.jl b/src/ProfileProcessing.jl index 135a8f995..522862a9e 100644 --- a/src/ProfileProcessing.jl +++ b/src/ProfileProcessing.jl @@ -319,6 +319,7 @@ function create_profile_volume!(Profile::ProfileData, VolData::NamedTuple; DimsV cross_tmp = cross_section(VolData[ivol], dims = DimsVolCross, Start = Profile.start_lonlat, End = Profile.end_lonlat, Depth_extent = Depth_extent) # create the cross section # add new fields to the cross_add structure, which already contains the data from the previous datasets names_fields = String.(keys(cross_tmp.fields)) + println("fields: $(names_fields)") for ifield in eachindex(names_fields) println("dataset: $(datasetnames[ivol]), field: $(names_fields[ifield])") name_new_field = datasetnames[ivol] * "_" * names_fields[ifield] # name of new field includes name of dataset From a2e795792998bce9fc0679f913d11b1f2c023dec Mon Sep 17 00:00:00 2001 From: Marcel Thielmann Date: Wed, 23 Sep 2026 23:02:56 +0200 Subject: [PATCH 08/25] potential fix --- src/ProfileProcessing.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ProfileProcessing.jl b/src/ProfileProcessing.jl index 522862a9e..5528d19dd 100644 --- a/src/ProfileProcessing.jl +++ b/src/ProfileProcessing.jl @@ -297,7 +297,7 @@ function create_profile_volume!(Profile::ProfileData, VolData::NamedTuple; DimsV datasetnames = String.(keys(VolData)) # get the names of the datasets if Profile.vertical # take a vertical cross section - for ivol in eachindex(VolData) # loop over the different datasets and create a cross section through each of them + for ivol in eachindex(datasetnames) # loop over the different datasets and create a cross section through each of them if ivol == 1 cross_tmp = cross_section(VolData[1], dims = DimsVolCross, Start = Profile.start_lonlat, End = Profile.end_lonlat, Depth_extent = Depth_extent) # create the cross section # flatten cross section and add this data to the structure From 4163cd0c1b68f57d19a23ff9c6500bba223c0f78 Mon Sep 17 00:00:00 2001 From: Marcel Thielmann Date: Wed, 23 Sep 2026 23:05:47 +0200 Subject: [PATCH 09/25] next bugfix --- src/ProfileProcessing.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ProfileProcessing.jl b/src/ProfileProcessing.jl index 5528d19dd..ce33ef651 100644 --- a/src/ProfileProcessing.jl +++ b/src/ProfileProcessing.jl @@ -307,7 +307,7 @@ function create_profile_volume!(Profile::ProfileData, VolData::NamedTuple; DimsV # the issue is now that the fields do not contain any information about the originating dataset, so we add the name of the dataset to the field names # we now do this by creating a new data structure named cross_add, which is then built up in the first iteration, and then merged with the next datasets in the following iterations - cross_add = GeoData(cross_tmp.lon.val, cross_tmp.lat.val, cross_tmp.depth.val, (x_profile = cross_tmp.x_profile.val,)) # create a new GeoData structure with the x_profile field + cross_add = GeoData(cross_tmp.lon.val, cross_tmp.lat.val, cross_tmp.depth.val, (x_profile = cross_tmp.fields.x_profile.val,)) # create a new GeoData structure with the x_profile field names_fields = String.(keys(cross_tmp.fields)) for ifield in eachindex(names_fields) From b51a5be24f6e5f83766a24c95cce3f2805693be2 Mon Sep 17 00:00:00 2001 From: Marcel Thielmann Date: Wed, 23 Sep 2026 23:09:45 +0200 Subject: [PATCH 10/25] another bugfix --- src/ProfileProcessing.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ProfileProcessing.jl b/src/ProfileProcessing.jl index ce33ef651..1b83e925e 100644 --- a/src/ProfileProcessing.jl +++ b/src/ProfileProcessing.jl @@ -307,7 +307,7 @@ function create_profile_volume!(Profile::ProfileData, VolData::NamedTuple; DimsV # the issue is now that the fields do not contain any information about the originating dataset, so we add the name of the dataset to the field names # we now do this by creating a new data structure named cross_add, which is then built up in the first iteration, and then merged with the next datasets in the following iterations - cross_add = GeoData(cross_tmp.lon.val, cross_tmp.lat.val, cross_tmp.depth.val, (x_profile = cross_tmp.fields.x_profile.val,)) # create a new GeoData structure with the x_profile field + cross_add = GeoData(cross_tmp.lon.val, cross_tmp.lat.val, cross_tmp.depth.val, (x_profile = cross_tmp.fields.x_profile,)) # create a new GeoData structure with the x_profile field names_fields = String.(keys(cross_tmp.fields)) for ifield in eachindex(names_fields) From 3803d54dbdd099642f7b2f7c247fc104e6f47b8d Mon Sep 17 00:00:00 2001 From: Marcel Thielmann Date: Wed, 23 Sep 2026 23:12:00 +0200 Subject: [PATCH 11/25] cleanup --- src/ProfileProcessing.jl | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/ProfileProcessing.jl b/src/ProfileProcessing.jl index 1b83e925e..bff64b13f 100644 --- a/src/ProfileProcessing.jl +++ b/src/ProfileProcessing.jl @@ -319,9 +319,7 @@ function create_profile_volume!(Profile::ProfileData, VolData::NamedTuple; DimsV cross_tmp = cross_section(VolData[ivol], dims = DimsVolCross, Start = Profile.start_lonlat, End = Profile.end_lonlat, Depth_extent = Depth_extent) # create the cross section # add new fields to the cross_add structure, which already contains the data from the previous datasets names_fields = String.(keys(cross_tmp.fields)) - println("fields: $(names_fields)") - for ifield in eachindex(names_fields) - println("dataset: $(datasetnames[ivol]), field: $(names_fields[ifield])") + for ifield in eachindex(names_fields) name_new_field = datasetnames[ivol] * "_" * names_fields[ifield] # name of new field includes name of dataset cross_add = addfield(cross_add, name_new_field, ustrip.(cross_tmp.fields[ifield])) # Note: we use ustrip here, and thereby remove the values, as the cross-section routine made problems end From e56d175ec32987caf73e608e4f71a2576129193f Mon Sep 17 00:00:00 2001 From: Marcel Thielmann Date: Wed, 23 Sep 2026 23:52:53 +0200 Subject: [PATCH 12/25] added a topo data field and updated the corresponding functions --- src/ProfileProcessing.jl | 60 ++++++++++++++++++++++++++++++++++------ 1 file changed, 52 insertions(+), 8 deletions(-) diff --git a/src/ProfileProcessing.jl b/src/ProfileProcessing.jl index bff64b13f..854213b32 100644 --- a/src/ProfileProcessing.jl +++ b/src/ProfileProcessing.jl @@ -16,6 +16,7 @@ Structure that holds profile data (interpolated/projected on the profile) depth :: Union{Nothing,Float64} VolData :: GeophysicalModelGenerator.GeoData SurfData :: Union{Nothing, NamedTuple} + TopoData :: Union{Nothing, NamedTuple} PointData :: Union{Nothing, NamedTuple} ScreenshotData :: Union{Nothing, NamedTuple} end @@ -29,6 +30,7 @@ mutable struct ProfileData depth::Union{Nothing, Float64} VolData::Union{Nothing, GeophysicalModelGenerator.GeoData} SurfData::Union{Nothing, NamedTuple} + TopoData::Union{Nothing, NamedTuple} PointData::Union{Nothing, NamedTuple} ScreenshotData::Union{Nothing, NamedTuple} @@ -67,6 +69,9 @@ function show(io::IO, g::ProfileData) if !isnothing(g.SurfData) println(io, " SurfData : $(keys(g.SurfData)) ") end + if !isnothing(g.TopoData) + println(io, " TopoData : $(keys(g.TopoData)) ") + end if !isnothing(g.PointData) println(io, " PointData: $(keys(g.PointData)) ") end @@ -349,8 +354,6 @@ function create_profile_volume!(Profile::ProfileData, VolData::NamedTuple; DimsV 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) @@ -411,6 +414,41 @@ function create_profile_surface!(Profile::ProfileData, DataSet::NamedTuple; Dims end +### internal function to process topography 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 +### this is essentially the same as create_profile_surface!, but we save the data in Profile.TopoData instead of Profile.SurfData +function create_profile_topo!(Profile::ProfileData, DataSet::NamedTuple; DimsSurfCross = (100,)) + 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 + # take a vertical cross section + data = cross_section_surface(data_tmp, dims = DimsSurfCross, Start = Profile.start_lonlat, End = Profile.end_lonlat) # create the cross section + + # flatten cross section and add this data to the structure + x_profile = flatten_cross_section(data, Start = Profile.start_lonlat) + data = addfield(data, "x_profile", x_profile) + + # add the data set as a NamedTuple + data_NT = NamedTuple{(DataSetName[idata],)}((data,)) + tmp = merge(tmp, data_NT) + + else + # we do not have this implemented + #error("horizontal profiles not yet implemented") + end + end + + Profile.TopoData = tmp # assign to profile data structure + return +end + + ### function to process point data - contrary to the volume data, we here have to save lon/lat/depth pairs for every point data set function create_profile_point!(Profile::ProfileData, DataSet::NamedTuple; section_width = 50km) num_datasets = length(DataSet) @@ -460,19 +498,22 @@ end Extracts data along a vertical or horizontal profile. Allows VolData to be passed as a NamedTuple. """ -function extract_ProfileData!(Profile::ProfileData, VolData::NamedTuple = NamedTuple(), SurfData::NamedTuple = NamedTuple(), PointData::NamedTuple = NamedTuple(); DimsVolCross = (100, 100), Depth_extent = nothing, DimsSurfCross = (100,), section_width = 50km, ScreenshotData = nothing) +function extract_ProfileData!(Profile::ProfileData, VolData::NamedTuple = NamedTuple(), SurfData::NamedTuple = NamedTuple(), PointData::NamedTuple = NamedTuple(); DimsVolCross = (100, 100), Depth_extent = nothing, DimsSurfCross = (100,), section_width = 50km, ScreenshotData = nothing, TopoData::NamedTuple = NamedTuple()) - return extract_ProfileData!(Profile, VolData, SurfData, PointData, ScreenshotData; DimsVolCross = DimsVolCross, Depth_extent = Depth_extent, DimsSurfCross = DimsSurfCross, section_width = section_width) + return extract_ProfileData!(Profile, VolData, SurfData, PointData, TopoData, ScreenshotData; DimsVolCross = DimsVolCross, Depth_extent = Depth_extent, DimsSurfCross = DimsSurfCross, section_width = section_width) end # Internal method - called by the main method with ScreenshotData as positional argument, allows VolData as NamedTuple -function extract_ProfileData!(Profile::ProfileData, VolData::NamedTuple, SurfData::NamedTuple, PointData::NamedTuple, ScreenshotData::Union{Nothing, NamedTuple}; DimsVolCross = (100, 100), Depth_extent = nothing, DimsSurfCross = (100,), section_width = 50km) +function extract_ProfileData!(Profile::ProfileData, VolData::NamedTuple, SurfData::NamedTuple, PointData::NamedTuple, TopoData:: Union{Nothing, NamedTuple},ScreenshotData::Union{Nothing, NamedTuple}; DimsVolCross = (100, 100), Depth_extent = nothing, DimsSurfCross = (100,), section_width = 50km) if !isempty(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 !isnothing(TopoData) + create_profile_topo!(Profile, TopoData, DimsSurfCross = DimsSurfCross*5) # we use a larger number of points for the topography, as it is often more detailed than the surface data + end if !isnothing(ScreenshotData) create_profile_screenshot!(Profile, ScreenshotData) end @@ -485,19 +526,22 @@ end Extracts data along a vertical or horizontal profile """ -function extract_ProfileData!(Profile::ProfileData, VolData::GeoData, SurfData::NamedTuple = NamedTuple(), PointData::NamedTuple = NamedTuple(); DimsVolCross = (100, 100), Depth_extent = nothing, DimsSurfCross = (100,), section_width = 50km, ScreenshotData = nothing) +function extract_ProfileData!(Profile::ProfileData, VolData::GeoData, SurfData::NamedTuple = NamedTuple(), PointData::NamedTuple = NamedTuple(); DimsVolCross = (100, 100), Depth_extent = nothing, DimsSurfCross = (100,), section_width = 50km, ScreenshotData = nothing, TopoData::NamedTuple = NamedTuple()) - return extract_ProfileData!(Profile, VolData, SurfData, PointData, ScreenshotData; DimsVolCross = DimsVolCross, Depth_extent = Depth_extent, DimsSurfCross = DimsSurfCross, section_width = section_width) + return extract_ProfileData!(Profile, VolData, SurfData, PointData, TopoData, ScreenshotData; DimsVolCross = DimsVolCross, Depth_extent = Depth_extent, DimsSurfCross = DimsSurfCross, section_width = section_width) end # Internal method - called by the main method with ScreenshotData as positional argument -function extract_ProfileData!(Profile::ProfileData, VolData::GeoData, SurfData::NamedTuple, PointData::NamedTuple, ScreenshotData::Union{Nothing, NamedTuple}; DimsVolCross = (100, 100), Depth_extent = nothing, DimsSurfCross = (100,), section_width = 50km) +function extract_ProfileData!(Profile::ProfileData, VolData::GeoData, SurfData::NamedTuple, PointData::NamedTuple, TopoData::Union{Nothing, NamedTuple}, ScreenshotData::Union{Nothing, 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 !isnothing(TopoData) + create_profile_topo!(Profile, TopoData, DimsSurfCross = DimsSurfCross*5) # we use a larger number of points for the topography, as it is often more detailed than the surface data + end if !isnothing(ScreenshotData) create_profile_screenshot!(Profile, ScreenshotData) end From 041023c4d8ad33817fa1cfa906ba36bbaae1597e Mon Sep 17 00:00:00 2001 From: Marcel Thielmann Date: Wed, 23 Sep 2026 23:55:44 +0200 Subject: [PATCH 13/25] bugfix --- src/ProfileProcessing.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/ProfileProcessing.jl b/src/ProfileProcessing.jl index 854213b32..ec922e844 100644 --- a/src/ProfileProcessing.jl +++ b/src/ProfileProcessing.jl @@ -512,7 +512,7 @@ function extract_ProfileData!(Profile::ProfileData, VolData::NamedTuple, SurfDat create_profile_surface!(Profile, SurfData, DimsSurfCross = DimsSurfCross) create_profile_point!(Profile, PointData, section_width = section_width) if !isnothing(TopoData) - create_profile_topo!(Profile, TopoData, DimsSurfCross = DimsSurfCross*5) # we use a larger number of points for the topography, as it is often more detailed than the surface data + create_profile_topo!(Profile, TopoData, DimsSurfCross = 5 .* DimsSurfCross) # we use a larger number of points for the topography, as it is often more detailed than the surface data end if !isnothing(ScreenshotData) create_profile_screenshot!(Profile, ScreenshotData) @@ -540,7 +540,7 @@ function extract_ProfileData!(Profile::ProfileData, VolData::GeoData, SurfData:: create_profile_surface!(Profile, SurfData, DimsSurfCross = DimsSurfCross) create_profile_point!(Profile, PointData, section_width = section_width) if !isnothing(TopoData) - create_profile_topo!(Profile, TopoData, DimsSurfCross = DimsSurfCross*5) # we use a larger number of points for the topography, as it is often more detailed than the surface data + create_profile_topo!(Profile, TopoData, DimsSurfCross = 5 .* DimsSurfCross) # we use a larger number of points for the topography, as it is often more detailed than the surface data end if !isnothing(ScreenshotData) create_profile_screenshot!(Profile, ScreenshotData) From 8a43f3dab9485963cd5108c9607353bb28dc8285 Mon Sep 17 00:00:00 2001 From: Marcel Thielmann Date: Thu, 24 Sep 2026 00:12:08 +0200 Subject: [PATCH 14/25] some fixes to the function documentation --- src/ProfileProcessing.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/ProfileProcessing.jl b/src/ProfileProcessing.jl index ec922e844..37f704852 100644 --- a/src/ProfileProcessing.jl +++ b/src/ProfileProcessing.jl @@ -494,7 +494,7 @@ end """ - extract_ProfileData!(Profile::ProfileData,VolData::NamedTuple, SurfData::NamedTuple, PointData::NamedTuple; DimsVolCross=(100,100),Depth_extent=nothing,DimsSurfCross=(100,),section_width=50, ScreenshotData=nothing) + extract_ProfileData!(Profile::ProfileData,VolData::NamedTuple, SurfData::NamedTuple, PointData::NamedTuple; DimsVolCross=(100,100),Depth_extent=nothing,DimsSurfCross=(100,),section_width=50, ScreenshotData=nothing, TopoData=NamedTuple()) Extracts data along a vertical or horizontal profile. Allows VolData to be passed as a NamedTuple. """ @@ -522,7 +522,7 @@ end """ - extract_ProfileData!(Profile::ProfileData,VolData::GeoData, SurfData::NamedTuple, PointData::NamedTuple; DimsVolCross=(100,100),Depth_extent=nothing,DimsSurfCross=(100,),section_width=50, ScreenshotData=nothing) + extract_ProfileData!(Profile::ProfileData,VolData::GeoData, SurfData::NamedTuple, PointData::NamedTuple; DimsVolCross=(100,100),Depth_extent=nothing,DimsSurfCross=(100,),section_width=50, ScreenshotData=nothing, TopoData=NamedTuple()) Extracts data along a vertical or horizontal profile """ From 0581e27c28e060e2ae2778e7c8695d50b101a0f2 Mon Sep 17 00:00:00 2001 From: Marcel Thielmann Date: Thu, 24 Sep 2026 07:31:25 +0200 Subject: [PATCH 15/25] started to add tests --- src/ProfileProcessing.jl | 36 ++++++++++++++++++---------------- test/gmt.history | 4 ++-- test/test_ProfileProcessing.jl | 10 ++++++++++ 3 files changed, 31 insertions(+), 19 deletions(-) diff --git a/src/ProfileProcessing.jl b/src/ProfileProcessing.jl index 37f704852..2ba7824e7 100644 --- a/src/ProfileProcessing.jl +++ b/src/ProfileProcessing.jl @@ -317,7 +317,7 @@ function create_profile_volume!(Profile::ProfileData, VolData::NamedTuple; DimsV names_fields = String.(keys(cross_tmp.fields)) for ifield in eachindex(names_fields) name_new_field = datasetnames[1] * "_" * names_fields[ifield] # name of new field includes name of dataset - cross_add = addfield(cross_add, name_new_field, ustrip.(cross_tmp.fields[ifield])) # Note: we use ustrip here, and thereby remove the values, as the cross-section routine made problems + cross_add = addfield(cross_add, name_new_field, ustrip.(cross_tmp.fields[ifield])) # Note: we use ustrip here, and thereby remove the units, as the cross-section routine made problems end else @@ -326,27 +326,29 @@ function create_profile_volume!(Profile::ProfileData, VolData::NamedTuple; DimsV names_fields = String.(keys(cross_tmp.fields)) for ifield in eachindex(names_fields) name_new_field = datasetnames[ivol] * "_" * names_fields[ifield] # name of new field includes name of dataset - cross_add = addfield(cross_add, name_new_field, ustrip.(cross_tmp.fields[ifield])) # Note: we use ustrip here, and thereby remove the values, as the cross-section routine made problems + cross_add = addfield(cross_add, name_new_field, ustrip.(cross_tmp.fields[ifield])) # Note: we use ustrip here, and thereby remove the units, as the cross-section routine made problems end end end - else # take a horizontal cross section - - if ivol == 1 - cross_tmp = cross_section(VolData, Depth_level = Profile.depth, Interpolate = true, dims = DimsVolCross) # create a horizontal cross section - cross_add = GeoData(cross_tmp.lon.val, cross_tmp.lat.val, cross_tmp.depth.val, (FlatCrossSection = cross_tmp.fields.FlatCrossSection,)) # create a basic cross section structure with the FlatCrossSection field - names_fields = String.(keys(cross_tmp.fields)) - for ifield in eachindex(names_fields) - name_new_field = datasetnames[ivol] * "_" * names_fields[ifield] # name of new field includes name of dataset - cross_add = addfield(cross_add, name_new_field, ustrip.(cross_tmp.fields[ifield])) # Note: we use ustrip here, and thereby remove the values, as the cross-section routine made problems - end - else - cross_tmp = cross_section(VolData[ivol], Depth_level = Profile.depth, Interpolate = true, dims = DimsVolCross) # create a horizontal cross section - names_fields = String.(keys(cross_tmp.fields)) - for ifield in eachindex(names_fields) - name_new_field = datasetnames[ivol] * "_" * names_fields[ifield] # name of new field includes name of dataset - cross_add = addfield(cross_add, name_new_field, ustrip.(cross_tmp.fields[ifield])) # Note: we use ustrip here, and thereby remove the values, as the cross-section routine made problems + else # take a horizontal cross section + for ivol in eachindex(datasetnames) # loop over the different datasets and create a cross section through each of them - + if ivol == 1 + cross_tmp = cross_section(VolData, Depth_level = Profile.depth, Interpolate = true, dims = DimsVolCross) # create a horizontal cross section + cross_add = GeoData(cross_tmp.lon.val, cross_tmp.lat.val, cross_tmp.depth.val, (FlatCrossSection = cross_tmp.fields.FlatCrossSection,)) # create a basic cross section structure with the FlatCrossSection field + names_fields = String.(keys(cross_tmp.fields)) + for ifield in eachindex(names_fields) + name_new_field = datasetnames[ivol] * "_" * names_fields[ifield] # name of new field includes name of dataset + cross_add = addfield(cross_add, name_new_field, ustrip.(cross_tmp.fields[ifield])) # Note: we use ustrip here, and thereby remove the units, as the cross-section routine made problems + end + else + cross_tmp = cross_section(VolData[ivol], Depth_level = Profile.depth, Interpolate = true, dims = DimsVolCross) # create a horizontal cross section + names_fields = String.(keys(cross_tmp.fields)) + for ifield in eachindex(names_fields) + name_new_field = datasetnames[ivol] * "_" * names_fields[ifield] # name of new field includes name of dataset + cross_add = addfield(cross_add, name_new_field, ustrip.(cross_tmp.fields[ifield])) # Note: we use ustrip here, and thereby remove the units, as the cross-section routine made problems + end end end end diff --git a/test/gmt.history b/test/gmt.history index ebcde8711..9ffd9e24a 100644 --- a/test/gmt.history +++ b/test/gmt.history @@ -1,4 +1,4 @@ # GMT 6 Session common arguments shelf -BEGIN GMT 6.6.0 -R 6.5/7.3/50.2/50.6 +BEGIN GMT 6.8.0 +R 341.8/342.5/28.4/29.0 END diff --git a/test/test_ProfileProcessing.jl b/test/test_ProfileProcessing.jl index 0a27c1738..040198d8a 100644 --- a/test/test_ProfileProcessing.jl +++ b/test/test_ProfileProcessing.jl @@ -70,6 +70,16 @@ GeophysicalModelGenerator.create_profile_volume!(prof2, VolData_combined1) GeophysicalModelGenerator.create_profile_volume!(prof1, VolData_combined1, Depth_extent = (-300, -100)) @test extrema(prof1.VolData.depth.val) == (-300.0, -100.0) +# test routines with volumetric data, but with NamedTuple instead of GeoData +GeophysicalModelGenerator.create_profile_volume!(prof1, Data.Volume) +@test prof1.VolData.fields.Hua2017_Vp[30, 40] ≈ 9.141520976523731 + +GeophysicalModelGenerator.create_profile_volume!(prof2, Data.Volume) +@test prof2.VolData.fields.Hua2017_Vp[30, 40] ≈ 8.177263544536272 + +GeophysicalModelGenerator.create_profile_volume!(prof1, Data.Volume, Depth_extent = (-300, -100)) +@test extrema(prof1.VolData.depth.val) == (-300.0, -100.0) + # Intersect surface data: GeophysicalModelGenerator.create_profile_surface!(prof1, Data.Surface) @test prof1.SurfData[1].fields.MohoDepth[80] ≈ -37.58791461075397km From b34bf08bb4a60bc1a11a879315bc6643faa1f27c Mon Sep 17 00:00:00 2001 From: Marcel Thielmann Date: Thu, 24 Sep 2026 08:16:10 +0200 Subject: [PATCH 16/25] fixed create_profile_volume! --- src/ProfileProcessing.jl | 78 ++++++++++++++++++++-------------------- 1 file changed, 39 insertions(+), 39 deletions(-) diff --git a/src/ProfileProcessing.jl b/src/ProfileProcessing.jl index 2ba7824e7..94d1f0500 100644 --- a/src/ProfileProcessing.jl +++ b/src/ProfileProcessing.jl @@ -302,53 +302,53 @@ function create_profile_volume!(Profile::ProfileData, VolData::NamedTuple; DimsV datasetnames = String.(keys(VolData)) # get the names of the datasets if Profile.vertical # take a vertical cross section - for ivol in eachindex(datasetnames) # loop over the different datasets and create a cross section through each of them - if ivol == 1 - cross_tmp = cross_section(VolData[1], dims = DimsVolCross, Start = Profile.start_lonlat, End = Profile.end_lonlat, Depth_extent = Depth_extent) # create the cross section - # flatten cross section and add this data to the structure - x_profile = flatten_cross_section(cross_tmp, Start = Profile.start_lonlat) # in the first iteration, we create the x_profile field, which is the same for all datasets, so we only need to do this once - cross_tmp = addfield(cross_tmp, "x_profile", x_profile) - # the issue is now that the fields do not contain any information about the originating dataset, so we add the name of the dataset to the field names - # we now do this by creating a new data structure named cross_add, which is then built up in the first iteration, and then merged with the next datasets in the following iterations + # process the first dataset, and create the cross_add structure + #----------------------------------------------------------------------- + cross_tmp = cross_section(VolData[1], dims = DimsVolCross, Start = Profile.start_lonlat, End = Profile.end_lonlat, Depth_extent = Depth_extent) # create the cross section + # flatten cross section and add this data to the structure + x_profile = flatten_cross_section(cross_tmp, Start = Profile.start_lonlat) # in the first iteration, we create the x_profile field, which is the same for all datasets, so we only need to do this once + cross_tmp = addfield(cross_tmp, "x_profile", x_profile) + + # the issue is now that the fields do not contain any information about the originating dataset, so we add the name of the dataset to the field names + # we now do this by creating a new data structure named cross_add, which is then built up in the first iteration, and then merged with the next datasets in the following iterations - cross_add = GeoData(cross_tmp.lon.val, cross_tmp.lat.val, cross_tmp.depth.val, (x_profile = cross_tmp.fields.x_profile,)) # create a new GeoData structure with the x_profile field + cross_add = GeoData(cross_tmp.lon.val, cross_tmp.lat.val, cross_tmp.depth.val, (x_profile = cross_tmp.fields.x_profile,)) # create a new GeoData structure with the x_profile field - names_fields = String.(keys(cross_tmp.fields)) - for ifield in eachindex(names_fields) - name_new_field = datasetnames[1] * "_" * names_fields[ifield] # name of new field includes name of dataset - cross_add = addfield(cross_add, name_new_field, ustrip.(cross_tmp.fields[ifield])) # Note: we use ustrip here, and thereby remove the units, as the cross-section routine made problems - end + names_fields = String.(keys(cross_tmp.fields)) + for ifield in eachindex(names_fields) + name_new_field = datasetnames[1] * "_" * names_fields[ifield] # name of new field includes name of dataset + cross_add = addfield(cross_add, name_new_field, ustrip.(cross_tmp.fields[ifield])) # Note: we use ustrip here, and thereby remove the units, as the cross-section routine made problems + end - else - cross_tmp = cross_section(VolData[ivol], dims = DimsVolCross, Start = Profile.start_lonlat, End = Profile.end_lonlat, Depth_extent = Depth_extent) # create the cross section - # add new fields to the cross_add structure, which already contains the data from the previous datasets - names_fields = String.(keys(cross_tmp.fields)) - for ifield in eachindex(names_fields) - name_new_field = datasetnames[ivol] * "_" * names_fields[ifield] # name of new field includes name of dataset - cross_add = addfield(cross_add, name_new_field, ustrip.(cross_tmp.fields[ifield])) # Note: we use ustrip here, and thereby remove the units, as the cross-section routine made problems - end + # loop over the remaining datasets and add the result to the cross_add structure + for ivol in eachindex(datasetnames)[2:end] + cross_tmp = cross_section(VolData[ivol], dims = DimsVolCross, Start = Profile.start_lonlat, End = Profile.end_lonlat, Depth_extent = Depth_extent) # create the cross section + # add new fields to the cross_add structure, which already contains the data from the previous datasets + names_fields = String.(keys(cross_tmp.fields)) + for ifield in eachindex(names_fields) + name_new_field = datasetnames[ivol] * "_" * names_fields[ifield] # name of new field includes name of dataset + cross_add = addfield(cross_add, name_new_field, ustrip.(cross_tmp.fields[ifield])) # Note: we use ustrip here, and thereby remove the units, as the cross-section routine made problems end - end else # take a horizontal cross section - for ivol in eachindex(datasetnames) # loop over the different datasets and create a cross section through each of them - - if ivol == 1 - cross_tmp = cross_section(VolData, Depth_level = Profile.depth, Interpolate = true, dims = DimsVolCross) # create a horizontal cross section - cross_add = GeoData(cross_tmp.lon.val, cross_tmp.lat.val, cross_tmp.depth.val, (FlatCrossSection = cross_tmp.fields.FlatCrossSection,)) # create a basic cross section structure with the FlatCrossSection field - names_fields = String.(keys(cross_tmp.fields)) - for ifield in eachindex(names_fields) - name_new_field = datasetnames[ivol] * "_" * names_fields[ifield] # name of new field includes name of dataset - cross_add = addfield(cross_add, name_new_field, ustrip.(cross_tmp.fields[ifield])) # Note: we use ustrip here, and thereby remove the units, as the cross-section routine made problems - end - else - cross_tmp = cross_section(VolData[ivol], Depth_level = Profile.depth, Interpolate = true, dims = DimsVolCross) # create a horizontal cross section - names_fields = String.(keys(cross_tmp.fields)) - for ifield in eachindex(names_fields) - name_new_field = datasetnames[ivol] * "_" * names_fields[ifield] # name of new field includes name of dataset - cross_add = addfield(cross_add, name_new_field, ustrip.(cross_tmp.fields[ifield])) # Note: we use ustrip here, and thereby remove the units, as the cross-section routine made problems - end + + # process the first dataset, and create the cross_add structure + cross_tmp = cross_section(VolData[1], Depth_level = Profile.depth, Interpolate = true, dims = DimsVolCross) # create a horizontal cross section + cross_add = GeoData(cross_tmp.lon.val, cross_tmp.lat.val, cross_tmp.depth.val, (FlatCrossSection = cross_tmp.fields.FlatCrossSection,)) # create a basic cross section structure with the FlatCrossSection field + names_fields = String.(keys(cross_tmp.fields)) + for ifield in eachindex(names_fields) + name_new_field = datasetnames[1] * "_" * names_fields[ifield] # name of new field includes name of dataset + cross_add = addfield(cross_add, name_new_field, ustrip.(cross_tmp.fields[ifield])) # Note: we use ustrip here, and thereby remove the units, as the cross-section routine made problems + end + + for ivol in eachindex(datasetnames)[2:end] # loop over the remaining datasets and create a cross section through each of them + cross_tmp = cross_section(VolData[ivol], Depth_level = Profile.depth, Interpolate = true, dims = DimsVolCross) # create a horizontal cross section + names_fields = String.(keys(cross_tmp.fields)) + for ifield in eachindex(names_fields) + name_new_field = datasetnames[ivol] * "_" * names_fields[ifield] # name of new field includes name of dataset + cross_add = addfield(cross_add, name_new_field, ustrip.(cross_tmp.fields[ifield])) # Note: we use ustrip here, and thereby remove the units, as the cross-section routine made problems end end end From e0e087fa2e1e4755995977f05b46c0a13f80757d Mon Sep 17 00:00:00 2001 From: Marcel Thielmann Date: Thu, 24 Sep 2026 08:44:24 +0200 Subject: [PATCH 17/25] deleted gmt.history --- .gitignore | 1 + test/gmt.history | 4 ---- 2 files changed, 1 insertion(+), 4 deletions(-) delete mode 100644 test/gmt.history diff --git a/.gitignore b/.gitignore index 91c96fe2b..9670f7a89 100644 --- a/.gitignore +++ b/.gitignore @@ -27,3 +27,4 @@ tutorial/test.tiff *.tiff *.nc !/test/test_files/ISCTest.xml +test/gmt.history diff --git a/test/gmt.history b/test/gmt.history deleted file mode 100644 index 9ffd9e24a..000000000 --- a/test/gmt.history +++ /dev/null @@ -1,4 +0,0 @@ -# GMT 6 Session common arguments shelf -BEGIN GMT 6.8.0 -R 341.8/342.5/28.4/29.0 -END From 6aca216cb36b1ae38397aa0931e7f6af445414eb Mon Sep 17 00:00:00 2001 From: Marcel Thielmann Date: Fri, 2 Oct 2026 10:10:25 +0200 Subject: [PATCH 18/25] modified read_picked_profiles to allow for the import of horizontal profile specifications --- src/ProfileProcessing.jl | 25 +++++++++++++++++++------ 1 file changed, 19 insertions(+), 6 deletions(-) diff --git a/src/ProfileProcessing.jl b/src/ProfileProcessing.jl index 94d1f0500..b28495462 100644 --- a/src/ProfileProcessing.jl +++ b/src/ProfileProcessing.jl @@ -562,13 +562,26 @@ function read_picked_profiles(ProfileCoordFile::String) profiles = Vector{ProfileData}() profile_data = readdlm(ProfileCoordFile, skipstart = 1, ',') - for i in 1:size(profile_data, 1) - start_lonlat = (profile_data[i, 2:3]...,) - end_lonlat = (profile_data[i, 4:5]...,) - profile = ProfileData(start_lonlat = start_lonlat, end_lonlat = end_lonlat) - push!(profiles, profile) + # check the number of columns in the profile data + if size(profile_data, 2) == 2 + # we have a horizontal profile, so we only have the depth in the second column + for i in 1:size(profile_data, 1) + depth = profile_data[i, 2] + profile = ProfileData(depth = depth) + push!(profiles, profile) + end + elseif size(profile_data, 2) == 5 + # we have a vertical profile, so we have the start and end lon/lat in + for i in 1:size(profile_data, 1) + start_lonlat = (profile_data[i, 2:3]...,) + end_lonlat = (profile_data[i, 4:5]...,) + profile = ProfileData(start_lonlat = start_lonlat, end_lonlat = end_lonlat) + push!(profiles, profile) + end + else + error("ProfileCoordFile should have either 2 columns (for horizontal profiles) or 5 columns (for vertical profiles).") end - + return profiles end From c215f0612c0997633be09342bddb06f7e803c1d6 Mon Sep 17 00:00:00 2001 From: Marcel Thielmann Date: Fri, 2 Oct 2026 12:14:15 +0200 Subject: [PATCH 19/25] additional modification to read_picked_profiles to allow for mixed profile files --- src/ProfileProcessing.jl | 23 ++++++++++++++++------- 1 file changed, 16 insertions(+), 7 deletions(-) diff --git a/src/ProfileProcessing.jl b/src/ProfileProcessing.jl index b28495462..c4795b774 100644 --- a/src/ProfileProcessing.jl +++ b/src/ProfileProcessing.jl @@ -564,24 +564,33 @@ function read_picked_profiles(ProfileCoordFile::String) # check the number of columns in the profile data if size(profile_data, 2) == 2 - # we have a horizontal profile, so we only have the depth in the second column + # we have horizontal profiles only, so we only have the depth in the second column for i in 1:size(profile_data, 1) depth = profile_data[i, 2] profile = ProfileData(depth = depth) push!(profiles, profile) end elseif size(profile_data, 2) == 5 - # we have a vertical profile, so we have the start and end lon/lat in + # we have mixed horizontal/vertical profiles or vertical profiles only for i in 1:size(profile_data, 1) - start_lonlat = (profile_data[i, 2:3]...,) - end_lonlat = (profile_data[i, 4:5]...,) - profile = ProfileData(start_lonlat = start_lonlat, end_lonlat = end_lonlat) - push!(profiles, profile) + # check whether the current profile is horizontal or vertical + if all(isempty.(bla[1,3:5])) # there are only two entries: profile number and depth of the profile, so this is a horizontal profile + depth = profile_data[i, 2] + profile = ProfileData(depth = depth) + push!(profiles, profile) + elseif !all(isempty.(bla[7,1:5])) # there are five entries: profile number, start lon/lat, end lon/lat, so this is a vertical profile + start_lonlat = (profile_data[i, 2:3]...,) + end_lonlat = (profile_data[i, 4:5]...,) + profile = ProfileData(start_lonlat = start_lonlat, end_lonlat = end_lonlat) + push!(profiles, profile) + else + error("ProfileCoordFile should have either 2 columns (for horizontal profiles) or 5 columns (for vertical profiles).") + end end else error("ProfileCoordFile should have either 2 columns (for horizontal profiles) or 5 columns (for vertical profiles).") end - + return profiles end From 557e1c2aa36bceb52df8a494207e60d808013b45 Mon Sep 17 00:00:00 2001 From: Marcel Thielmann Date: Fri, 2 Oct 2026 14:40:51 +0200 Subject: [PATCH 20/25] fixed the issue in create_profile_screenshot --- src/ProfileProcessing.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/ProfileProcessing.jl b/src/ProfileProcessing.jl index c4795b774..ae7c609e3 100644 --- a/src/ProfileProcessing.jl +++ b/src/ProfileProcessing.jl @@ -363,9 +363,9 @@ function create_profile_screenshot!(Profile::ProfileData, DataSet::NamedTuple) tmp = NamedTuple() # initialize empty one DataSetName = keys(DataSet) # Names of the datasets - for idata in 1:num_datasets + for idata in eachindex(DataSetName) # 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] + data_tmp = DataSet[DataSetName[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) From bc05e729ba7ef7e1c92c03db27b34b22acd6a3d0 Mon Sep 17 00:00:00 2001 From: Marcel Thielmann Date: Fri, 2 Oct 2026 14:57:30 +0200 Subject: [PATCH 21/25] Trying to fix the error in the test --- test/test_ProfileProcessing.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_ProfileProcessing.jl b/test/test_ProfileProcessing.jl index 040198d8a..d4806c755 100644 --- a/test/test_ProfileProcessing.jl +++ b/test/test_ProfileProcessing.jl @@ -110,7 +110,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) +extract_ProfileData!(prof5, VolData_combined3, Data.Surface, Data.Point;ScreenshotData=Data.Screenshot) # Test that it works if only EQ's are provided: From fd0d085df6de0a704ee490cc6662509600ac493a Mon Sep 17 00:00:00 2001 From: Marcel Thielmann Date: Fri, 2 Oct 2026 23:43:14 +0200 Subject: [PATCH 22/25] Accept nothing/empty tuples in extract_ProfileData!, fix read_picked_profiles Co-Authored-By: Claude Sonnet 5.5 --- src/ProfileProcessing.jl | 56 +++++++++------------------------- test/test_ProfileProcessing.jl | 10 ++++++ 2 files changed, 24 insertions(+), 42 deletions(-) diff --git a/src/ProfileProcessing.jl b/src/ProfileProcessing.jl index ae7c609e3..b9d9347aa 100644 --- a/src/ProfileProcessing.jl +++ b/src/ProfileProcessing.jl @@ -496,59 +496,31 @@ end """ - extract_ProfileData!(Profile::ProfileData,VolData::NamedTuple, SurfData::NamedTuple, PointData::NamedTuple; DimsVolCross=(100,100),Depth_extent=nothing,DimsSurfCross=(100,),section_width=50, ScreenshotData=nothing, TopoData=NamedTuple()) + extract_ProfileData!(Profile::ProfileData, VolData=nothing, SurfData=nothing, PointData=nothing; DimsVolCross=(100,100), Depth_extent=nothing, DimsSurfCross=(100,), section_width=50km, ScreenshotData=nothing, TopoData=nothing) -Extracts data along a vertical or horizontal profile. Allows VolData to be passed as a NamedTuple. +Extracts data along a vertical or horizontal profile. `VolData` can be a `GeoData`/`AbstractGeneralGrid` or a `NamedTuple` of those; the other data sets are `NamedTuple`s. +Any data argument can be `nothing` or empty (`NamedTuple()` or `()`), in which case it is skipped. """ -function extract_ProfileData!(Profile::ProfileData, VolData::NamedTuple = NamedTuple(), SurfData::NamedTuple = NamedTuple(), PointData::NamedTuple = NamedTuple(); DimsVolCross = (100, 100), Depth_extent = nothing, DimsSurfCross = (100,), section_width = 50km, ScreenshotData = nothing, TopoData::NamedTuple = NamedTuple()) +function extract_ProfileData!(Profile::ProfileData, VolData = nothing, SurfData = nothing, PointData = nothing; DimsVolCross = (100, 100), Depth_extent = nothing, DimsSurfCross = (100,), section_width = 50km, ScreenshotData = nothing, TopoData = nothing) - return extract_ProfileData!(Profile, VolData, SurfData, PointData, TopoData, ScreenshotData; DimsVolCross = DimsVolCross, Depth_extent = Depth_extent, DimsSurfCross = DimsSurfCross, section_width = section_width) -end - -# Internal method - called by the main method with ScreenshotData as positional argument, allows VolData as NamedTuple -function extract_ProfileData!(Profile::ProfileData, VolData::NamedTuple, SurfData::NamedTuple, PointData::NamedTuple, TopoData:: Union{Nothing, NamedTuple},ScreenshotData::Union{Nothing, NamedTuple}; DimsVolCross = (100, 100), Depth_extent = nothing, DimsSurfCross = (100,), section_width = 50km) - - if !isempty(VolData) + if !_isnodata(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 !isnothing(TopoData) + create_profile_surface!(Profile, _as_namedtuple(SurfData), DimsSurfCross = DimsSurfCross) + create_profile_point!(Profile, _as_namedtuple(PointData), section_width = section_width) + if !_isnodata(TopoData) create_profile_topo!(Profile, TopoData, DimsSurfCross = 5 .* DimsSurfCross) # we use a larger number of points for the topography, as it is often more detailed than the surface data end - if !isnothing(ScreenshotData) + if !_isnodata(ScreenshotData) create_profile_screenshot!(Profile, ScreenshotData) end return nothing end +# helpers: treat `nothing` and empty containers (NamedTuple(), ()) as "no data" +_isnodata(x) = isnothing(x) || (x isa Union{Tuple, NamedTuple} && isempty(x)) +_as_namedtuple(x) = _isnodata(x) ? NamedTuple() : x -""" - extract_ProfileData!(Profile::ProfileData,VolData::GeoData, SurfData::NamedTuple, PointData::NamedTuple; DimsVolCross=(100,100),Depth_extent=nothing,DimsSurfCross=(100,),section_width=50, ScreenshotData=nothing, TopoData=NamedTuple()) - -Extracts data along a vertical or horizontal profile -""" -function extract_ProfileData!(Profile::ProfileData, VolData::GeoData, SurfData::NamedTuple = NamedTuple(), PointData::NamedTuple = NamedTuple(); DimsVolCross = (100, 100), Depth_extent = nothing, DimsSurfCross = (100,), section_width = 50km, ScreenshotData = nothing, TopoData::NamedTuple = NamedTuple()) - - return extract_ProfileData!(Profile, VolData, SurfData, PointData, TopoData, ScreenshotData; DimsVolCross = DimsVolCross, Depth_extent = Depth_extent, DimsSurfCross = DimsSurfCross, section_width = section_width) -end - -# Internal method - called by the main method with ScreenshotData as positional argument -function extract_ProfileData!(Profile::ProfileData, VolData::GeoData, SurfData::NamedTuple, PointData::NamedTuple, TopoData::Union{Nothing, NamedTuple}, ScreenshotData::Union{Nothing, 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 !isnothing(TopoData) - create_profile_topo!(Profile, TopoData, DimsSurfCross = 5 .* DimsSurfCross) # we use a larger number of points for the topography, as it is often more detailed than the surface data - end - if !isnothing(ScreenshotData) - create_profile_screenshot!(Profile, ScreenshotData) - end - return nothing -end @@ -574,11 +546,11 @@ function read_picked_profiles(ProfileCoordFile::String) # we have mixed horizontal/vertical profiles or vertical profiles only for i in 1:size(profile_data, 1) # check whether the current profile is horizontal or vertical - if all(isempty.(bla[1,3:5])) # there are only two entries: profile number and depth of the profile, so this is a horizontal profile + if all(isempty.(profile_data[i,3:5])) # there are only two entries: profile number and depth of the profile, so this is a horizontal profile depth = profile_data[i, 2] profile = ProfileData(depth = depth) push!(profiles, profile) - elseif !all(isempty.(bla[7,1:5])) # there are five entries: profile number, start lon/lat, end lon/lat, so this is a vertical profile + elseif !any(isempty.(profile_data[i,2:5])) # there are five entries: profile number, start lon/lat, end lon/lat, so this is a vertical profile start_lonlat = (profile_data[i, 2:3]...,) end_lonlat = (profile_data[i, 4:5]...,) profile = ProfileData(start_lonlat = start_lonlat, end_lonlat = end_lonlat) diff --git a/test/test_ProfileProcessing.jl b/test/test_ProfileProcessing.jl index d4806c755..5d2c5293f 100644 --- a/test/test_ProfileProcessing.jl +++ b/test/test_ProfileProcessing.jl @@ -120,6 +120,16 @@ extract_ProfileData!(prof4, nothing, NamedTuple(), Data.Point) @test isempty(prof4.SurfData) @test length(prof4.PointData[1].depth) == 3280 +# nothing / empty tuples are accepted for all data arguments +prof6 = ProfileData(depth = -20) +extract_ProfileData!(prof6, (), nothing, Data.Point; TopoData=(), ScreenshotData=nothing) +@test isnothing(prof6.VolData) +@test isempty(prof6.SurfData) +@test length(prof6.PointData[1].depth) == 3280 +prof7 = ProfileData(depth = -20) +extract_ProfileData!(prof7) +@test isempty(prof7.SurfData) + @test prof1.SurfData[1].fields[1][80] ≈ -37.58791461075397km @test isempty(prof2.SurfData) @test isnan(prof3.SurfData[1].fields[1][80]) From 27dcf4f60a611baf01e6498d43964cfea5d39f4e Mon Sep 17 00:00:00 2001 From: Marcel Thielmann Date: Mon, 5 Oct 2026 09:05:21 +0200 Subject: [PATCH 23/25] changed handling of input arguments in profile processing, fixed a bug --- src/ProfileProcessing.jl | 22 ++++++++++++++++++---- test/test_ProfileProcessing.jl | 13 ++++++++++++- 2 files changed, 30 insertions(+), 5 deletions(-) diff --git a/src/ProfileProcessing.jl b/src/ProfileProcessing.jl index b9d9347aa..28a0a65cd 100644 --- a/src/ProfileProcessing.jl +++ b/src/ProfileProcessing.jl @@ -299,6 +299,11 @@ Creates a cross-section through a volumetric 3D dataset `VolData` with the data """ function create_profile_volume!(Profile::ProfileData, VolData::NamedTuple; DimsVolCross::NTuple = (100, 100), Depth_extent = nothing) + _isnodata(VolData) && return nothing # empty NamedTuple: nothing to intersect + for (name, data) in pairs(VolData) + data isa AbstractGeneralGrid || throw(ArgumentError("VolData.$name must be a GeoData/AbstractGeneralGrid, got $(typeof(data))")) + end + datasetnames = String.(keys(VolData)) # get the names of the datasets if Profile.vertical # take a vertical cross section @@ -356,6 +361,17 @@ function create_profile_volume!(Profile::ProfileData, VolData::NamedTuple; DimsV return nothing end +""" + create_profile_volume!(Profile::ProfileData, VolData; DimsVolCross::NTuple=(100,100), Depth_extent=nothing) + +Fallback that normalizes `VolData`: `nothing` and empty containers (`NamedTuple()`, `()`) are treated as "no volume data" and leave `Profile` unchanged. +Any other unsupported type throws an `ArgumentError`. +""" +function create_profile_volume!(Profile::ProfileData, VolData; DimsVolCross::NTuple = (100, 100), Depth_extent = nothing) + _isnodata(VolData) && return nothing + throw(ArgumentError("VolData must be a GeoData/AbstractGeneralGrid, a NamedTuple of those, or nothing; got $(typeof(VolData))")) +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) @@ -378,7 +394,7 @@ function create_profile_screenshot!(Profile::ProfileData, DataSet::NamedTuple) #error("horizontal profiles not yet implemented") end end - Profile.SurfData = tmp # assign to profile data structure + Profile.ScreenshotData = tmp # assign to profile data structure return end @@ -503,9 +519,7 @@ Any data argument can be `nothing` or empty (`NamedTuple()` or `()`), in which c """ function extract_ProfileData!(Profile::ProfileData, VolData = nothing, SurfData = nothing, PointData = nothing; DimsVolCross = (100, 100), Depth_extent = nothing, DimsSurfCross = (100,), section_width = 50km, ScreenshotData = nothing, TopoData = nothing) - if !_isnodata(VolData) - create_profile_volume!(Profile, VolData; DimsVolCross = DimsVolCross, Depth_extent = Depth_extent) - end + create_profile_volume!(Profile, VolData; DimsVolCross = DimsVolCross, Depth_extent = Depth_extent) create_profile_surface!(Profile, _as_namedtuple(SurfData), DimsSurfCross = DimsSurfCross) create_profile_point!(Profile, _as_namedtuple(PointData), section_width = section_width) if !_isnodata(TopoData) diff --git a/test/test_ProfileProcessing.jl b/test/test_ProfileProcessing.jl index 5d2c5293f..56fff6379 100644 --- a/test/test_ProfileProcessing.jl +++ b/test/test_ProfileProcessing.jl @@ -92,7 +92,7 @@ GeophysicalModelGenerator.create_profile_point!(prof4, Data.Point, section_width # test screenshot data GeophysicalModelGenerator.create_profile_screenshot!(prof5, Data.Screenshot) -@test prof5.SurfData[1].fields.x_profile[1,1,1] == 0 +@test prof5.ScreenshotData[1].fields.x_profile[1,1,1] == 0 # Test the main profile extraction routines: @@ -130,6 +130,17 @@ prof7 = ProfileData(depth = -20) extract_ProfileData!(prof7) @test isempty(prof7.SurfData) +# create_profile_volume! normalizes VolData: nothing / empty containers are no-ops, bad types error +prof8 = ProfileData(depth = -100) +for novol in (nothing, NamedTuple(), ()) + GeophysicalModelGenerator.create_profile_volume!(prof8, novol) + @test isnothing(prof8.VolData) +end +GeophysicalModelGenerator.create_profile_volume!(prof8, (Hua2017 = Data.Volume[1],)) +@test haskey(prof8.VolData.fields, :Hua2017_Vp) +@test_throws ArgumentError GeophysicalModelGenerator.create_profile_volume!(prof8, 42) +@test_throws ArgumentError GeophysicalModelGenerator.create_profile_volume!(prof8, (a = 42,)) + @test prof1.SurfData[1].fields[1][80] ≈ -37.58791461075397km @test isempty(prof2.SurfData) @test isnan(prof3.SurfData[1].fields[1][80]) From 794266791b4f09a9eae0fc170823fd8f54140a6a Mon Sep 17 00:00:00 2001 From: Marcel Thielmann Date: Mon, 5 Oct 2026 09:59:10 +0200 Subject: [PATCH 24/25] Added a tutorial page on how to create profiles for the AdriaArrayGeometryPicker --- docs/make.jl | 3 +- .../man/tutorial_AdriaArrayGeometryPicker.md | 213 ++++++++++++++++++ src/ProfileProcessing.jl | 2 +- 3 files changed, 216 insertions(+), 2 deletions(-) create mode 100644 docs/src/man/tutorial_AdriaArrayGeometryPicker.md diff --git a/docs/make.jl b/docs/make.jl index 154676eed..05d38c616 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -106,7 +106,8 @@ makedocs(; "20 - 2D model setups" => "man/Tutorial_NumericalModel_2D.md", "21 - 3D model setups" => "man/Tutorial_NumericalModel_3D.md", "22 - 3D Volcano setup" => "man/Tutorial_VolcanoModel_3D.md", - "23 - Build geometry from polygons" => "man/tutorial_Polygon_structures.md" + "23 - Build geometry from polygons" => "man/tutorial_Polygon_structures.md", + "24 - Profiles for the AdriaArrayGeometryPicker" => "man/tutorial_AdriaArrayGeometryPicker.md" ], "User Guide" => Any[ "Installation" => "man/installation.md", diff --git a/docs/src/man/tutorial_AdriaArrayGeometryPicker.md b/docs/src/man/tutorial_AdriaArrayGeometryPicker.md new file mode 100644 index 000000000..c74c47431 --- /dev/null +++ b/docs/src/man/tutorial_AdriaArrayGeometryPicker.md @@ -0,0 +1,213 @@ +# Create profiles for the AdriaArrayGeometryPicker + +## Goal +The [AdriaArrayGeometryPicker](https://github.com/JuliaGeodynamics/AdriaArrayGeometryPicker.jl) is a graphical user interface to compare different geophysical datasets along profiles and to pick structures (e.g., slabs or the Moho) in them. It does not read the original datasets itself, but works with *profiles* that have been created beforehand with GeophysicalModelGenerator: every dataset is projected onto the profile, and the result is saved as a `.pgmg` file that is then opened in the picker. + +This tutorial shows how to create these `.pgmg` files for a whole set of vertical and horizontal profiles at once. You need two text files: + +1. a **profile file** that lists the profiles (e.g. `profiles_mixed.txt`), +2. a **dataset file** that lists the datasets to be projected onto the profiles (e.g. `Datasets_Adria.txt`). + +Both are described in detail below. Some general background on profile processing is given in the [Profile Processing](profile_processing.md) section. + +!!! note + The AdriaArrayGeometryPicker currently relies on the `profile_processing_mt` branch of GeophysicalModelGenerator. If you install the picker, this branch is installed automatically. + +## 1. The profile file +The profile file is a comma-separated text file. The first line is a header and is skipped. Every following line defines one profile, which can be either horizontal or vertical: + +| Profile type | Columns | Meaning | +|:--|:--|:--| +| horizontal | `Num, Depth` | profile number and depth of the horizontal slice in km (negative below the surface) | +| vertical | `Num, LonStart, LatStart, LonEnd, LatEnd` | profile number and the longitude/latitude (in degrees) of the start and end point | + +Horizontal and vertical profiles can be mixed in one file. This is the file `profiles_mixed.txt`, which contains two horizontal slices (at 200 and 300 km depth) and five vertical profiles across the Adriatic region: + +``` +Num,Depth +1,-200 +2,-300 +3,9.777223439354351,42.936916884179524,15.465423125301268,47.79656290565586 +4,10.160585214657562,42.71827704599679,15.860636558746865,47.55888190360081 +5,10.541291181660043,42.49821887191029,16.252277335405516,47.31986389797744 +6,10.919360031561707,42.27677067866074,16.64038623036512,47.07954242435165 +7,11.294811014000167,42.05396063315099,17.0250043232569,46.837950655777426 +``` + +A few remarks: +- The header is only skipped, its content does not matter. For a file with only vertical profiles you can, e.g., use `Num,LonStart,LatStart,LonEnd,LatEnd`. +- A line with one value after `Num` is interpreted as a horizontal profile, a line with four values as a vertical profile. Any other number of values results in an error. +- The profiles are processed in the order in which they appear in the file; `Num` itself is not used. +- A vertical profile starts at (`LonStart`,`LatStart`). The distance along the profile (`x_profile`, in km) is measured from this point. + +The file is read with `read_picked_profiles`, which returns a vector of `ProfileData` objects: +```julia +julia> using GeophysicalModelGenerator +julia> profiles = read_picked_profiles("profiles_mixed.txt"); +julia> profiles[1] +Horizontal ProfileData + depth : -200.0 +julia> profiles[3] +Vertical ProfileData + lon/lat : (9.777223439354351, 42.936916884179524)-(15.465423125301268, 47.79656290565586) +``` + +## 2. The dataset file +The dataset file is also a comma-separated text file, with one line per dataset. Again, the first line is a header and is skipped. Each line has the columns + +| Column | Meaning | +|:--|:--| +| `Name` | Name of the dataset. It appears in the picker and is used as prefix of the field names in the profile (e.g. `Giacomuzzi_dVp_perc`). Avoid spaces and special characters. | +| `Location` | Location of the dataset: a local `*.jld2` file (absolute path, or relative to the directory in which you run julia), or a url starting with `http` from which the file is downloaded. The `.jld2` file must contain a GMG data structure, such as one saved with `save_GMG`. | +| `Type` | One of `Volume`, `Surface`, `Point`, `Topography` or `Screenshot` (see below). | +| `Active` | Optional, `true` or `false`. Inactive datasets are not loaded. If the column is missing or empty, `true` is assumed. | + +The meaning of the different types is: +- `Volume`: 3D data, such as seismic tomography models. A cross-section is computed through every volume dataset. +- `Surface`: surfaces such as the Moho depth. For vertical profiles, the intersection of the surface with the profile is computed, which is shown as a line. +- `Point`: point data, such as earthquake hypocenters. All points within a band of width `section_width` around the profile are projected onto it. +- `Topography`: the topography of the region. It is shown above vertical profiles and in the map overview of the picker. The extent of the (first) topography dataset also determines the resolution of the profiles in the script below, so it should cover the region of interest. +- `Screenshot`: screenshots of figures from publications (see [Import screenshots](tutorial_Screenshot_To_Paraview.md)). These are not used in this tutorial. + +This is the file `Datasets_Adria.txt` used in this tutorial (with the paths shortened to `DataSet/...`, a directory relative to where julia is started): +``` +Name,Location,Type, [Active] +Giacomuzzi,DataSet/Giacomuzzi.jld2,Volume,true +Timko2023,DataSet/CaPaREA2023_Mantle_Timko.jld2,Volume,true +Friederich2025,DataSet/Friederich2025_Pwave_Alps.jld2,Volume,true +Koulakov,DataSet/Koulakov_Europe.jld2,Volume,true +ElSharkawy2020,DataSet/MeRE2020_El-Sharkawy.jld2,Volume,true +MIT08_Li,DataSet/MIT08_Pwave.jld2,Volume,true +Piromallo2003,DataSet/Piromallo2003.jld2,Volume,true +REVEAL,DataSet/REVEAL_rel_avg_AK135.jld2,Volume,true +SAVANI,DataSet/SAVANI_Auer2014.jld2,Volume,true +UU07_Amaru,DataSet/UU07_Pwave.jld2,Volume,true +Zhu2015,DataSet/Zhu2015.jld2,Volume,true +ETOPO1,DataSet/etopo1.jld2,Topography,true +Grad2009,DataSet/Grad2009_EU_Amr.jld2,Surface,true +ISC,DataSet/isc24_MedS.jld2,Point,true +``` +It contains eleven tomographic models, the ETOPO1 topography, the Moho depth of Grad et al. (2009) and earthquake locations from the ISC catalogue. To exclude a dataset temporarily (for example, a large one while you test your profiles), set its last column to `false`. + +!!! warning + Since the file is comma-separated, the paths and urls must not contain commas. + +The dataset file is read with `load_dataset_file`, and `load_GMG` then loads all active datasets and sorts them by type: +```julia +julia> datasets = load_dataset_file("Datasets_Adria.txt"); +julia> data = load_GMG(datasets); +julia> keys(data) +(:Volume, :Surface, :Point, :Screenshot, :Topography) +julia> keys(data.Volume) +(:Giacomuzzi, :Timko2023, :Friederich2025, :Koulakov, :ElSharkawy2020, :MIT08_Li, :Piromallo2003, :REVEAL, :SAVANI, :UU07_Amaru, :Zhu2015) +``` +Every entry of `data` is a `NamedTuple` with the datasets of that type, named as in the `Name` column. + +## 3. Create the profiles +The following script projects all datasets onto all profiles and saves every profile as a `.pgmg` file. Besides GeophysicalModelGenerator, it needs the [JLD2](https://github.com/JuliaIO/JLD2.jl) package (`] add JLD2`). + +```julia +using GeophysicalModelGenerator, JLD2 + +file_profiles = "profiles_mixed.txt" +file_datasets = "Datasets_Adria.txt" + +# read the profiles and load all active datasets +profiles = read_picked_profiles(file_profiles) +datasets = load_dataset_file(file_datasets) +data = load_GMG(datasets) +``` + +Next, we determine the resolution of the profiles. We take the lon/lat extent of the topography and the finest grid spacing of all volume datasets, so that no information of the highest-resolution model is lost: +```julia +lon_range = extrema(data.Topography[1].lon.val) +lat_range = extrema(data.Topography[1].lat.val) + +res_lon = minimum(minimum(diff(vol.lon.val, dims = 1)) for vol in data.Volume) +res_lat = minimum(minimum(diff(vol.lat.val, dims = 2)) for vol in data.Volume) + +nlon = round(Int, (lon_range[2] - lon_range[1]) / res_lon) +nlat = round(Int, (lat_range[2] - lat_range[1]) / res_lat) +``` + +Now we loop over the profiles, project the data onto every profile with `extract_ProfileData!`, and save the result: +```julia +for (i, profile) in enumerate(profiles) + if profile.vertical + # vertical profile: (points along the profile, points in depth) + dims_vol = (nlon, 200) + prefix = "Profile_vertical" + else + # horizontal profile: (points in lon, points in lat) + dims_vol = (nlon, nlat) + prefix = "Profile_horizontal" + end + + extract_ProfileData!(profile, data.Volume, data.Surface, data.Point; + TopoData = data.Topography, + DimsVolCross = dims_vol, + Depth_extent = (-400, 0), + DimsSurfCross = (nlon,), + section_width = 50km) + + # surfaces and topography are only intersected with vertical profiles, + # so for horizontal profiles we add the full datasets for the map view + if !profile.vertical + profile.SurfData = data.Surface + profile.TopoData = data.Topography + end + + # save the profile for the picker + jldsave("$(prefix)$(i).pgmg"; profile) +end +``` + +The options of `extract_ProfileData!` are: +- `DimsVolCross`: number of grid points of the cross-sections through the volume data (see above). +- `Depth_extent`: minimum and maximum depth (in km) of vertical profiles. Here, we use the uppermost 400 km. +- `DimsSurfCross`: number of points along the profile at which surfaces are intersected. The topography is sampled with 5 times as many points. +- `section_width`: width of the band around the profile from which point data (earthquakes) are projected onto the profile. + +You can also call `extract_ProfileData!` with `nothing` (or an empty `NamedTuple()`) for any dataset type that you do not have, e.g. `extract_ProfileData!(profile, data.Volume, nothing, data.Point)`. + +After running the script, you have one file per profile: `Profile_horizontal1.pgmg`, `Profile_horizontal2.pgmg` and `Profile_vertical3.pgmg` to `Profile_vertical7.pgmg`. A processed vertical profile contains: +```julia +julia> profiles[3] +Vertical ProfileData + lon/lat : (9.777223439354351, 42.936916884179524)-(15.465423125301268, 47.79656290565586) + VolData : (:x_profile, :Giacomuzzi_dVp_perc, :Giacomuzzi_dVs_perc, ..., :Piromallo2003_Vp, :Piromallo2003_dVp_perc, ...) + SurfData : (:Grad2009,) + TopoData : (:ETOPO1,) + PointData: (:ISC,) +``` +All volume datasets are combined into a single `GeoData` structure, in which the name of each field is prefixed with the name of its dataset. `x_profile` is the distance along the profile in km. + +!!! tip + Processing many large tomographic models can take a while and needs a fair amount of memory, as all datasets are loaded at the same time. Test your profiles first with only a few active datasets. + +## 4. The `.pgmg` file format +A `.pgmg` file is a [JLD2](https://github.com/JuliaIO/JLD2.jl) file that contains a single `ProfileData` object; the name under which it is stored does not matter. The picker uses the following fields: + +| Field | Content | +|:--|:--| +| `vertical` | `true` for a vertical profile, `false` for a horizontal slice | +| `start_lonlat`, `end_lonlat` | start and end point of a vertical profile | +| `depth` | depth of a horizontal slice | +| `VolData` | the cross-sections through all volume datasets, combined in one `GeoData` | +| `SurfData` | `NamedTuple` of `GeoData`, one per surface dataset | +| `PointData` | `NamedTuple` of `GeoData`, one per point dataset | +| `TopoData` | `NamedTuple` of `GeoData`, one per topography dataset | + +You can load a profile back into julia with +```julia +julia> using JLD2 +julia> profile = load("Profile_vertical3.pgmg", "profile"); +``` + +## 5. Open the profiles in the AdriaArrayGeometryPicker +Install the picker as described in its [README](https://github.com/JuliaGeodynamics/AdriaArrayGeometryPicker.jl), and open a profile with +```julia +julia> using AdriaArrayGeometryPicker +julia> geometry_picker("Profile_vertical3.pgmg") +``` +Alternatively, start `geometry_picker()` without arguments and use **Load Profile** in the menu. How to pick structures and save the picks is described in the documentation of the picker. diff --git a/src/ProfileProcessing.jl b/src/ProfileProcessing.jl index 28a0a65cd..0764c3cf4 100644 --- a/src/ProfileProcessing.jl +++ b/src/ProfileProcessing.jl @@ -35,7 +35,7 @@ mutable struct ProfileData 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) + K = new(true, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing) for (key, value) in kwargs # make sure that start and end point are given as tuples of Float64 if key == Symbol("start_lonlat") From d1b634a8037acbf442e674ee051967083b8e3453 Mon Sep 17 00:00:00 2001 From: Marcel Thielmann Date: Mon, 5 Oct 2026 11:41:56 +0200 Subject: [PATCH 25/25] set rendering of the documentation in pull requests to true --- docs/make.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/make.jl b/docs/make.jl index 05d38c616..13c1ac5fd 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -146,5 +146,5 @@ deploydocs(; devbranch = "main", devurl = "dev", forcepush=true, - push_preview = false, + push_preview = true, )