Merge pull request #171 from jfrasunk/main

new function to import seismic tomographies from NetCDF files

Author: boriskaus

docs/src/man/dataimport.md

@@ -8,4 +8,5 @@ GeophysicalModelGenerator.screenshot_to_CartData
 GeophysicalModelGenerator.screenshot_to_UTMData
 GeophysicalModelGenerator.import_topo
 GeophysicalModelGenerator.import_GeoTIFF
+GeophysicalModelGenerator.tomo_2_GeoData
 ```

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)
+
+    # create GeoData struct
+    Tomo_data = GeoData(Lon, Lat, Depth, (vel=vel,))
+
+    return Tomo_data
+end

test/test_GMT.jl

@@ -1,11 +1,11 @@
 using Test
 using GeophysicalModelGenerator, GMT
 
-Topo = import_topo(lat = [30, 31], lon = [50, 51])
-@test sum(Topo.depth.val) ≈ 2777.5705
+Topo = import_topo(lon = [8, 9], lat = [50, 51])
+@test sum(Topo.depth.val) ≈ 1076.7045 rtol = 1e-1
 
-Topo = import_topo([50, 51, 30, 31]);
-@test sum(Topo.depth.val) ≈ 2777.5705
+Topo = import_topo([8, 9, 50, 51]);
+@test sum(Topo.depth.val) ≈ 1076.7045 rtol = 1e-1
 
 test_fwd = import_GeoTIFF("test_files/length_fwd.tif", fieldname = :forward)
 @test  maximum(test_fwd.fields.forward) ≈ 33.17775km

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