Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

fix GeoData export #134

Merged
merged 2 commits into from
Jul 29, 2024
Merged

fix GeoData export #134

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
e4f1977
fix GeoData export
aelligp Jul 24, 2024
4094505
add test
aelligp Jul 29, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
33 changes: 19 additions & 14 deletions src/utils.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -70,7 +70,7 @@

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 @@



"""
"""
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 @@
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 @@

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")

Check warning on line 360 in src/utils.jl

View check run for this annotation

Codecov / codecov/patch

src/utils.jl#L360 

Added line #L360 was not covered by tests
end
return Data_profile
end

Expand Down Expand Up @@ -545,9 +550,9 @@
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 @@
```
"""
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 @@
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 @@

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 @@
end
end
return inside
end
end
42 changes: 28 additions & 14 deletions test/test_utils.jl
View file
Open in desktop
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
Loading