Skip to content

Additional Alpine datasets and updated GMG version #5

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 9 commits into
base: main
Choose a base branch
from
+87 −52
Open

Additional Alpine datasets and updated GMG version #5

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
95d4c2c
update version numbers
boriskaus Mar 18, 2025
c401628
add Moho surfaces
boriskaus Mar 18, 2025
22210f5
update routines to new GMG version
boriskaus Mar 18, 2025
2753690
update functions due to update GMG version
boriskaus Mar 18, 2025
b00935b
add Vp,Vs, seismic and topo data
boriskaus Mar 18, 2025
cf9044d
add SAVANI_Medit_Vs
boriskaus Mar 19, 2025
378c410
added SAVANI_Alps_Vs
boriskaus Mar 19, 2025
af64fb1
synthetic dataset
boriskaus May 19, 2025
770092d
update
boriskaus May 19, 2025
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
6 changes: 3 additions & 3 deletions Project.toml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -26,11 +26,11 @@ UrlDownload = "856ac37a-3032-4c1c-9122-f86d88358c8b"
Dash = "1.3.0"
DashBootstrapComponents = "1.4.1"
DelimitedFiles = "1"
GeophysicalModelGenerator = "0.5.4"
GeophysicalModelGenerator = "0.6, 0.7"
HTTP = "1.0"
JLD2 = "0.4"
JLD2 = "0.4, 0.5"
JSON3 = "1"
Meshes = "0.30"
Meshes = "0.30 - 0.40"
PlotlyJS = "0.18"
Statistics = "1.9"
UrlDownload = "1.0.1"
Expand Down
6 changes: 3 additions & 3 deletions playground/Dash/utils.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -57,12 +57,12 @@ Extracts a cross-section from a tomographic dataset and returns this as cartesia
function get_cross_section(DataAlps::GeoData, start_value=(10,41), end_value=(10,49), field=:dVp_paf21)

# retrieve the cross-section in GeoData format
cross = CrossSection(DataAlps, Start=start_value, End=end_value, Interpolate=true)
cross = cross_section(DataAlps, Start=start_value, End=end_value, Interpolate=true)

# transfer it to cartesian data
p = ProjectionPoint(Lon=minimum(cross.lon.val),Lat=minimum(cross.lat.val));
cross_cart = Convert2CartData(cross,p)
x_cross = FlattenCrossSection(cross_cart);
cross_cart = convert2CartData(cross,p)
x_cross = flatten_cross_section(cross_cart);
x_cart = x_cross[:,1];
z_cart = cross_cart.z.val[1,:,1]

Expand Down
16 changes: 8 additions & 8 deletions playground/Marcel/JuliaGTKMakie/ProfileProcessing.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@
# this is ProfileProcessing.jl
# It contains functions and type definitions to gather selected data for given profiles

export ProfileData, ExtractProfileData
export ProfileData, extract_ProfileData

# load packages
using GeophysicalModelGenerator
Expand Down Expand Up @@ -63,7 +63,7 @@ function CreateProfileVolume!(Profile,DataSetName,DataSetFile,DimsVolCross,Depth
tmp_load = load(DataSetFile[idata]) # this gives us a dict with a key that is the name if the data set and the values as the GeoData structure
tmp_load = collect(values(tmp_load)) # this gives us a vector with a single GeoData entry
data_tmp = tmp_load[1] # and now we extract that entry...
cross_tmp = CrossSection(data_tmp,dims=DimsVolCross,Start=Profile.start_point,End=Profile.end_point,Depth_extent=DepthVol) # create the cross section
cross_tmp = cross_section(data_tmp,dims=DimsVolCross,Start=Profile.start_point,End=Profile.end_point,Depth_extent=DepthVol) # create the cross section

# store profile coordinates and field data on first go
if idata==1
Expand Down Expand Up @@ -103,7 +103,7 @@ function CreateProfileVolume!(Profile,DataSetName,DataSetFile,DimsVolCross,Depth

tmp = GeophysicalModelGenerator.GeoData(lon_vol,lat_vol,depth_vol,fields_vol)
# flatten cross section and add this data to the structure
x_profile = FlattenCrossSection(tmp,Start=Profile.start_point)
x_profile = flatten_cross_section(tmp,Start=Profile.start_point)
tmp = AddField(tmp,"x_profile",x_profile)

Profile.VolData = tmp # assign to Profile data structure
Expand All @@ -123,9 +123,9 @@ function CreateProfileSurface!(Profile,DataSetName,DataSetFile,DimsSurfCross)
tmp_load = load(DataSetFile[idata]) # this gives us a dict with a key that is the name if the data set and the values as the GeoData structure
tmp_load = collect(values(tmp_load)) # this gives us a vector with a single GeoData entry
data_tmp = tmp_load[1] # and now we extract that entry...
tmp[idata] = CrossSection(data_tmp, dims=DimsSurfCross,Start=Profile.start_point,End=Profile.end_point) # create the cross section
tmp[idata] = cross_section(data_tmp, dims=DimsSurfCross,Start=Profile.start_point,End=Profile.end_point) # create the cross section
# flatten cross section and add this data to the structure
x_profile = FlattenCrossSection(tmp[idata],Start=Profile.start_point)
x_profile = flatten_cross_section(tmp[idata],Start=Profile.start_point)
tmp[idata] = AddField(tmp[idata],"x_profile",x_profile)
# add the data set name as an attribute (not required if there is proper metadata, but odds are that there is not)
tmp[idata].atts["dataset"] = DataSetName[idata]
Expand All @@ -147,9 +147,9 @@ function CreateProfilePoint!(Profile,DataSetName,DataSetFile,WidthPointProfile)
tmp_load = load(DataSetFile[idata]) # this gives us a dict with a key that is the name if the data set and the values as the GeoData structure
tmp_load = collect(values(tmp_load)) # this gives us a vector with a single GeoData entry
data_tmp = tmp_load[1] # and now we extract that entry...
tmp[idata] = CrossSection(data_tmp,Start=Profile.start_point,End=Profile.end_point,section_width = WidthPointProfile) # create the cross section
tmp[idata] = cross_section(data_tmp,Start=Profile.start_point,End=Profile.end_point,section_width = WidthPointProfile) # create the cross section
# flatten cross section and add this data to the structure
x_profile = FlattenCrossSection(tmp[idata],Start=Profile.start_point)
x_profile = flatten_cross_section(tmp[idata],Start=Profile.start_point)
tmp[idata] = AddField(tmp[idata],"x_profile",x_profile)
# add the data set name as an attribute (not required if there is proper metadata, but odds are that there is not)
tmp[idata].atts["dataset"] = DataSetName[idata]
Expand All @@ -161,7 +161,7 @@ function CreateProfilePoint!(Profile,DataSetName,DataSetFile,WidthPointProfile)
end

### wrapper function to process everything
function ExtractProfileData(ProfileCoordFile,ProfileNumber,DataSetName,DataSetFile,DataSetType,DimsVolCross,DepthVol,DimsSurfCross,WidthPointProfile)
function extract_ProfileData(ProfileCoordFile,ProfileNumber,DataSetName,DataSetFile,DataSetType,DimsVolCross,DepthVol,DimsSurfCross,WidthPointProfile)

# start and end points are saved in a text file
profile_data = readdlm(ProfileCoordFile,skipstart=1,',')
Expand Down
2 changes: 1 addition & 1 deletion playground/Marcel/PickingWorkflow.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -67,7 +67,7 @@ WidthPointProfile = 20km

# 2. process the profiles
include("ProfileProcessing.jl")
ExtractedData = ExtractProfileData(ProfileCoordFile,ProfileNumber,DataSetName,DataSetFile,DataSetType,DimsVolCross,DimsSurfCross,WidthPointProfile)
ExtractedData = extract_ProfileData(ProfileCoordFile,ProfileNumber,DataSetName,DataSetFile,DataSetType,DimsVolCross,DimsSurfCross,WidthPointProfile)

# 3. save data as MATLAB
fn = "Profile"*string(ProfileNumber)
Expand Down
16 changes: 8 additions & 8 deletions playground/Marcel/ProfileProcessing.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@
# this is ProfileProcessing.jl
# It contains functions and type definitions to gather selected data for given profiles

export ProfileData, ExtractProfileData
export ProfileData, extract_ProfileData

# load packages
using GeophysicalModelGenerator
Expand Down Expand Up @@ -60,7 +60,7 @@ function CreateProfileVolume!(Profile,DataSetName,DataSetFile,DimsVolCross)
tmp_load = load(DataSetFile[idata]) # this gives us a dict with a key that is the name if the data set and the values as the GeoData structure
tmp_load = collect(values(tmp_load)) # this gives us a vector with a single GeoData entry
data_tmp = tmp_load[1] # and now we extract that entry...
cross_tmp = CrossSection(data_tmp,dims=DimsVolCross,Start=Profile.start_point,End=Profile.end_point) # create the cross section
cross_tmp = cross_section(data_tmp,dims=DimsVolCross,Start=Profile.start_point,End=Profile.end_point) # create the cross section

# store profile coordinates and field data on first go
if idata==1
Expand Down Expand Up @@ -100,7 +100,7 @@ function CreateProfileVolume!(Profile,DataSetName,DataSetFile,DimsVolCross)

tmp = GeoData(lon_vol,lat_vol,depth_vol,fields_vol)
# flatten cross section and add this data to the structure
x_profile = FlattenCrossSection(tmp,Start=Profile.start_point)
x_profile = flatten_cross_section(tmp,Start=Profile.start_point)
tmp = AddField(tmp,"x_profile",x_profile)

Profile.VolData = tmp # assign to Profile data structure
Expand All @@ -117,9 +117,9 @@ function CreateProfileSurface!(Profile,DataSetName,DataSetFile,DimsSurfCross)
tmp_load = load(DataSetFile[idata]) # this gives us a dict with a key that is the name if the data set and the values as the GeoData structure
tmp_load = collect(values(tmp_load)) # this gives us a vector with a single GeoData entry
data_tmp = tmp_load[1] # and now we extract that entry...
tmp[idata] = CrossSection(data_tmp, dims=DimsSurfCross,Start=Profile.start_point,End=Profile.end_point) # create the cross section
tmp[idata] = cross_section(data_tmp, dims=DimsSurfCross,Start=Profile.start_point,End=Profile.end_point) # create the cross section
# flatten cross section and add this data to the structure
x_profile = FlattenCrossSection(tmp[idata],Start=Profile.start_point)
x_profile = flatten_cross_section(tmp[idata],Start=Profile.start_point)
tmp[idata] = AddField(tmp[idata],"x_profile",x_profile)
# add the data set name as an attribute (not required if there is proper metadata, but odds are that there is not)
tmp[idata].atts["dataset"] = DataSetName[idata]
Expand All @@ -141,9 +141,9 @@ function CreateProfilePoint!(Profile,DataSetName,DataSetFile,WidthPointProfile)
tmp_load = load(DataSetFile[idata]) # this gives us a dict with a key that is the name if the data set and the values as the GeoData structure
tmp_load = collect(values(tmp_load)) # this gives us a vector with a single GeoData entry
data_tmp = tmp_load[1] # and now we extract that entry...
tmp[idata] = CrossSection(data_tmp,Start=Profile.start_point,End=Profile.end_point,section_width = WidthPointProfile) # create the cross section
tmp[idata] = cross_section(data_tmp,Start=Profile.start_point,End=Profile.end_point,section_width = WidthPointProfile) # create the cross section
# flatten cross section and add this data to the structure
x_profile = FlattenCrossSection(tmp[idata],Start=Profile.start_point)
x_profile = flatten_cross_section(tmp[idata],Start=Profile.start_point)
tmp[idata] = AddField(tmp[idata],"x_profile",x_profile)
# add the data set name as an attribute (not required if there is proper metadata, but odds are that there is not)
tmp[idata].atts["dataset"] = DataSetName[idata]
Expand All @@ -155,7 +155,7 @@ function CreateProfilePoint!(Profile,DataSetName,DataSetFile,WidthPointProfile)
end

### wrapper function to process everything
function ExtractProfileData(ProfileCoordFile,ProfileNumber,DataSetName,DataSetFile,DataSetType,DimsVolCross,DimsSurfCross,WidthPointProfile)
function extract_ProfileData(ProfileCoordFile,ProfileNumber,DataSetName,DataSetFile,DataSetType,DimsVolCross,DimsSurfCross,WidthPointProfile)

# start and end points are saved in a text file
profile_data = readdlm(ProfileCoordFile,skipstart=1,',')
Expand Down
14 changes: 7 additions & 7 deletions playground/ex3_CrossSection.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,10 +6,10 @@ using GLMakie, ColorSchemes, JLD2, GeophysicalModelGenerator, Interpolations, Li
# a toggle selectes either vertical or horizontal crossection
# the postion of the vertical crossection can be selected in the Topograpgy plot in the upper right corner
# --> press "v" and then left mouseclick to select the startpoint of the cross section and "b" + left mouseclick for the End point (first press the key, then to the mouseclick while the key is still pressed!)
# the depth of the horizontal crosssection is changed by pressing "h" and "up arrow" key or "down arrow" key,respectively
# the depth of the horizontal cross_section is changed by pressing "h" and "up arrow" key or "down arrow" key,respectively


# this functions ouptuts all the information needed to plot a crosssection (LonLat values of Crossection, interpolated data, Topography)
# this functions ouptuts all the information needed to plot a cross_section (LonLat values of Crossection, interpolated data, Topography)
function PrepCrossSection(lonvec, latvec, depths,liftlons,liftlone,liftlats,liftlate,data,lonT,latT,dataTopo)

lons = lonvec[1]
Expand All @@ -19,15 +19,15 @@ function PrepCrossSection(lonvec, latvec, depths,liftlons,liftlone,liftlats,lift
zs = depths[1]
ze = depths[end]

# create crosssection
# create cross_section
n = 100
X,Y,Z = XYZGrid(LinRange(lons, lone, n), LinRange(lats, late, n), LinRange(ze, zs, n));
X,Y,Z = xyz_grid(LinRange(lons, lone, n), LinRange(lats, late, n), LinRange(ze, zs, n));
cross = GeoData(X, Y, Z, (ran = zeros(size(Z)),));

STup = @lift ($liftlons,$liftlats)
ETup = @lift ($liftlone,$liftlate)

surf = @lift CrossSection(cross, Start = $STup, End = $ETup,Interpolate=true)
surf = @lift cross_section(cross, Start = $STup, End = $ETup,Interpolate=true)
t1 = @lift NumValue($surf.lon); t1 = @lift $t1[:,:,1];
t2 = @lift NumValue($surf.lat); t2 = @lift $t2[:,:,1];
t3 = @lift NumValue($surf.depth); t3 = @lift $t3[:,:,1];
Expand Down Expand Up @@ -169,12 +169,12 @@ function make_plot()
length_sli = @lift length($depthr)

#################################
#### horizontal crosssection ####
#### horizontal cross_section ####
#################################

# toggle
toggle1 = Toggle(fig, active = false)
lab1 = Label(fig, lift(xtog -> xtog ? "Horizontal Crosssection visible" : "Vertical Crosssection invisible", toggle1.active))
lab1 = Label(fig, lift(xtog -> xtog ? "Horizontal cross_section visible" : "Vertical cross_section invisible", toggle1.active))

# update location of point
Hup = (Keyboard.h, Keyboard.up)
Expand Down
Loading