new function to import seismic tomographies from NetCDF files

src/data_import.jl

@@ -6,7 +6,7 @@
 
 using LightXML
 
-export screenshot_to_GeoData, screenshot_to_CartData, screenshot_to_UTMData, getlonlatdepthmag_QuakeML
+export screenshot_to_GeoData, screenshot_to_CartData, screenshot_to_UTMData, getlonlatdepthmag_QuakeML, tomo_2_GeoData
 
 # import CSV data using standard library functions
 # 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)
     Data_ISC = GeoData(lon, lat, -1 * depth / 1.0e3, (Magnitude = mag, Depth = -1 * depth / 1.0e3 * km))
     return Data_ISC
 end
+
+"""
+    Read_TomoData(filename::String)
+
+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).
+The function assumes that the NetCDF file contains variables for depth, longitude, latitude, and the specified velocity type.
+
+tomodata = tomo_2_GeoData("path/to/tomo_data.nc")
+"""
+function tomo_2_GeoData(filename::String; vel_type::String = "vs")
+
+    # Open the NetCDF file
+    data = NCDataset(filename)
+
+    # Extract the variables
+    depth = data["depth"][:]
+    lon   = data["longitude"][:]
+    lat   = data["latitude"][:]
+    vel   = data[vel_type][:,:,:]
+
+    # create lon, lat, depth grid
+    Lon, Lat, Depth = lonlatdepth_grid(lon, lat, .- depth)
+
+    # greate GeoDate struct
+    Tomo_data = GeoData(Lon, Lat, Depth, (vel=vel,))
+
+    return Tomo_data
+end

test/test_data_import.jl

@@ -87,3 +87,16 @@ data_Image = screenshot_to_UTMData(filename, Corner_LowerLeft, Corner_UpperRight
 Data_ISC = getlonlatdepthmag_QuakeML("test_files/ISCTest.xml");
 @test Value(Data_ISC.depth[1]) == -13.0km
 @test Data_ISC.fields.Magnitude[1] == 5.8
+
+# test the import of a seismic tomography model from a NetCDF file
+tomo_file = "test_files/test_tomo_data.nc"
+Tomo_data = tomo_2_GeoData(tomo_file)
+test_array = [i for i in 0:1:10]
+
+@test Tomo_data.lon.val[:,1,1] ≈ test_array
+@test Tomo_data.lat.val[1,:,1] ≈ test_array
+@test Tomo_data.depth.val[1,1,:] ≈ .- test_array
+
+for i in eachindex(test_array)
+    @test Tomo_data.fields.vel[i,i,i] ≈ i
+end