Skip to content

Commit

Permalink
Merge pull request #134 from JuliaGeodynamics/pa-cross
Browse files Browse the repository at this point in the history 
fix GeoData export
  • Loading branch information
@boriskaus
boriskaus authored Jul 29, 2024
2 parents 23d9d4a + 4094505 commit fdd372a
Show file tree
Hide file tree
Showing 2 changed files with 47 additions and 28 deletions.
33 changes: 19 additions & 14 deletions src/utils.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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.
Expand Down Expand Up @@ -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.
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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]
Expand Down Expand Up @@ -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


Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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)

Expand Down Expand Up @@ -1815,4 +1820,4 @@ function inpoly_fast(PolyX::Vector{T}, PolyY::Vector{T}, x::T, y::T) where T <:
end
end
return inside
end
end
42 changes: 28 additions & 14 deletions test/test_utils.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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)
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down

0 comments on commit fdd372a

Please sign in to comment.