create surfaces function; addition and subtraction included; rename D…

…rapeOnTopo -> drape_on_topo; add tests

docs/src/man/Tutorial_Jura.md

@@ -47,7 +47,7 @@ Geology  = Screenshot_To_GeoData("SchoriM_Encl_01_Jura-map_A1.png", lowerleft, u
 You can "drape" this image on the topographic map with
 
 ```julia
-TopoGeology = DrapeOnTopo(Topo, Geology)
+TopoGeology = drape_on_topo(Topo, Geology)
 ```
 
 ```julia

docs/src/man/tools.md

@@ -18,7 +18,7 @@ GeophysicalModelGenerator.PointData2NearestGrid
 GeophysicalModelGenerator.Convert2UTMzone
 GeophysicalModelGenerator.Convert2CartData
 GeophysicalModelGenerator.ProjectCartData
-GeophysicalModelGenerator.DrapeOnTopo
+GeophysicalModelGenerator.drape_on_topo
 GeophysicalModelGenerator.LithostaticPressure!
 GeophysicalModelGenerator.CountMap
 ```

src/GeophysicalModelGenerator.jl

@@ -43,6 +43,7 @@ include("stl.jl")
 include("ProfileProcessing.jl")
 include("IO.jl")
 include("event_counts.jl")
+include("surface_functions.jl")
 
 # Add optional routines (only activated when the packages are loaded)
 

src/data_types.jl

@@ -1,7 +1,7 @@
 # This is data_types.jl
 # contains type definitions to be used in GeophysicalModelGenerator
 
-import Base: show
+import Base: show, size
 
 export  GeoData, ParaviewData, UTMData, CartData,
         LonLatDepthGrid, XYZGrid, Velocity_SphericalToCartesian!,
@@ -216,6 +216,7 @@ struct GeoData <: AbstractGeneralGrid
      end
 
 end
+size(d::GeoData) = size(d.lon.val)
 
 # Print an overview of the Geodata struct:
 function Base.show(io::IO, d::GeoData)
@@ -247,6 +248,23 @@ function Base.show(io::IO, d::GeoData)
 end
 
 
+"""
+    GeoData(lld::Tuple{Array,Array,Array})
+
+This creates a `GeoData` struct if you have a Tuple with 3D coordinates as input.
+# Example
+```julia
+julia> data = GeoData(LonLatDepthGrid(-10:10,-5:5,0))
+GeoData 
+  size      : (21, 11, 1)
+  lon       ϵ [ -10.0 : 10.0]
+  lat       ϵ [ -5.0 : 5.0]
+  depth     ϵ [ 0.0 : 0.0]
+  fields    : (:Z,)
+```
+"""
+GeoData(lld::Tuple) = GeoData(lld[1],lld[2],lld[3],(Z=lld[3],))
+
 
 """
     ParaviewData(x::GeoUnit, y::GeoUnit, z::GeoUnit, values::NamedTuple)
@@ -265,6 +283,7 @@ mutable struct ParaviewData <: AbstractGeneralGrid
     z       ::  GeoUnit
     fields  ::  NamedTuple
 end
+size(d::ParaviewData) = size(d.x.val)
 
 # Print an overview of the ParaviewData struct:
 function Base.show(io::IO, d::ParaviewData)
@@ -464,6 +483,7 @@ struct UTMData <: AbstractGeneralGrid
      end
 
 end
+size(d::UTMData) = size(d.EW.val)
 
 # Print an overview of the UTMData struct:
 function Base.show(io::IO, d::UTMData)
@@ -748,6 +768,7 @@ struct CartData <: AbstractGeneralGrid
      end
 
 end
+size(d::CartData) = size(d.x.val)
 
 # Print an overview of the UTMData struct:
 function Base.show(io::IO, d::CartData)
@@ -1032,6 +1053,7 @@ struct CartGrid{FT, D} <: AbstractGeneralGrid
     coord1D     :: NTuple{D,StepRangeLen{FT}}   # Tuple with 1D vectors in all directions
     coord1D_cen :: NTuple{D,StepRangeLen{FT}}   # Tuple with 1D vectors of center points in all directions
 end
+size(d::CartGrid) = d.N
 
 
 """
@@ -1149,8 +1171,6 @@ function CreateCartGrid(;
 
 end
 
-
-
 # view grid object
 function show(io::IO, g::CartGrid{FT, DIM}) where {FT, DIM}
 

src/surface_functions.jl

@@ -0,0 +1,160 @@
+# This contains a number of routines that deal with surfaces
+using NearestNeighbors
+
+export remove_NaN_Surface!, drape_on_topo, is_surface
+import Base: +,-
+
+"""
+    issurf = is_surface(surf::AbstractGeneralGrid)
+
+Returns true if `surf` is a horizontal 3D surface.
+"""
+function is_surface(surf::AbstractGeneralGrid)
+    if size(surf)[3] == 1
+        issurf = true
+    else
+        issurf = false
+    end
+    return issurf
+end
+
+
+function  +(a::_T, b::_T) where _T<:AbstractGeneralGrid 
+    @assert size(a) == size(b)
+    return _addSurfaces(a,b)
+end
+
+function  -(a::_T, b::_T) where _T<:AbstractGeneralGrid 
+    @assert size(a) == size(b)
+    return _subtractSurfaces(a,b)
+end
+
+# Internal routines
+_addSurfaces(a::_T, b::_T) where _T<:GeoData        = GeoData(a.lon.val, a.lat.val, a.depth.val + b.depth.val, merge(a.fields,b.fields))
+_addSurfaces(a::_T, b::_T) where _T<:UTMData        = UTMData(a.EW.val, a.NS.val, a.depth.val + b.depth.val, merge(a.fields,b.fields))
+_addSurfaces(a::_T, b::_T) where _T<:CartData       = CartData(a.x.val, a.y.val, a.z.val + b.z.val, merge(a.fields,b.fields))
+_addSurfaces(a::_T, b::_T) where _T<:ParaviewData   = ParaviewData(a.x.val, a.y.val, a.z.val + b.z.val, merge(a.fields,b.fields))
+
+_subtractSurfaces(a::_T, b::_T) where _T<:GeoData        = GeoData(a.lon.val, a.lat.val, a.depth.val - b.depth.val, merge(a.fields,b.fields))
+_subtractSurfaces(a::_T, b::_T) where _T<:UTMData        = UTMData(a.EW.val, a.NS.val, a.depth.val - b.depth.val, merge(a.fields,b.fields))
+_subtractSurfaces(a::_T, b::_T) where _T<:CartData       = CartData(a.x.val, a.y.val, a.z.val - b.z.val, merge(a.fields,b.fields))
+_subtractSurfaces(a::_T, b::_T) where _T<:ParaviewData   = ParaviewData(a.x.val, a.y.val, a.z.val - b.z.val, merge(a.fields,b.fields))
+
+"""
+    remove_NaN_Surface!(Z::Array,X::Array,Y::Array)
+
+Removes NaN's from a grid `Z` by taking the closest points as specified by `X` and `Y`.
+"""
+function remove_NaN_Surface!(Z,X,Y)
+    @assert size(Z) == size(X) == size(Y)
+
+    # use nearest neighbour to interpolate data
+    id      = findall(isnan.(Z) .== false)
+    id_NaN  = findall(isnan.(Z))
+
+    coord   =   [X[id]'; Y[id]'];
+    kdtree  =   KDTree(coord; leafsize = 10);
+
+    points    = [X[id_NaN]'; Y[id_NaN]'];
+    idx,dist  = nn(kdtree, points);
+
+    Z[id_NaN] = Z[id[idx]]
+
+    return nothing
+end
+
+
+"""
+    Topo = drape_on_topo(Topo::GeoData, Data::GeoData)
+
+This drapes fields of a data set `Data` on the topography `Topo`
+
+"""
+function drape_on_topo(Topo::GeoData, Data::GeoData)
+
+    Lon,Lat,_   =   LonLatDepthGrid( Topo.lon.val[:,1,1], Topo.lat.val[1,:,1],Topo.depth.val[1,1,:]);
+
+    # use nearest neighbour to interpolate data
+    coord       =   [vec(Data.lon.val)'; vec(Data.lat.val)'];
+    kdtree      =   KDTree(coord; leafsize = 10);
+    points      =   [vec(Lon)';vec(Lat)'];
+    idx,_       =   nn(kdtree, points);
+
+    idx_out     = findall(  (Lon .<  minimum(Data.lon.val)) .| (Lon .>  maximum(Data.lon.val)) .|
+                            (Lat .<  minimum(Data.lat.val)) .| (Lat .>  maximum(Data.lat.val)) )
+
+    fields_new  = Topo.fields;
+    field_names = keys(Data.fields);
+
+    for i = 1:length(Data.fields)
+
+        if typeof(Data.fields[i]) <: Tuple
+
+            # vector or anything that contains more than 1 field
+            data_tuple = Data.fields[i]      # we have a tuple (likely a vector field), so we have to loop
+
+            data_array = zeros(typeof(data_tuple[1][1]),size(Topo.lon.val,1),size(Topo.lon.val,2),size(Topo.lon.val,3),length(Data.fields[i]));     # create a 3D array that holds the 2D interpolated values
+            unit_array = zeros(size(data_array));
+
+            for j=1:length(data_tuple)
+                data_field           =   data_tuple[j];
+                tmp                  =   data_array[:,:,:,1];
+                tmp                  =   data_field[idx]
+                data_array[:,:,:,j]  =   tmp
+            end
+
+            data_new    = tuple([data_array[:,:,:,c] for c in 1:size(data_array,4)]...)       # transform 4D matrix to tuple
+
+            # remove points outside domain
+            for j=1:length(data_tuple)
+                tmp           =   data_new[j];
+                tmp[idx_out] .= NaN
+                data_array[:,:,:,j]  =   tmp
+            end
+            data_new    = tuple([data_array[:,:,:,c] for c in 1:size(data_array,4)]...)       # transform 4D matrix to tuple
+
+        else
+
+            # scalar field
+            data_new        =   zeros(typeof(Data.fields[i][1]), size(Topo.lon.val,1),size(Topo.lon.val,2),size(Topo.lon.val,3));
+            data_new        =   Data.fields[i][idx]                                 # interpolate data field
+
+        end
+
+        # replace the one
+        new_field   =   NamedTuple{(field_names[i],)}((data_new,))                          # Create a tuple with same name
+        fields_new  =   merge(fields_new, new_field);                                       # replace the field in fields_new
+
+    end
+
+    Topo_new        =   GeoData(Topo.lon.val,Topo.lat.val,Topo.depth.val, fields_new)
+
+    return Topo_new
+end
+
+
+"""
+    drape_on_topo(Topo::CartData, Data::CartData)
+
+Drapes Cartesian Data on topography
+"""
+function drape_on_topo(Topo::CartData, Data::CartData)
+    Topo_lonlat = GeoData(ustrip.(Topo.x.val),ustrip.(Topo.y.val), ustrip.(Topo.z.val), Topo.fields )
+    Data_lonlat = GeoData(ustrip.(Data.x.val),ustrip.(Data.y.val), ustrip.(Data.z.val), Data.fields )
+
+    Topo_new_lonlat = drape_on_topo(Topo_lonlat, Data_lonlat)
+
+    Topo_new = CartData(Topo_new_lonlat.lon.val, Topo_new_lonlat.lat.val, Topo_new_lonlat.depth.val, Topo_new_lonlat.fields)
+
+    return Topo_new
+end
+
+
+"""
+
+This projects 3D points described by the vectors `lon_pt`,`lat_pt`, `depth_pt` to a GeoData surface `surf`
+
+"""
+function project_points_to_surface(surf::GeoData, lon_pt::Vector, lat_pt::Vector, depth_pt::Vector)
+
+end