diff --git a/src/utils.jl b/src/utils.jl index fc276bf6..7609bbec 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -70,7 +70,7 @@ addfield(V::GeoData,new_fields::NamedTuple) = GeoData(V.lon.val, V.lat.val, V.de Add `new_fields` fields to a `Q1Data` dataset; set `cellfield` to `true` if the field is a cell field; otherwise it is a vertex field """ -function addfield(V::Q1Data,new_fields::NamedTuple; cellfield=false) +function addfield(V::Q1Data,new_fields::NamedTuple; cellfield=false) if cellfield return Q1Data(V.x.val, V.y.val, V.z.val, V.fields, merge(V.cellfields, new_fields)) else @@ -122,7 +122,7 @@ end -""" +""" cross_section_volume(Volume::AbstractGeneralGrid; dims=(100,100), Interpolate=false, Depth_level=nothing; Lat_level=nothing; Lon_level=nothing; Start=nothing, End=nothing, Depth_extent=nothing ) Creates a cross-section through a volumetric (3D) `GeoData` object. @@ -259,7 +259,7 @@ cross_section_surface(Surface::GeoData; dims=(100,), Interpolate=false, Depth_le Creates a cross-section through a surface (2D) `GeoData` object. - Cross-sections can be horizontal (map view at a given depth), if `Depth_level` is specified -- They can also be vertical, either by specifying `Lon_level` or `Lat_level` (for a fixed lon/lat), or by defining both `Start=(lon,lat)` & `End=(lon,lat)` points. +- They can also be vertical, either by specifying `Lon_level` or `Lat_level` (for a fixed lon/lat), or by defining both `Start=(lon,lat)` & `End=(lon,lat)` points. Start and End points will be in km! - IMPORTANT: The surface to be extracted has to be given as a gridded GeoData object. It may also contain NaNs where it is not defined. Any points lying outside of the defined surface will be considered NaN. @@ -351,9 +351,14 @@ function cross_section_surface(S::AbstractGeneralGrid; dims=(100,), Interpolate= end - # create GeoData structure with the interpolated points - Data_profile = GeoData(Lon, Lat, depth_intp, (fields_new)); - + # create GeoData/CartData structure with the interpolated points + if isa(S,GeoData) + Data_profile = GeoData(Lon, Lat, depth_intp, fields_new); + elseif isa(S,CartData) + Data_profile = CartData(Lon, Lat, depth_intp, fields_new); + else + error("still to be implemented") + end return Data_profile end @@ -545,9 +550,9 @@ Creates a cross-section through a `GeoData` object. julia> Lon,Lat,Depth = lonlatdepth_grid(10:20,30:40,(-300:25:0)km); julia> Data = Depth*2; # some data julia> Vx,Vy,Vz = ustrip(Data*3),ustrip(Data*4),ustrip(Data*5); -julia> Data_set3D = GeoData(Lon,Lat,Depth,(Depthdata=Data,LonData=Lon, Velocity=(Vx,Vy,Vz))); -julia> Data_cross = cross_section(Data_set3D, Depth_level=-100km) -GeoData +julia> Data_set3D = GeoData(Lon,Lat,Depth,(Depthdata=Data,LonData=Lon, Velocity=(Vx,Vy,Vz))); +julia> Data_cross = cross_section(Data_set3D, Depth_level=-100km) +GeoData size : (11, 11, 1) lon ϵ [ 10.0 : 20.0] lat ϵ [ 30.0 : 40.0] @@ -601,7 +606,7 @@ CartData ``` """ function flatten_cross_section(V::CartData) - + x_new = sqrt.((V.x.val.-V.x.val[1,1,1]).^2 .+ (V.y.val.-V.y.val[1,1,1]).^2) # NOTE: the result is in km, as V.x and V.y are stored in km @@ -1244,11 +1249,11 @@ function interpolate_datafields_2D(Original::GeoData, New::GeoData; Rotate=0.0, return GeoData(New.lon.val,New.lat.val,Znew, fields_new) end -""" +""" Surf_interp = interpolate_datafields_2D(V::GeoData, x::AbstractRange, y::AbstractRange; Lat::Number, Lon::Number) Interpolates a 3D data set `V` with a projection point `proj=(Lat, Lon)` on a plane defined by `x` and `y`, where `x` and `y` are uniformly spaced. -Returns the 2D array `Surf_interp`. +Returns the 2D array `Surf_interp`. """ function interpolate_datafields_2D(V::GeoData, x::AbstractRange, y::AbstractRange; Lat=49.9929, Lon=8.2473) # Default: Lat=49.9929, Lon=8.2473 => Mainz (center of universe) @@ -1340,7 +1345,7 @@ function InterpolateDataFields2D_vecs(EW_vec, NS_vec, depth, fields_new, EW, NS) return depth_new, fields_new end - + # Extracts a sub-data set using indices function ExtractDataSets(V::AbstractGeneralGrid, iLon, iLat, iDepth) @@ -1815,4 +1820,4 @@ function inpoly_fast(PolyX::Vector{T}, PolyY::Vector{T}, x::T, y::T) where T <: end end return inside -end \ No newline at end of file +end 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