Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
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
53 changes: 49 additions & 4 deletions src/ProfileProcessing.jl
Original file line number Diff line number Diff line change
Expand Up @@ -17,6 +17,7 @@ Structure that holds profile data (interpolated/projected on the profile)
VolData :: GeophysicalModelGenerator.GeoData
SurfData :: Union{Nothing, NamedTuple}
PointData :: Union{Nothing, NamedTuple}
ScreenshotData :: Union{Nothing, NamedTuple}
end

Structure to store cross section data
Expand All @@ -29,6 +30,7 @@ mutable struct ProfileData
VolData::Union{Nothing, GeophysicalModelGenerator.GeoData}
SurfData::Union{Nothing, NamedTuple}
PointData::Union{Nothing, NamedTuple}
ScreenshotData::Union{Nothing, NamedTuple}

function ProfileData(; kwargs...) # this constructor allows to define only certain fields and leave the others blank
K = new(true, nothing, nothing, nothing, nothing, nothing, nothing)
Expand All @@ -51,7 +53,6 @@ mutable struct ProfileData
end
end


function show(io::IO, g::ProfileData)
if g.vertical
println(io, "Vertical ProfileData")
Expand All @@ -69,7 +70,9 @@ function show(io::IO, g::ProfileData)
if !isnothing(g.PointData)
println(io, " PointData: $(keys(g.PointData)) ")
end

if !isnothing(g.ScreenshotData)
println(io, " ScreenshotData: $(keys(g.ScreenshotData)) ")
end
return nothing
end

Expand Down Expand Up @@ -282,6 +285,31 @@ function create_profile_volume!(Profile::ProfileData, VolData::AbstractGeneralGr
return nothing
end

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

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

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

# add the data set as a NamedTuple
data_NT = NamedTuple{(DataSetName[idata],)}((data_tmp,))
tmp = merge(tmp, data_NT)
else
# we do not have this implemented
#error("horizontal profiles not yet implemented")
end
end
Profile.SurfData = tmp # assign to profile data structure
return
end

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


"""
extract_ProfileData!(Profile::ProfileData,VolData::GeoData, SurfData::NamedTuple, PointData::NamedTuple; DimsVolCross=(100,100),Depth_extent=nothing,DimsSurfCross=(100,),section_width=50)
extract_ProfileData!(Profile::ProfileData,VolData::GeoData, SurfData::NamedTuple, PointData::NamedTuple, ScreenshotData::NamedTuple; DimsVolCross=(100,100),Depth_extent=nothing,DimsSurfCross=(100,),section_width=50)

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

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

"""
extract_ProfileData!(Profile::ProfileData,VolData::GeoData, SurfData::NamedTuple, PointData::NamedTuple; DimsVolCross=(100,100),Depth_extent=nothing,DimsSurfCross=(100,),section_width=50)

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

# call the actual function to extract the profile data
extract_ProfileData!(Profile, VolData, SurfData, PointData,NamedTuple(); DimsVolCross = DimsVolCross, Depth_extent = Depth_extent, DimsSurfCross = DimsSurfCross, section_width = section_width)
return nothing
end




"""
This reads the picked profiles from disk and returns a vector of ProfileData
"""
Expand Down
7 changes: 6 additions & 1 deletion test/test_ProfileProcessing.jl
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,6 @@ data_Vol2 = GMG_Dataset("Plomerova2022", "Volume", "https://seafile.rlp.net/f/ab
#data_Vol2 = GMG_Dataset("Zhao2016","Volume","https://seafile.rlp.net/f/e81a6d075f6746609973/?dl=1", true)

# Now load these datasets into NamedTuples

SurfData = load_GMG(data_Surf)
PointData = load_GMG(data_EQ)
ScreenshotData = load_GMG(data_SS)
Expand Down Expand Up @@ -59,6 +58,7 @@ prof1 = ProfileData(start_lonlat = (5, 45), end_lonlat = (15, 49))
prof2 = ProfileData(depth = -100)
prof3 = ProfileData(start_lonlat = (5, 45), end_lonlat = (5, 49))
prof4 = ProfileData(depth = -20)
prof5 = ProfileData(start_lonlat = (12, 49), end_lonlat = (12, 45))

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

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


# Test the main profile extraction routines:
extract_ProfileData!(prof1, VolData_combined1, Data.Surface, Data.Point)
Expand All @@ -96,6 +100,7 @@ extract_ProfileData!(prof1, VolData_combined3, Data.Surface, Data.Point)
extract_ProfileData!(prof2, VolData_combined3, Data.Surface, Data.Point)
extract_ProfileData!(prof3, VolData_combined3, Data.Surface, Data.Point)
extract_ProfileData!(prof4, VolData_combined3, Data.Surface, Data.Point)
extract_ProfileData!(prof5, VolData_combined3, Data.Surface, Data.Point, Data.Screenshot)


# Test that it works if only EQ's are provided:
Expand Down
Loading