diff --git a/src/ProfileProcessing.jl b/src/ProfileProcessing.jl index 6fe0aff9..94d1f050 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 @@ -264,7 +269,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) @@ -285,6 +291,71 @@ 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 + + # 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 + + 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 + + # 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 + + # 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 + 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) @@ -345,6 +416,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) @@ -390,23 +496,54 @@ 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::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. +""" +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, 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) + 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 + + +""" + 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::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, 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::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, 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 @@ -414,6 +551,9 @@ function extract_ProfileData!(Profile::ProfileData, VolData::Union{Nothing, GeoD end + + + """ This reads the picked profiles from disk and returns a vector of ProfileData """ diff --git a/test/gmt.history b/test/gmt.history index ebcde871..9ffd9e24 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 0a27c173..040198d8 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