From 4094505111852ec56f9a4cae8efb2d59abfe30a4 Mon Sep 17 00:00:00 2001 From: Pascal Aellig Date: Mon, 29 Jul 2024 12:03:25 +0200 Subject: [PATCH] add test --- test/test_utils.jl | 42 ++++++++++++++++++++++++++++-------------- 1 file changed, 28 insertions(+), 14 deletions(-) diff --git a/test/test_utils.jl b/test/test_utils.jl index b5afdd31..519c0a7e 100644 --- a/test/test_utils.jl +++ b/test/test_utils.jl @@ -7,8 +7,8 @@ using GeophysicalModelGenerator Lon,Lat,Depth = lonlatdepth_grid(10:20,30:40,-50km); Data1 = Depth*2; # some data Vx1,Vy1,Vz1 = Data1*3,Data1*4,Data1*5 -Data_set2D = GeoData(Lon,Lat,Depth,(Depthdata=Data1,LonData1=Lon, Velocity=(Vx1,Vy1,Vz1))) -Data_set2D0 = GeoData(Lon,Lat,Depth,(Depthdata=Data1,LonData1=Lon)) +Data_set2D = GeoData(Lon,Lat,Depth,(Depthdata=Data1,LonData1=Lon, Velocity=(Vx1,Vy1,Vz1))) +Data_set2D0 = GeoData(Lon,Lat,Depth,(Depthdata=Data1,LonData1=Lon)) @test_throws ErrorException cross_section(Data_set2D, Depth_level=-10) # Test interpolation of depth to a given cartesian XY-plane @@ -22,7 +22,7 @@ plane2 = interpolate_datafields_2D(Data_set2D, proj, x, y) Lon1,Lat1,Depth1 = lonlatdepth_grid(12:18,33:39,-50km); Data2 = Depth1*2; # some data Vx1,Vy1,Vz1 = Data2*3,Data2*4,Data2*5 -Data_set2D_1 = GeoData(Lon1,Lat1,Depth1,(Depthdata1=Data2,LonData2=Lon1)) +Data_set2D_1 = GeoData(Lon1,Lat1,Depth1,(Depthdata1=Data2,LonData2=Lon1)) plane3 = interpolate_datafields_2D(Data_set2D0, Data_set2D_1) @test sum(plane3.fields.Depthdata) ≈ -4900.0km @@ -35,23 +35,23 @@ plane3 = interpolate_datafields_2D(Data_set2D0, Data_set2D_1) Lon,Lat,Depth = lonlatdepth_grid(10:20,30:40,(-300:25:0)km); Data = Depth*2; # some data Vx,Vy,Vz = ustrip(Data*3)*km/s,ustrip(Data*4)*km/s,ustrip(Data*5)*km/s; -Data_set3D = GeoData(Lon,Lat,Depth,(Depthdata=Data,LonData=Lon, Velocity=(Vx,Vy,Vz))) +Data_set3D = GeoData(Lon,Lat,Depth,(Depthdata=Data,LonData=Lon, Velocity=(Vx,Vy,Vz))) # Test addfield -Data_set3D = addfield(Data_set3D,"Lat", Lat) +Data_set3D = addfield(Data_set3D,"Lat", Lat) @test keys(Data_set3D.fields) == (:Depthdata, :LonData, :Velocity, :Lat) -Data_set3D = addfield(Data_set3D,(;Lat, Lon)) +Data_set3D = addfield(Data_set3D,(;Lat, Lon)) @test keys(Data_set3D.fields) == (:Depthdata, :LonData, :Velocity, :Lat, :Lon) # Create 3D cartesian dataset -Data_setCart3D = CartData(Lon,Lat,Depth,(Depthdata=Data,LonData=Lon, Velocity=(Vx,Vy,Vz))) +Data_setCart3D = CartData(Lon,Lat,Depth,(Depthdata=Data,LonData=Lon, Velocity=(Vx,Vy,Vz))) # Create 3D volume with some fake data Lon,Lat,Depth = lonlatdepth_grid(10:20,30:40,(0:-25:-300)km); Data = Depth*2; # some data Vx,Vy,Vz = ustrip(Data*3)*km/s,ustrip(Data*4)*km/s,ustrip(Data*5)*km/s; -Data_set3D_reverse = GeoData(Lon,Lat,Depth,(Depthdata=Data,LonData=Lon, Velocity=(Vx,Vy,Vz))) +Data_set3D_reverse = GeoData(Lon,Lat,Depth,(Depthdata=Data,LonData=Lon, Velocity=(Vx,Vy,Vz))) # Create cross-sections in various directions (no interpolation which is default) test_cross = cross_section(Data_set3D, Depth_level=-100km) @@ -99,6 +99,20 @@ test_cross = cross_section(Data_set3D, Start=(10,30), End=(20,40), dims=( #@test size(test_cross_rev.fields[3][2])==(50,100,1) #@test write_paraview(test_cross_rev, "profile_test_rev")[1]=="profile_test_rev.vts" +# Cross section of a topography +depth_values = [rand(0:0.1:3.5)] +Lon, Lat, Depth =lonlatdepth_grid(10:20, 30:40, depth_values[:]); +Data_Topo = GeoData(Lon, Lat, Depth, (Depthdata=Depth,)) +Data_Topo_geo= cross_section(Data_Topo, Start=(10,30), End=(20,40), dims=(50,100), Interpolate=true) +@test Data_Topo_geo isa GeoData + +Lon,Lat,Depth = lonlatdepth_grid(5:25,20:50,0); +Depth = cos.(Lon/5).*sin.(Lat)*10; +Data_surf = GeoData(Lon,Lat,Depth,(Z=Depth,)); +Data_surf_cart = convert2CartData(Data_surf, proj); +Data_surf_cross = cross_section(Data_surf_cart, Start=(-1693,2500), End=(-1000,3650), dims=(50,100), Interpolate=true) +@test Data_surf_cross isa CartData + # Cross-section with cartesian data test_cross = cross_section(Data_setCart3D, Lon_level=15, dims=(50,100), Interpolate=true) @test size(test_cross.fields[3][2])==(1,50,100) @@ -173,16 +187,16 @@ Data_pert = subtract_horizontalmean(Data, Percentage=true) # 3D w @test Data_pert[1000] == 0.0 Data2D = Data[:,1,:]; -Data_pert = subtract_horizontalmean(Data2D, Percentage=true) # 2D version with units [dp the same along a vertical profile] +Data_pert = subtract_horizontalmean(Data2D, Percentage=true) # 2D version with units [dp the same along a vertical profile] -Data_set2D = GeoData(Lon,Lat,Depth,(Depthdata=Data,LonData=Lon,Pertdata=Data_pert ,Velocity=(Vx,Vy,Vz))) +Data_set2D = GeoData(Lon,Lat,Depth,(Depthdata=Data,LonData=Lon,Pertdata=Data_pert ,Velocity=(Vx,Vy,Vz))) @test Data_set2D.fields[3][10,8,1] == 0 # Create surface ("Moho") Lon,Lat,Depth = lonlatdepth_grid(10:20,30:40,-40km); Depth = Depth + Lon*km; # some fake topography on Moho -Data_Moho = GeoData(Lon,Lat,Depth,(MohoDepth=Depth,LonData=Lon,TestData=(Depth,Depth,Depth))) +Data_Moho = GeoData(Lon,Lat,Depth,(MohoDepth=Depth,LonData=Lon,TestData=(Depth,Depth,Depth))) # Test intersecting a surface with 2D or 3D data sets @@ -218,7 +232,7 @@ Data_VoteMap = votemap(Data_set3D_reverse, "Depthdata<-560",dims=(10,10,10)) @test Data_VoteMap.fields[:votemap][101]==0 @test Data_VoteMap.fields[:votemap][100]==1 -# Combine 2 datasets +# Combine 2 datasets Data_VoteMap = votemap([Data_set3D_reverse, Data_set3D], ["Depthdata<-560","LonData>19"],dims=(10,10,10)) @test Data_VoteMap.fields[:votemap][10,9,1]==2 @test Data_VoteMap.fields[:votemap][9 ,9,1]==1 @@ -253,9 +267,9 @@ cross_tmp = cross_section(Data_EQ,Lat_level=36.2,section_width=10km) @test cross_tmp.fields.lat_proj[10]==36.2 # check if the projected latitude level is the chosen one cross_tmp = cross_section(Data_EQ,Lon_level=16.4,section_width=10km) -@test cross_tmp.fields.lon_proj[10]==16.4 # check if the projected longitude level is the chosen one +@test cross_tmp.fields.lon_proj[10]==16.4 # check if the projected longitude level is the chosen one cross_tmp = cross_section(Data_EQ,Start=(15.0,35.0),End=(17.0,37.0),section_width=10km) -@test cross_tmp.fields.lon_proj[20] ==15.314329874961091 +@test cross_tmp.fields.lon_proj[20] ==15.314329874961091 @test cross_tmp.fields.lat_proj[20] == 35.323420618580585 # test inPolygon