Merge pull request #134 from JuliaGeodynamics/pa-cross

fix GeoData export

Author: boriskaus

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
+end

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