Skip to content

Commit 41c8e38

Browse files
authored
Merge pull request #1 from jfrasunk/jf-readtomo
new function to import tomographies from NetCDF with tests
2 parents a9ecd79 + c328123 commit 41c8e38

File tree

2 files changed

+42
-1
lines changed

2 files changed

+42
-1
lines changed

src/data_import.jl

Lines changed: 29 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -6,7 +6,7 @@
66

77
using LightXML
88

9-
export screenshot_to_GeoData, screenshot_to_CartData, screenshot_to_UTMData, getlonlatdepthmag_QuakeML
9+
export screenshot_to_GeoData, screenshot_to_CartData, screenshot_to_UTMData, getlonlatdepthmag_QuakeML, tomo_2_GeoData
1010

1111
# import CSV data using standard library functions
1212
# here we assume that the data is indeed comma separated and that comments are preceded with a "#"
@@ -350,3 +350,31 @@ function getlonlatdepthmag_QuakeML(filename::String)
350350
Data_ISC = GeoData(lon, lat, -1 * depth / 1.0e3, (Magnitude = mag, Depth = -1 * depth / 1.0e3 * km))
351351
return Data_ISC
352352
end
353+
354+
"""
355+
Read_TomoData(filename::String)
356+
357+
Reads a seismic tomography dataset from a NetCDF file as a GeoData object. The keyword argument `vel_type::String` allows you to specify the type of velocity data to extract (default is "vs" for shear wave velocity).
358+
The function assumes that the NetCDF file contains variables for depth, longitude, latitude, and the specified velocity type.
359+
360+
tomodata = tomo_2_GeoData("path/to/tomo_data.nc")
361+
"""
362+
function tomo_2_GeoData(filename::String; vel_type::String = "vs")
363+
364+
# Open the NetCDF file
365+
data = NCDataset(filename)
366+
367+
# Extract the variables
368+
depth = data["depth"][:]
369+
lon = data["longitude"][:]
370+
lat = data["latitude"][:]
371+
vel = data[vel_type][:,:,:]
372+
373+
# create lon, lat, depth grid
374+
Lon, Lat, Depth = lonlatdepth_grid(lon, lat, .- depth)
375+
376+
# greate GeoDate struct
377+
Tomo_data = GeoData(Lon, Lat, Depth, (vel=vel,))
378+
379+
return Tomo_data
380+
end

test/test_data_import.jl

Lines changed: 13 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -87,3 +87,16 @@ data_Image = screenshot_to_UTMData(filename, Corner_LowerLeft, Corner_UpperRight
8787
Data_ISC = getlonlatdepthmag_QuakeML("test_files/ISCTest.xml");
8888
@test Value(Data_ISC.depth[1]) == -13.0km
8989
@test Data_ISC.fields.Magnitude[1] == 5.8
90+
91+
# test the import of a seismic tomography model from a NetCDF file
92+
tomo_file = "test_files/test_tomo_data.nc"
93+
Tomo_data = tomo_2_GeoData(tomo_file)
94+
test_array = [i for i in 0:1:10]
95+
96+
@test Tomo_data.lon.val[:,1,1] test_array
97+
@test Tomo_data.lat.val[1,:,1] test_array
98+
@test Tomo_data.depth.val[1,1,:] .- test_array
99+
100+
for i in eachindex(test_array)
101+
@test Tomo_data.fields.vel[i,i,i] i
102+
end

0 commit comments

Comments
 (0)