Skip to content
Merged
Show file tree
Hide file tree
Changes from 1 commit
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
34 changes: 31 additions & 3 deletions src/ProfileProcessing.jl
Original file line number Diff line number Diff line change
Expand Up @@ -29,6 +29,7 @@
VolData::Union{Nothing, GeophysicalModelGenerator.GeoData}
SurfData::Union{Nothing, NamedTuple}
PointData::Union{Nothing, NamedTuple}
ScreenshotData::Union{Nothing, NamedTuple}

function ProfileData(; kwargs...) # this constructor allows to define only certain fields and leave the others blank
K = new(true, nothing, nothing, nothing, nothing, nothing, nothing)
Expand Down Expand Up @@ -282,6 +283,31 @@
return nothing
end

### internal function to process screenshot data - contrary to the volume data, we here have to save lon/lat/depth pairs for every screenshot data set, so we create a NamedTuple of GeoData data sets
function create_profile_screenshot!(Profile::ProfileData, DataSet::NamedTuple)
num_datasets = length(DataSet)

tmp = NamedTuple() # initialize empty one
DataSetName = keys(DataSet) # Names of the datasets

for idata in 1:num_datasets
# load data set --> each data set is a single GeoData structure, so we'll only have to get the respective key to load the correct type
data_tmp = DataSet[idata]
if Profile.vertical
x_profile = flatten_cross_section(data_tmp, Start = Profile.start_lonlat) # compute the distance along the profile
data_tmp = addfield(data_tmp, "x_profile", x_profile)

# add the data set as a NamedTuple
data_NT = NamedTuple{(DataSetName[idata],)}((data_tmp,))
tmp = merge(tmp, data_NT)
else

Check warning on line 303 in src/ProfileProcessing.jl

View check run for this annotation

Codecov / codecov/patch

src/ProfileProcessing.jl#L303

Added line #L303 was not covered by tests
# we do not have this implemented
#error("horizontal profiles not yet implemented")
end

Check warning on line 306 in src/ProfileProcessing.jl

View check run for this annotation

Codecov / codecov/patch

src/ProfileProcessing.jl#L306

Added line #L306 was not covered by tests
end
Profile.SurfData = tmp # assign to profile data structure
return
end

### internal function to process surface data - contrary to the volume data, we here have to save lon/lat/depth pairs for every surface data set, so we create a NamedTuple of GeoData data sets
function create_profile_surface!(Profile::ProfileData, DataSet::NamedTuple; DimsSurfCross = (100,))
Expand Down Expand Up @@ -362,18 +388,20 @@


"""
extract_ProfileData!(Profile::ProfileData,VolData::GeoData, SurfData::NamedTuple, PointData::NamedTuple; DimsVolCross=(100,100),Depth_extent=nothing,DimsSurfCross=(100,),section_width=50)
extract_ProfileData!(Profile::ProfileData,VolData::GeoData, SurfData::NamedTuple, PointData::NamedTuple, ScreenshotData::NamedTuple; DimsVolCross=(100,100),Depth_extent=nothing,DimsSurfCross=(100,),section_width=50)
Extracts data along a vertical or horizontal profile
"""
function extract_ProfileData!(Profile::ProfileData, VolData::Union{Nothing, GeoData} = nothing, SurfData::NamedTuple = NamedTuple(), PointData::NamedTuple = NamedTuple(); DimsVolCross = (100, 100), Depth_extent = nothing, DimsSurfCross = (100,), section_width = 50km)
function extract_ProfileData!(Profile::ProfileData, VolData::Union{Nothing, GeoData} = nothing, SurfData::NamedTuple = NamedTuple(), PointData::NamedTuple = NamedTuple(),ScreenshotData::NamedTuple = NamedTuple(); DimsVolCross = (100, 100), Depth_extent = nothing, DimsSurfCross = (100,), section_width = 50km)

if !isnothing(VolData)
create_profile_volume!(Profile, VolData; DimsVolCross = DimsVolCross, Depth_extent = Depth_extent)
end
create_profile_surface!(Profile, SurfData, DimsSurfCross = DimsSurfCross)
create_profile_point!(Profile, PointData, section_width = section_width)

if !isempty(ScreenshotData)
create_profile_screenshot!(Profile, ScreenshotData)
end
return nothing
end

Expand Down
32 changes: 19 additions & 13 deletions test/test_ProfileProcessing.jl
Original file line number Diff line number Diff line change
Expand Up @@ -59,6 +59,7 @@ prof1 = ProfileData(start_lonlat = (5, 45), end_lonlat = (15, 49))
prof2 = ProfileData(depth = -100)
prof3 = ProfileData(start_lonlat = (5, 45), end_lonlat = (5, 49))
prof4 = ProfileData(depth = -20)
prof5 = ProfileData(start_lonlat = (12, 49), end_lonlat = (12, 45))

# test internal routines to intersect profile with volumetric data:
GeophysicalModelGenerator.create_profile_volume!(prof1, VolData_combined1)
Expand All @@ -80,27 +81,32 @@ GeophysicalModelGenerator.create_profile_point!(prof4, Data.Point, section_width
@test length(prof1.PointData[1].lon) == 13
@test length(prof4.PointData[1].lon) == 445

# test screenshot data
GeophysicalModelGenerator.create_profile_screenshot!(prof5, Data.Screenshot)
@test prof5.SurfData[1].fields.x_profile[1,1,1] == 0


# Test the main profile extraction routines:
extract_ProfileData!(prof1, VolData_combined1, Data.Surface, Data.Point)
extract_ProfileData!(prof2, VolData_combined1, Data.Surface, Data.Point)
extract_ProfileData!(prof3, VolData_combined1, Data.Surface, Data.Point)
extract_ProfileData!(prof4, VolData_combined1, Data.Surface, Data.Point)
extract_ProfileData!(prof1, VolData_combined1, Data.Surface, Data.Point, NamedTuple())
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

see remark above

extract_ProfileData!(prof2, VolData_combined1, Data.Surface, Data.Point, NamedTuple())
extract_ProfileData!(prof3, VolData_combined1, Data.Surface, Data.Point, NamedTuple())
extract_ProfileData!(prof4, VolData_combined1, Data.Surface, Data.Point, NamedTuple())

extract_ProfileData!(prof1, VolData_combined2, Data.Surface, Data.Point)
extract_ProfileData!(prof2, VolData_combined2, Data.Surface, Data.Point)
extract_ProfileData!(prof3, VolData_combined2, Data.Surface, Data.Point)
extract_ProfileData!(prof4, VolData_combined2, Data.Surface, Data.Point)
extract_ProfileData!(prof1, VolData_combined2, Data.Surface, Data.Point, NamedTuple())
extract_ProfileData!(prof2, VolData_combined2, Data.Surface, Data.Point, NamedTuple())
extract_ProfileData!(prof3, VolData_combined2, Data.Surface, Data.Point, NamedTuple())
extract_ProfileData!(prof4, VolData_combined2, Data.Surface, Data.Point, NamedTuple())

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!(prof1, VolData_combined3, Data.Surface, Data.Point, NamedTuple())
extract_ProfileData!(prof2, VolData_combined3, Data.Surface, Data.Point, NamedTuple())
extract_ProfileData!(prof3, VolData_combined3, Data.Surface, Data.Point, NamedTuple())
extract_ProfileData!(prof4, VolData_combined3, Data.Surface, Data.Point, NamedTuple())
extract_ProfileData!(prof5, VolData_combined3, Data.Surface, Data.Point, Data.Screenshot)


# Test that it works if only EQ's are provided:
prof4 = ProfileData(depth = -20)
extract_ProfileData!(prof4, nothing, NamedTuple(), Data.Point)
extract_ProfileData!(prof4, nothing, NamedTuple(), Data.Point, NamedTuple())
@test isnothing(prof4.VolData)
@test isempty(prof4.SurfData)
@test length(prof4.PointData[1].depth) == 3280
Expand Down
Loading