GitHub

Surfaces

We have a number of functions to deal with horizontal surfaces, which are defined as GeoData or CartData structures that have 3 dimensional flat coordinate array. Flat implies that the resolution is n by m by 1, which is required to correctly display them in paraview. The most straightforward surface is the topography (which you can obtain with import_topo). Surfaces can be added and subtracted.

GeophysicalModelGenerator.fit_surface_to_points — Function
surf_new = fit_surface_to_points(surf::GeoData, lon_pt::Vector, lat_pt::Vector, depth_pt::Vector)

This fits the depth values of the surface surf to the depth value of the closest-by-points in (lon_pt,lat_pt, depth_pt)

source
surf_new = fit_surface_to_points(surf::CartData, lon_pt::Vector, lat_pt::Vector, depth_pt::Vector)

This fits the depth values of the surface surf to the depth value of the closest-by-points in (lon_pt,lat_pt, depth_pt)

source
GeophysicalModelGenerator.above_surface — Function
Above = above_surface(Data_LaMEM::LaMEM_grid, DataSurface_Cart::CartData)

Determines if points within the 3D LaMEM_grid structure are above the Cartesian surface DataSurface_Cart

source
above_surface(Data::GeoData, DataSurface::GeoData; above=true)

Returns a boolean array of size(Data.Lon), which is true for points that are above the surface DataSurface (or for points below if above=false).

This can be used, for example, to mask points above/below the Moho in a volumetric dataset or in a profile.

Example

First we create a 3D data set and a 2D surface:

julia> Lon,Lat,Depth   =   lonlatdepth_grid(10:20,30:40,(-300:25:0)km);
julia> Data            =   Depth*2;
julia> Data_set3D      =   GeoData(Lon,Lat,Depth,(Depthdata=Data,LonData=Lon))
GeoData
  size  : (11, 11, 13)
  lon   ϵ [ 10.0 : 20.0]
  lat   ϵ [ 30.0 : 40.0]
  depth ϵ [ -300.0 km : 0.0 km]
  fields: (:Depthdata, :LonData)
julia> Lon,Lat,Depth   =   lonlatdepth_grid(10:20,30:40,-40km);
julia> Data_Moho       =   GeoData(Lon,Lat,Depth+Lon*km, (MohoDepth=Depth,))
  GeoData
    size  : (11, 11, 1)
    lon   ϵ [ 10.0 : 20.0]
    lat   ϵ [ 30.0 : 40.0]
    depth ϵ [ -30.0 km : -20.0 km]
    fields: (:MohoDepth,)

Next, we intersect the surface with the data set:

julia> Above       =   above_surface(Data_set3D, Data_Moho);

Now, Above is a boolean array that is true for points above the surface and false for points below and at the surface.

source
Above = above_surface(Data_Cart::ParaviewData, DataSurface_Cart::ParaviewData; above=true)

Determines if points within the 3D Data_Cart structure are above the Cartesian surface DataSurface_Cart

source
Above = above_surface(Data_Cart::Union{Q1Data,CartData}, DataSurface_Cart::CartData; above=true)

Determines if points within the 3D Data_Cart structure are above the Cartesian surface DataSurface_Cart

source
Above = above_surface(Grid::CartGrid, DataSurface_Cart::CartData; above=true, cell=false)

Determines if points described by the Grid CartGrid structure are above the Cartesian surface DataSurface_Cart

source
GeophysicalModelGenerator.below_surface — Function
Below = below_surface(Data_LaMEM::LaMEM_grid, DataSurface_Cart::CartData)

Determines if points within the 3D LaMEM_grid structure are below the Cartesian surface DataSurface_Cart

source
Below = below_surface(Data::GeoData, DataSurface::GeoData)

Determines if points within the 3D Data structure are below the GeoData surface DataSurface

source
Below = below_surface(Grid::CartGrid, DataSurface_Cart::CartData)

Determines if points described by the `Grid` CartGrid structure are above the Cartesian surface `DataSurface_Cart`
source
Below = below_surface(Data_Cart::ParaviewData, DataSurface_Cart::ParaviewData)

Determines if points within the 3D DataCart structure are below the Cartesian surface DataSurfaceCart

source
Below = below_surface(Data_Cart::Union{CartData,Q1Data}, DataSurface_Cart::CartData, cell=false)

Determines if points within the 3D Data_Cart structure are below the Cartesian surface DataSurface_Cart

source
GeophysicalModelGenerator.interpolate_data_surface — Function
Surf_interp = interpolate_data_surface(V::ParaviewData, Surf::ParaviewData)

Interpolates a 3D data set V on a surface defined by Surf.

Example

julia> Data
ParaviewData
  size  : (33, 33, 33)
  x     ϵ [ -3.0 : 3.0]
  y     ϵ [ -2.0 : 2.0]
  z     ϵ [ -2.0 : 0.0]
  fields: (:phase, :density, :visc_total, :visc_creep, :velocity, :pressure, :temperature, :dev_stress, :strain_rate, :j2_dev_stress, :j2_strain_rate, :plast_strain, :plast_dissip, :tot_displ, :yield, :moment_res, :cont_res)
julia> surf
ParaviewData
  size  : (96, 96, 1)
  x     ϵ [ -2.9671875 : 3.2671875]
  y     ϵ [ -1.9791666666666667 : 1.9791666666666667]
  z     ϵ [ -1.5353766679763794 : -0.69925457239151]
  fields: (:Depth,)
julia> Surf_interp = interpolate_data_surface(Data, surf)
  ParaviewData
    size  : (96, 96, 1)
    x     ϵ [ -2.9671875 : 3.2671875]
    y     ϵ [ -1.9791666666666667 : 1.9791666666666667]
    z     ϵ [ -1.5353766679763794 : -0.69925457239151]
    fields: (:phase, :density, :visc_total, :visc_creep, :velocity, :pressure, :temperature, :dev_stress, :strain_rate, :j2_dev_stress, :j2_strain_rate, :plast_strain, :plast_dissip, :tot_displ, :yield, :moment_res, :cont_res)
source
Surf_interp = interpolate_data_surface(V::GeoData, Surf::GeoData)

Interpolates a 3D data set V on a surface defined by Surf

source