List of all functions

Data = CartData(Grid::CartGrid, fields::NamedTuple; y_val=0.0)

Returns a CartData set given a cartesian grid Grid and fields defined on that grid.

CartData(Grid::LaMEM_grid, fields::NamedTuple)

Creates a CartData struct from a LaMEM grid and from fields stored on that grid. Note that one needs to have a field Phases and optionally a field Temp to create LaMEM marker files.

CartData(xyz::Tuple{Array,Array,Array})

This creates a CartData struct if you have a Tuple with 3D coordinates as input.

Example

julia> data = CartData(xyz_grid(-10:10,-5:5,0))
CartData
    size    : (21, 11, 1)
    x       ϵ [ -10.0 km : 10.0 km]
    y       ϵ [ -5.0 km : 5.0 km]
    z       ϵ [ 0.0 km : 0.0 km]
    fields  : (:Z,)
  attributes: ["note"]

Structure that holds data for an orthogonal cartesian grid, which can be described with 1D vectors

Structure that stores info about a GMG Dataset, which is useful to collect a wide variety of datasets.

• Name :: String # Name of the dataset
• Type :: String # Volumetric, Surface, Point, Screenshot
• DirName :: String # Directory name or url of dataset
• active :: Bool # should this data be loaded or not?

GeoData(lld::Tuple{Array,Array,Array})

This creates a GeoData struct if you have a Tuple with 3D coordinates as input.

Example

julia> data = GeoData(lonlatdepth_grid(-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,)

ParaviewData(Grid::LaMEM_grid, fields::NamedTuple)

Creates a ParaviewData struct from a LaMEM grid and from fields stored on that grid. Note that one needs to have a field Phases and optionally a field Temp to create LaMEM marker files.

Structure that holds profile data (interpolated/projected on the profile)

struct ProfileData
    vertical        ::  Bool # vertical:true, horizontal:false
    start_lonlat    ::  Union{Nothing,Tuple{Float64,Float64}}
    end_lonlat      ::  Union{Nothing,Tuple{Float64,Float64}}
    depth           ::  Union{Nothing,Float64}
    VolData         ::  GeophysicalModelGenerator.GeoData
    SurfData        ::  Union{Nothing, NamedTuple}
    PointData       ::  Union{Nothing, NamedTuple}
end

Structure to store cross section data

ProjectionPoint(EW::Float64, NS::Float64, Zone::Int64, isnorth::Bool)

Defines a projection point used for map projections, by specifying UTM coordinates (EW/NS), UTM Zone and whether you are on the northern hemisphere

ProjectionPoint(; Lat=49.9929, Lon=8.2473)

Defines a projection point used for map projections, by specifying latitude and longitude

Q1Data(xyz::Tuple{Array,Array,Array})

This creates a Q1Data struct if you have a Tuple with 3D coordinates as input.

Example

julia> data = Q1Data(xyz_grid(-10:10,-5:5,0))
CartData
    size    : (21, 11, 1)
    x       ϵ [ -10.0 km : 10.0 km]
    y       ϵ [ -5.0 km : 5.0 km]
    z       ϵ [ 0.0 km : 0.0 km]
    fields  : (:Z,)
  attributes: ["note"]

Trench structure

Structure that defines the geometry of the trench and the slab.

Parameters

• Start - Start of the trench (x,y) coordinates

• End - End of the trench (x,y) coordinates

• n_seg - The number of segment through which the slab is discretize along the dip

• Length - The length of the slab

• Thickness - The thickness of the slab

• Lb - Critical distance through which apply the bending angle functions Lb ∈ [0,Length];

• θ_max - maximum angle of bending ∈ [0°,90°].

• direction - the direction of the dip The rotation of the coordinate system is done as such that the new X is parallel to the segment. Since the rotation is anticlockwise the coordinate y has specific values: direction tells if the subduction is directed along the positive or negative direction of the new y coordinate system. In practice, it apply an additional transformation to y by multiplying it with -1 or +1;

• d_decoupling - depth at which the slab is fully submerged into the mantle.

• type_bending - is the type of bending angle of the slab [:Linear, :Ribe]. The angle of slab changes as a function of l (∈ [0,Length]). l is the actual distance along the slab length from the trench. In case: - :Linear math θ(l) = ((θ_max - 0.0)/(Lb-0))*l; - :Ribe math θ(l) = θ_max*l^2*((3*Lb-2*l))/(Lb^3); which is taken from Ribe 2010 [Bending mechanics and mode selection in free subduction: a thin-sheet analysis]

  For l>Lb, θ(l) = θ_max;

• WeakzoneThickness - Thickness of the weakzone [km]

• WeakzonePhase - Phase of the weakzone

Converts a UTMData structure to a GeoData structure

Converts a GeoData structure to a UTMData structure

Returns 1D coordinate vectors of grid points and of marker locations for a regular spacing

Returns 1D coordinate vectors of grid points and of marker locations for a regular spacing

InterpolateDataFields2D_vecs(x_vec, y_vec, depth, fields_new, X, Y)

Interpolates a data field V on a 2D grid defined by UTM. Typically used for horizontal surfaces

This parses a LaMEM command line argument string and checks if the keyword exists there

value = ParseValue_LaMEM_InputFile(file,keyword,type; args::String=nothing)

Extracts a certain keyword from a LaMEM input file and convert it to a certain type. Optionally, you can also pass command-line arguments which will override the value read from the input file.

Example 1:

julia> nmark_z = ParseValue_LaMEM_InputFile("SaltModels.dat","nmark_z",Int64)

Example 2:

julia> nmark_z = ParseValue_LaMEM_InputFile("SaltModels.dat","nmark_z",Int64, args="-nel_x 128 -coord_x -4,4")

PetscBinaryWrite_Vec(filename, A)

Writes a vector A to disk, such that it can be read with PetscBinaryRead (which assumes a Big Endian type)

xrot, yrot, zrot = Rot3D(X::Number,Y::Number,Z::Number, cosStrikeAngle, sindStrikeAngle, cosDipAngle, sinDipAngle)

Perform rotation for a point in 3D space

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

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

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.

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

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

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

add_box!(Phase, Temp, Grid::AbstractGeneralGrid; xlim=Tuple{2}, [ylim=Tuple{2}], zlim=Tuple{2},
        Origin=nothing, StrikeAngle=0, DipAngle=0,
        phase = ConstantPhase(1),
        T=nothing,
        cell=false )

Adds a box with phase & temperature structure to a 3D model setup. This simplifies creating model geometries in geodynamic models

Parameters

• Phase - Phase array (consistent with Grid)
• Temp - Temperature array (consistent with Grid)
• Grid - grid structure (can be any of the grid types in GMG)
• xlim - left/right coordinates of box
• ylim - front/back coordinates of box [optional; if not specified we use the whole box]
• zlim - bottom/top coordinates of box
• Origin - the origin, used to rotate the box around. Default is the left-front-top corner
• StrikeAngle - strike angle of slab
• DipAngle - dip angle of slab
• phase - specifies the phase of the box. See ConstantPhase(),LithosphericPhases()
• T - specifies the temperature of the box. See ConstantTemp(),LinearTemp(),HalfspaceCoolingTemp(),SpreadingRateTemp(),LithosphericTemp()
• cell - if true, Phase and Temp are defined on centers

Examples

Example 1) Box with constant phase and temperature & a dip angle of 10 degrees:

julia> Grid = read_LaMEM_inputfile("test_files/SaltModels.dat")
LaMEM Grid:
  nel         : (32, 32, 32)
  marker/cell : (3, 3, 3)
  markers     : (96, 96, 96)
  x           ϵ [-3.0 : 3.0]
  y           ϵ [-2.0 : 2.0]
  z           ϵ [-2.0 : 0.0]
julia> Phases = zeros(Int32,   size(Grid.X));
julia> Temp   = zeros(Float64, size(Grid.X));
julia> add_box!(Phases,Temp,Grid, xlim=(0,500), zlim=(-50,0), phase=ConstantPhase(3), DipAngle=10, T=ConstantTemp(1000))
julia> Model3D = ParaviewData(Grid, (Phases=Phases,Temp=Temp)); # Create Cartesian model
julia> write_paraview(Model3D,"LaMEM_ModelSetup")           # Save model to paraview
1-element Vector{String}:
 "LaMEM_ModelSetup.vts"

Example 2) Box with halfspace cooling profile

julia> Grid = CartData(xyz_grid(-1000:10:1000,0,-660:10:0))
julia> Phases = zeros(Int32,   size(Grid));
julia> Temp   = zeros(Float64, size(Grid));
julia> add_box!(Phases,Temp,Grid, xlim=(0,500), zlim=(-50,0), phase=ConstantPhase(3), DipAngle=10, T=HalfspaceCoolingTemp(Age=30))
julia> Grid = addfield(Grid, (;Phases,Temp));       # Add to Cartesian model
julia> write_paraview(Grid,"LaMEM_ModelSetup")  # Save model to paraview
1-element Vector{String}:
 "LaMEM_ModelSetup.vts"

add_cylinder!(Phase, Temp, Grid::AbstractGeneralGrid; base=Tuple{3}, cap=Tuple{3}, radius=Tuple{1},
        phase = ConstantPhase(1).
        T=nothing, cell=false )

Adds a cylinder with phase & temperature structure to a 3D model setup. This simplifies creating model geometries in geodynamic models

Parameters

• Phase - Phase array (consistent with Grid)
• Temp - Temperature array (consistent with Grid)
• Grid - Grid structure (usually obtained with readLaMEMinputfile)
• base - center coordinate of bottom of cylinder
• cap - center coordinate of top of cylinder
• radius - radius of the cylinder
• phase - specifies the phase of the box. See ConstantPhase(),LithosphericPhases()
• T - specifies the temperature of the box. See ConstantTemp(),LinearTemp(),HalfspaceCoolingTemp(),SpreadingRateTemp()
• cell - if true, Phase and Temp are defined on cell centers

Example

Cylinder with constant phase and temperature:

julia> Grid = read_LaMEM_inputfile("test_files/SaltModels.dat")
LaMEM Grid:
  nel         : (32, 32, 32)
  marker/cell : (3, 3, 3)
  markers     : (96, 96, 96)
  x           ϵ [-3.0 : 3.0]
  y           ϵ [-2.0 : 2.0]
  z           ϵ [-2.0 : 0.0]
julia> Phases = zeros(Int32,   size(Grid.X));
julia> Temp   = zeros(Float64, size(Grid.X));
julia> add_cylinder!(Phases,Temp,Grid, base=(-1,-1,-1.5), cap=(1,1,-0.5), radius=0.25, phase=ConstantPhase(4), T=ConstantTemp(400))
julia> Model3D = ParaviewData(Grid, (Phases=Phases,Temp=Temp)); # Create Cartesian model
julia> write_paraview(Model3D,"LaMEM_ModelSetup")           # Save model to paraview
1-element Vector{String}:
 "LaMEM_ModelSetup.vts"

add_ellipsoid!(Phase, Temp, Grid::AbstractGeneralGrid; cen=Tuple{3}, axes=Tuple{3},
        Origin=nothing, StrikeAngle=0, DipAngle=0,
        phase = ConstantPhase(1).
        T=nothing, cell=false )

Adds an Ellipsoid with phase & temperature structure to a 3D model setup. This simplifies creating model geometries in geodynamic models

Parameters

• Phase - Phase array (consistent with Grid)
• Temp - Temperature array (consistent with Grid)
• Grid - LaMEM grid structure (usually obtained with readLaMEMinputfile)
• cen - center coordinates of sphere
• axes - semi-axes of ellipsoid in X,Y,Z
• Origin - the origin, used to rotate the box around. Default is the left-front-top corner
• StrikeAngle - strike angle of slab
• DipAngle - dip angle of slab
• phase - specifies the phase of the box. See ConstantPhase(),LithosphericPhases()
• T - specifies the temperature of the box. See ConstantTemp(),LinearTemp(),HalfspaceCoolingTemp(),SpreadingRateTemp()
• cell - if true, Phase and Temp are defined on cell centers

Example

Ellipsoid with constant phase and temperature, rotated 90 degrees and tilted by 45 degrees:

julia> Grid = read_LaMEM_inputfile("test_files/SaltModels.dat")
LaMEM Grid:
  nel         : (32, 32, 32)
  marker/cell : (3, 3, 3)
  markers     : (96, 96, 96)
  x           ϵ [-3.0 : 3.0]
  y           ϵ [-2.0 : 2.0]
  z           ϵ [-2.0 : 0.0]
julia> Phases = zeros(Int32,   size(Grid.X));
julia> Temp   = zeros(Float64, size(Grid.X));
julia> add_ellipsoid!(Phases,Temp,Grid, cen=(-1,-1,-1), axes=(0.2,0.1,0.5), StrikeAngle=90, DipAngle=45, phase=ConstantPhase(3), T=ConstantTemp(600))
julia> Model3D = ParaviewData(Grid, (Phases=Phases,Temp=Temp)); # Create Cartesian model
julia> write_paraview(Model3D,"LaMEM_ModelSetup")           # Save model to paraview
1-element Vector{String}:
 "LaMEM_ModelSetup.vts"

add_layer!(Phase, Temp, Grid::AbstractGeneralGrid; xlim=Tuple{2}, [ylim=Tuple{2}], zlim=Tuple{2},
        phase = ConstantPhase(1),
        T=nothing, cell=false )

Adds a layer with phase & temperature structure to a 3D model setup. The most common use would be to add a lithospheric layer to a model setup. This simplifies creating model geometries in geodynamic models

Parameters

• Phase - Phase array (consistent with Grid)
• Temp - Temperature array (consistent with Grid)
• Grid - grid structure (usually obtained with readLaMEMinputfile, but can also be other grid types)
• xlim - left/right coordinates of box
• ylim - front/back coordinates of box
• zlim - bottom/top coordinates of box
• phase - specifies the phase of the box. See ConstantPhase(),LithosphericPhases()
• T - specifies the temperature of the box. See ConstantTemp(),LinearTemp(),HalfspaceCoolingTemp(),SpreadingRateTemp()

Examples

Example 1) Layer with constant phase and temperature

julia> Grid = read_LaMEM_inputfile("test_files/SaltModels.dat")
LaMEM Grid:
  nel         : (32, 32, 32)
  marker/cell : (3, 3, 3)
  markers     : (96, 96, 96)
  x           ϵ [-3.0 : 3.0]
  y           ϵ [-2.0 : 2.0]
  z           ϵ [-2.0 : 0.0]
julia> Phases = zeros(Int32,   size(Grid.X));
julia> Temp   = zeros(Float64, size(Grid.X));
julia> add_layer!(Phases,Temp,Grid, zlim=(-50,0), phase=ConstantPhase(3), T=ConstantTemp(1000))
julia> Model3D = ParaviewData(Grid, (Phases=Phases,Temp=Temp)); # Create Cartesian model
julia> write_paraview(Model3D,"LaMEM_ModelSetup")           # Save model to paraview
1-element Vector{String}:
 "LaMEM_ModelSetup.vts"

Example 2) Box with halfspace cooling profile

julia> Grid = read_LaMEM_inputfile("test_files/SaltModels.dat")
julia> Phases = zeros(Int32,   size(Grid.X));
julia> Temp   = zeros(Float64, size(Grid.X));
julia> add_layer!(Phases,Temp,Grid, zlim=(-50,0), phase=ConstantPhase(3), T=HalfspaceCoolingTemp())
julia> Model3D = ParaviewData(Grid, (Phases=Phases,Temp=Temp)); # Create Cartesian model
julia> write_paraview(Model3D,"LaMEM_ModelSetup")           # Save model to paraview
1-element Vector{String}:
 "LaMEM_ModelSetup.vts"

    add_polygon!(Phase, Temp, Grid::AbstractGeneralGrid; xlim::Vector(), ylim=Vector(2), zlim=Vector(), phase = ConstantPhase(1), T=nothing, cell=false )

Adds a polygon with phase & temperature structure to a 3D model setup. This simplifies creating model geometries in geodynamic models

Parameters

• Phase - Phase array (consistent with Grid)
• Temp - Temperature array (consistent with Grid)
• Grid - Grid structure (usually obtained with readLaMEMinputfile)
• xlim - x-coordinate of the polygon points, same ordering as zlim, number of points unlimited
• ylim - y-coordinate, limitation in length possible (two values (start and stop))
• zlim - z-coordinate of the polygon points, same ordering as xlim, number of points unlimited
• phase - specifies the phase of the box. See ConstantPhase()
• T - specifies the temperature of the box. See ConstantTemp(),LinearTemp(),HalfspaceCoolingTemp(),SpreadingRateTemp()
• cell - if true, Phase and Temp are defined on cell centers

Example

Polygon with constant phase and temperature:

julia> Grid = read_LaMEM_inputfile("test_files/SaltModels.dat")
LaMEM Grid:
  nel         : (32, 32, 32)
  marker/cell : (3, 3, 3)
  markers     : (96, 96, 96)
  x           ϵ [-3.0 : 3.0]
  y           ϵ [-2.0 : 2.0]
  z           ϵ [-2.0 : 0.0]
julia> Phases = zeros(Int32,   size(Grid.X));
julia> Temp   = zeros(Float64, size(Grid.X));
julia> add_polygon!(Phase, Temp, Cart; xlim=[0.0,0.0, 1.6, 2.0],ylim=[0.0,0.8], zlim=[0.0,-1.0,-2.0,0.0], phase = ConstantPhase(8), T=ConstantTemp(30))
julia> Model3D = ParaviewData(Grid, (Phases=Phases,Temp=Temp)); # Create Cartesian model
julia> write_paraview(Model3D,"LaMEM_ModelSetup")           # Save model to paraview
1-element Vector{String}:
 "LaMEM_ModelSetup.vts"

add_slab!(Phase, Temp, Grid::AbstractGeneralGrid,  trench::Trench; phase = ConstantPhase(1), T = nothing, cell=false)

Adds a curved slab with phase & temperature structure to a 3D model setup.

Parameters

• Phase - Phase array (consistent with Grid)
• Temp - Temperature array (consistent with Grid)
• Grid - grid structure (can be any of the grid types in GMG)
• trench - Trench structure
• phase - specifies the phase of the box. See ConstantPhase(),LithosphericPhases()
• T - specifies the temperature of the box. See ConstantTemp(),LinearTemp(),HalfspaceCoolingTemp(),SpreadingRateTemp(),LithosphericTemp()
• cell - if true, Phase and Temp are defined on cells

Examples

Example 1) Slab

julia> x     = LinRange(0.0,1200.0,128);
julia> y     = LinRange(0.0,1200.0,128);
julia> z     = LinRange(-660,50,128);
julia> Cart  = CartData(xyz_grid(x, y, z));
julia> Phase = ones(Int64,size(Cart));
julia> Temp  = fill(1350.0,size(Cart));
# Define the trench:
julia> trench= Trench(Start = (400.0,400.0), End = (800.0,800.0), θ_max = 45.0, direction = 1.0, n_seg = 50, Length = 600.0, Thickness = 80.0, Lb = 500.0, d_decoupling = 100.0, type_bending =:Ribe)
julia> phase = LithosphericPhases(Layers=[5 7 88], Phases = [2 3 4], Tlab=nothing)
julia> TsHC  = HalfspaceCoolingTemp(Tsurface=20.0, Tmantle=1350, Age=30, Adiabat=0.4)
julia> add_slab!(Phase, Temp, Cart, trench, phase = phase, T = TsHC)

add_sphere!(Phase, Temp, Grid::AbstractGeneralGrid; cen=Tuple{3}, radius=Tuple{1},
        phase = ConstantPhase(1).
        T=nothing, cell=false )

Adds a sphere with phase & temperature structure to a 3D model setup. This simplifies creating model geometries in geodynamic models

Parameters

• Phase - Phase array (consistent with Grid)
• Temp - Temperature array (consistent with Grid)
• Grid - LaMEM grid structure (usually obtained with readLaMEMinputfile)
• cen - center coordinates of sphere
• radius - radius of sphere
• phase - specifies the phase of the box. See ConstantPhase(),LithosphericPhases()
• T - specifies the temperature of the box. See ConstantTemp(),LinearTemp(),HalfspaceCoolingTemp(),SpreadingRateTemp()
• cell - if true, Phase and Temp are defined on cell centers

Example

Sphere with constant phase and temperature:

julia> Grid = read_LaMEM_inputfile("test_files/SaltModels.dat")
LaMEM Grid:
  nel         : (32, 32, 32)
  marker/cell : (3, 3, 3)
  markers     : (96, 96, 96)
  x           ϵ [-3.0 : 3.0]
  y           ϵ [-2.0 : 2.0]
  z           ϵ [-2.0 : 0.0]
julia> Phases = zeros(Int32,   size(Grid.X));
julia> Temp   = zeros(Float64, size(Grid.X));
julia> add_sphere!(Phases,Temp,Grid, cen=(0,0,-1), radius=0.5, phase=ConstantPhase(2), T=ConstantTemp(800))
julia> Model3D = ParaviewData(Grid, (Phases=Phases,Temp=Temp)); # Create Cartesian model
julia> write_paraview(Model3D,"LaMEM_ModelSetup")           # Save model to paraview
1-element Vector{String}:
 "LaMEM_ModelSetup.vts"

add_stripes!(Phase, Grid::AbstractGeneralGrid;
    stripAxes       = (1,1,0),
    stripeWidth     =  0.2,
    stripeSpacing   =  1,
    Origin          =  nothing,
    StrikeAngle     =  0,
    DipAngle        =  10,
    phase           =  ConstantPhase(3),
    stripePhase     =  ConstantPhase(4),
    cell            = false)

Adds stripes to a pre-defined phase (e.g. added using add_box!)

Parameters

• Phase - Phase array (consistent with Grid)
• Grid - grid structure (usually obtained with readLaMEMinputfile, but can also be other grid types)
• stripAxes - sets the axis for which we want the stripes. Default is (1,1,0) i.e. X, Y and not Z
• stripeWidth - width of the stripe
• stripeSpacing - space between two stripes
• Origin - the origin, used to rotate the box around. Default is the left-front-top corner
• StrikeAngle - strike angle
• DipAngle - dip angle
• phase - specifies the phase we want to apply stripes to
• stripePhase - specifies the stripe phase
• cell - if true, Phase and Temp are defined on centers

Example

Example: Box with striped phase and constant temperature & a dip angle of 10 degrees:

julia> Grid = read_LaMEM_inputfile("test_files/SaltModels.dat")
LaMEM Grid:
  nel         : (32, 32, 32)
  marker/cell : (3, 3, 3)
  markers     : (96, 96, 96)
  x           ϵ [-3.0 : 3.0]
  y           ϵ [-2.0 : 2.0]
  z           ϵ [-2.0 : 0.0]
julia> Phases = zeros(Int32,   size(Grid.X));
julia> Temp   = zeros(Float64, size(Grid.X));
julia> add_box!(Phases,Temp,Grid, xlim=(0,500), zlim=(-50,0), phase=ConstantPhase(3), DipAngle=10, T=ConstantTemp(1000))
julia> add_stripes!(Phases, Grid, stripAxes=(1,1,1), stripeWidth=0.2, stripeSpacing=1, Origin=nothing, StrikeAngle=0, DipAngle=10, phase=ConstantPhase(3), stripePhase=ConstantPhase(4))
julia> Model3D = ParaviewData(Grid, (Phases=Phases,Temp=Temp)); # Create Cartesian model
julia> write_paraview(Model3D,"LaMEM_ModelSetup")           # Save model to paraview
1-element Vector{String}:
 "LaMEM_ModelSetup.vts"

V = addfield(V::AbstractGeneralGrid,field_name::String,data::Any)

Add Fields Data to GeoData or CartData

V = addfield(V::CartData,new_fields::NamedTuple)

Add new_fields fields to a CartData dataset

V = addfield(V::FEData,new_fields::NamedTuple; cellfield=false)

Add new_fields fields to a FEData dataset; set cellfield to true if the field is a cell field; otherwise it is a vertex field

V = addfield(V::GeoData,new_fields::NamedTuple)

Add new_fields fields to a GeoData dataset

V = addfield(V::Q1Data,new_fields::NamedTuple; cellfield=false)

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

out = average_q1(d::Array)

3D linear averaging of a 3D array

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

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

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

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

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

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

Below = below_surface(Data_Cart::ParaviewData, DataSurface_Cart::ParaviewData)

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

tags = cell_tags_from_gmsh(mesh::GmshDiscreteModel)

Returns a list with integers that are the tags for each of the cells

VolData_combined = combine_vol_data(VolData::NamedTuple; lat=nothing, lon=nothing, depth=nothing, dims=(100,100,100), dataset_preferred = 1)

This takes different volumetric datasets (specified in VolData) & merges them into a single one. You need to either provide the "reference" dataset within the NamedTuple (dataset_preferred), or the lat/lon/depth and dimensions of the new dataset.

θ = compute_bending_angle(θ_max,Lb,l,type)

function that computes the bending angle θ as a function of length along the slab l.

Parameters

θ_max = maximum bending angle Lb = length at which the function of bending is applied (Lb<=Length) l = current position within the slab type = type of bending [:Ribe,:Linear]

Phase = compute_phase(Phase, Temp, X, Y, Z, s::LithosphericPhases, Ztop)

or

Phase = compute_phase(Phase, Temp, Grid::AbstractGeneralGrid, s::LithosphericPhases)

This copies the layered lithosphere onto the Phase matrix.

Parameters

• Phase - Phase array
• Temp - Temperature array
• X - x-coordinate array (consistent with Phase and Temp)
• Y - y-coordinate array (consistent with Phase and Temp)
• Z - Vertical coordinate array (consistent with Phase and Temp)
• s - LithosphericPhases
• Ztop - Vertical coordinate of top of model box
• Grid - Grid structure (usually obtained with readLaMEMinputfile)

Top, Bot = compute_slab_surface(trench::Trench)

Computes the (x,z) coordinates of the slab top, bottom surface using the mid surface of the slab as reference.

Parameters

• trench - Trench structure that contains the relevant parameters

Method

It computes it by discretizing the slab surface in n_seg segments, and computing the average bending angle (which is a function of the current length of the slab). Next, it compute the coordinates assuming that the trench is at 0.0, and assuming a positive θ_max angle.

compute_thermal_structure(Temp, X, Y, Z, Phase, s::LinearWeightedTemperature)

Weight average along distance Do a weight average between two field along a specified direction Given a distance {could be any array, from X,Y} -> it increase from the origin the weight of F1, while F2 decreases. This function has been conceived for averaging the solution of Mckenzie and half space cooling model, but in can be used to smooth the temperature field from continent ocean: -> Select the boundary to apply; -> transform the coordinate such that dist represent the perpendicular direction along which you want to apply this smoothening and in a such way that 0.0 is the point in which the weight of F1 is equal to 0.0; -> Select the points that belongs to this area -> compute the thermal fields {F1} {F2} -> then modify F.

compute_thermal_structure(Temp, X, Y, Z, Phase, s::McKenzie_subducting_slab)

Compute the temperature field of a McKenzie_subducting_slab. Uses the analytical solution of McKenzie (1969) ["Speculations on the consequences and causes of plate motions"]. The functions assumes that the bottom of the slab is the coordinate Z=0. Internally the function shifts the coordinate.

Parameters

============================= Temp Temperature array

• X X Array
• Y Y Array
• Z Z Array
• Phase Phase array
• s McKenziesubductingslab

convert2CartData(d::GeoData, proj::ProjectionPoint)

Converts a GeoData structure to a CartData structure, which essentially transfers the dimensions to km

convert2CartData(d::UTMData, proj::ProjectionPoint)

Converts a UTMData structure to a CartData structure, which essentially transfers the dimensions to km

fe_data::FEData = convert2FEData(d::Q1Data)

Creates a Q1 FEM mesh from the Q1Data data which holds the vertex coordinates and cell/vertex fields

convert2UTMzone(d::CartData, proj::ProjectionPoint)

This transfers a CartData dataset to a UTMData dataset, that has a single UTM zone. The point around which we project is ProjectionPoint

convert2UTMzone(d::GeoData, p::ProjectionPoint)

Converts a GeoData structure to fixed UTM zone, around a given ProjectionPoint This useful to use real data as input for a cartesian geodynamic model setup (such as in LaMEM). In that case, we need to project map coordinates to cartesian coordinates. One way to do this is by using UTM coordinates. Close to the ProjectionPoint the resulting coordinates will be rectilinear and distance in meters. The map distortion becomes larger the further you are away from the center.

X,Y,Z = coordinate_grids(Data::CartData; cell=false)

Returns 3D coordinate arrays

X,Y,Z = coordinate_grids(Data::CartGrid; cell=false)

Returns 3D coordinate arrays

LON,LAT,Z = coordinate_grids(Data::GeoData; cell=false)

Returns 3D coordinate arrays

X,Y,Z = coordinate_grids(Data::LaMEM_grid; cell=false)

Returns 3D coordinate arrays

X,Y,Z = coordinate_grids(Data::ParaviewData; cell=false)

Returns 3D coordinate arrays

X,Y,Z = coordinate_grids(Data::Q1Data; cell=false)

Returns 3D coordinate arrays

EW,NS,Z = coordinate_grids(Data::UTMData; cell=false)

Returns 3D coordinate arrays

DatasetcountMap = countmap(DataSet::GeoData,field::String,stepslon::Int64,stepslat::Int64)

Takes a 2D GeoData struct and counts entries of field per predefined control area. field should only consist of 1.0s and 0.0s. The control area is defined by steplon and steplat. steplon is the number of control areas in longitude direction and steplat the number of control areas in latitude direction. The counts per control area are normalized by the highest count.

julia> Data_Faults         = GeoData(Lon3D,Lat3D,Faults,(Faults=Faults,))
GeoData 
    size      : (375, 208, 1)
    lon       ϵ [ -9.932408319802885 : 34.93985125012068]
    lat       ϵ [ 35.086096468211394 : 59.919210145128545]
    depth     ϵ [ 0.0 : 1.0]
    fields    : (:Faults,)

julia> steplon  = 125
julia> steplat  = 70
julia> countmap = countmap(Data_Faults,"Faults",steplon,steplat)

GeoData 
    size      : (124, 69, 1)
    lon       ϵ [ -9.751471789279 : 34.75891471959677]
    lat       ϵ [ 35.26604656731949 : 59.73926004602028]
    depth     ϵ [ 0.0 : 1.0]
    fields    : (:countmap,)

julia

Grid = create_CartGrid(; size=(), x = nothing, z = nothing, y = nothing, extent = nothing, CharDim = nothing)

Creates a 1D, 2D or 3D cartesian grid of given size. Grid can be created by defining the size and either the extent (length) of the grid in all directions, or by defining start & end points (x,y,z). If you specify CharDim (a structure with characteristic dimensions created with GeoParams.jl), we will nondimensionalize the grd before creating the struct.

Spacing is assumed to be constant in a given direction

This can also be used for staggered grids, as we also create 1D vectors for the central points. The points you indicate in size are the corner points.

Note: since this is mostly for solid Earth geoscience applications, the second dimension is called z (vertical)

Examples

====

A basic case with non-dimensional units:

julia> Grid = create_CartGrid(size=(10,20),x=(0.,10), z=(2.,10))
Grid{Float64, 2}
           size: (10, 20)
         length: (10.0, 8.0)
         domain: x ∈ [0.0, 10.0], z ∈ [2.0, 10.0]
 grid spacing Δ: (1.1111111111111112, 0.42105263157894735)

An example with dimensional units:

julia> CharDim = GEO_units()
julia> Grid    = create_CartGrid(size=(10,20),x=(0.0km, 10km), z=(-20km, 10km), CharDim=CharDim)
CartGrid{Float64, 2}
           size: (10, 20)
         length: (0.01, 0.03)
         domain: x ∈ [0.0, 0.01], z ∈ [-0.02, 0.01]
 grid spacing Δ: (0.0011111111111111111, 0.0015789473684210528)

create_partitioning_file(LaMEM_input::String, NumProc::Int64; LaMEM_dir::String=pwd(), LaMEM_options::String="", MPI_dir="", verbose=true)

This executes LaMEM for the input file LaMEM_input & creates a parallel partitioning file for NumProc processors. The directory where the LaMEM binary is can be specified; if not it is assumed to be in the current directory. Likewise for the mpiexec directory (if not specified it is assumed to be available on the command line).

create_profile_volume!(Profile::ProfileData, VolData::AbstractGeneralGrid; DimsVolCross::NTuple=(100,100), Depth_extent=nothing)

Creates a cross-section through a volumetric 3D dataset VolData with the data supplied in Profile. Depth_extent can be the minimum & maximum depth for vertical profiles

cross_section(DataSet::AbstractGeneralGrid; dims=(100,100), Interpolate=false, Depth_level=nothing, Lat_level=nothing, Lon_level=nothing, Start=nothing, End=nothing, Depth_extent=nothing, section_width=50km)

Creates a cross-section through a 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.
• Depending on the type of input data (volume, surface or point data), cross sections will be created in a different manner:

1. Volume data: data will be interpolated or directly extracted from the data set.
2. Surface data: surface data will be interpolated or directly extracted from the data set
3. Point data: data will be projected to the chosen profile. Only data within a chosen distance (default is 50 km) will be used

• Interpolate indicates whether we want to simply extract the data from the data set (default) or whether we want to linearly interpolate it on a new grid, which has dimensions as specified in dims NOTE: THIS ONLY APPLIES TO VOLUMETRIC AND SURFACE DATA SETS
• 'section_width' indicates the maximal distance within which point data will be projected to the profile

Example:

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 
  size  : (11, 11, 1)
  lon   ϵ [ 10.0 : 20.0]
  lat   ϵ [ 30.0 : 40.0]
  depth ϵ [ -100.0 km : -100.0 km]
  fields: (:Depthdata, :LonData, :Velocity)

function cross_section_points(P::GeoData; Depth_level=nothing, Lat_level=nothing, Lon_level=nothing, Start=nothing, End=nothing, section_width=50 )

Creates a projection of separate points (saved as a GeoData object) onto a chosen plane. Only points with a maximum distance of section_width are taken into account

crosssectionsurface(Surface::GeoData; dims=(100,), Interpolate=false, Depthlevel=nothing; Latlevel=nothing; Lon_level=nothing; Start=nothing, End=nothing )

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.

• 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.

Example:

julia> Lon,Lat,Depth   =   lonlatdepth_grid(10:20,30:40,-50km);
julia> Data            =   Depth*2;                # some data
julia> Vx,Vy,Vz        =   ustrip(Data*3),ustrip(Data*4),ustrip(Data*5);
julia> Data_set2D      =   GeoData(Lon,Lat,Depth,(Depth=Depth,));
julia> Data_cross      =   cross_section_surface(Data_set2D, Lat_level =15)
GeoData
  size      : (100,)
  lon       ϵ [ 10.0 : 20.0]
  lat       ϵ [ 15.0 : 15.0]
  depth     ϵ [ NaN : NaN]
  fields    : (:Depth,)
  attributes: ["note"]

crosssectionvolume(Volume::AbstractGeneralGrid; dims=(100,100), Interpolate=false, Depthlevel=nothing; Latlevel=nothing; Lonlevel=nothing; Start=nothing, End=nothing, Depthextent=nothing )

Creates a cross-section through a volumetric (3D) 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.
• When both Start=(lon,lat) & End=(lon,lat) are given, one can also provide a the depth extent of the profile by providing Depthextent=(depthmin,depth_max)
• Interpolate indicates whether we want to simply extract the data from the 3D volume (default) or whether we want to linearly interpolate it on a new grid, which has dimensions as specified in dims
• Depth_extent is an optional parameter that can indicate the depth extent over which you want to interpolate the vertical cross-section. Default is the full vertical extent of the 3D dataset

Example:

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_volume(Data_set3D, Depth_level=-100km)
GeoData
  size  : (11, 11, 1)
  lon   ϵ [ 10.0 : 20.0]
  lat   ϵ [ 30.0 : 40.0]
  depth ϵ [ -100.0 km : -100.0 km]
  fields: (:Depthdata, :LonData, :Velocity)

distance_to_linesegment(p::NTuple{2,_T}, v::NTuple{2,_T}, w::NTuple{2,_T})

Computes the distance normal distance from a point p to a line segment defined by the points v and w.

download_data(url::String, local_filename="temp.dat"; dir=pwd(), maxattempts=5 )

Downloads a remote dataset with name url from a remote location and saves it to the current directory. If download fails, we make maxattempts attempts before giving up.

Example

julia> url  = "https://seafile.rlp.net/f/10f867e410bb4d95b3fe/?dl=1";
julia> download_data(url)
"/Users/kausb/.julia/dev/GeophysicalModelGenerator/temp.dat"

drape_on_topo(Topo::CartData, Data::CartData)

Drapes Cartesian Data on topography

Topo = drape_on_topo(Topo::GeoData, Data::GeoData)

This drapes fields of a data set Data on the topography Topo

extract_ProfileData!(Profile::ProfileData,VolData::GeoData, SurfData::NamedTuple, PointData::NamedTuple; DimsVolCross=(100,100),Depth_extent=nothing,DimsSurfCross=(100,),section_width=50)

Extracts data along a vertical or horizontal profile

extract_ProfileData(ProfileCoordFile::String,ProfileNumber::Int64,DataSetFile::String; DimsVolCross=(100,100),DepthVol=nothing,DimsSurfCross=(100,),WidthPointProfile=50km)

This is a convenience function (mostly for backwards compatibility with the MATLAB GUI) that loads the data from file & projects it onto a profile

extract_subvolume(V::CartData; Interpolate=false, X_level=nothing, Y_level=nothing, Z_level=nothing, dims=(50,50,50))

Extract or "cuts-out" a piece of a 2D or 3D GeoData set, defined by Lon, Lat and Depth coordinates.

This is useful if you are only interested in a part of a much bigger larger data set.

• Lon_level,Lat_level and Depth_level should be tuples that indicate (minimum_value, maximum_value) along the respective direction. If not specified we use the full range.
• By default, Interpolate=false and we find the closest indices within the data set (so your new data set will not go exactly from minimum to maximum).
• Alternatively, if Interpolate=true we interpolate the data onto a new grid that has dimensions dims. This can be useful to compare data sets that are originally given in different resolutions.

3D Example with no interpolation:

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)))
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, :Velocity)
julia> Data_extracted = extract_subvolume(Data_set3D,Lon_level=(10,12),Lat_level=(35,40))
GeoData
  size  : (3, 6, 13)
  lon   ϵ [ 10.0 : 12.0]
  lat   ϵ [ 35.0 : 40.0]
  depth ϵ [ -300.0 km : 0.0 km]
  fields: (:Depthdata, :LonData, :Velocity)

By default it extracts the data points closest to the area defined by Lonlevel/Latlevel/Depth_level.

2D Example along a cross-section through 3D data:

julia> X,Y,Z = xyz_grid(10:20,30:40,-300:25:0);
julia> Data = Z.*2
julia> Data_Int = Int64.(Data)
julia> DataSet_Cart = CartData(X,Y,Z,(Data=Data,Data_Int=Data_Int, Velocity=(X,Y,Z)))

julia> Data_cross = cross_section(DataSet_Cart, Start=(11.0,35), End=(19, 39.0))
CartData
    size    : (100, 100, 1)
    x       ϵ [ 11.0 : 19.0]
    y       ϵ [ 35.0 : 39.0]
    z       ϵ [ -300.0 : 0.0]
    fields  : (:Data, :Data_Int, :Velocity, :FlatCrossSection)
  attributes: ["note"]

julia> Data_extracted = extract_subvolume(Data_cross, X_level=(1,7), Z_level=(-200,-100))
  CartData
      size    : (50, 50, 1)
      x       ϵ [ 11.894427190999917 : 17.260990336999413]
      y       ϵ [ 35.44721359549995 : 38.130495168499706]
      z       ϵ [ -200.0 : -100.0]
      fields  : (:FlatCrossSection, :Data, :Data_Int, :Velocity)
    attributes: ["note"]
julia> typeof(Data_extracted.fields.Data_Int)
    Array{Int64, 3}

extract_subvolume(V::GeoData; Interpolate=false, Lon_level=nothing, Lat_level=nothing, Depth_level=nothing, dims=(50,50,50))

Extract or "cuts-out" a piece of a 2D or 3D GeoData set, defined by Lon, Lat and Depth coordinates.

This is useful if you are only interested in a part of a much bigger larger data set.

• Lon_level,Lat_level and Depth_level should be tuples that indicate (minimum_value, maximum_value) along the respective direction. If not specified we use the full range.
• By default, Interpolate=false and we find the closest indices within the data set (so your new data set will not go exactly from minimum to maximum).
• Alternatively, if Interpolate=true we interpolate the data onto a new grid that has dimensions dims. This can be useful to compare data sets that are originally given in different resolutions.

3D Example with no interpolation:

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)))
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, :Velocity)
julia> Data_extracted = extract_subvolume(Data_set3D,Lon_level=(10,12),Lat_level=(35,40))
GeoData
  size  : (3, 6, 13)
  lon   ϵ [ 10.0 : 12.0]
  lat   ϵ [ 35.0 : 40.0]
  depth ϵ [ -300.0 km : 0.0 km]
  fields: (:Depthdata, :LonData, :Velocity)

By default it extracts the data points closest to the area defined by Lonlevel/Latlevel/Depth_level.

3D Example with interpolation:

Alternatively, you can also interpolate the data onto a new grid:

julia> Data_extracted = extract_subvolume(Data_set3D,Lon_level=(10,12),Lat_level=(35,40), Interpolate=true, dims=(50,51,52))
GeoData
  size  : (50, 51, 52)
  lon   ϵ [ 10.0 : 12.0]
  lat   ϵ [ 35.0 : 40.0]
  depth ϵ [ -300.0 km : 0.0 km]
  fields: (:Depthdata, :LonData, :Velocity)

find_slab_distance!(ls, d, X,Y,Z, trench::Trench)

Function that finds the perpendicular distance to the top and bottom of the slab d, and the current length of the slab l.

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)

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)

flatten_cross_section(V::CartData)

Takes a diagonal 3D crosssection and flattens it to be converted to a 2D Grid by createCartGrid

Example

Grid                    = create_CartGrid(size=(100,100,100), x=(0.0km, 99.9km), y=(-10.0km, 20.0km), z=(-40km,4km));
X,Y,Z                   = xyz_grid(Grid.coord1D...);
DataSet                 = CartData(X,Y,Z,(Depthdata=Z,));

Data_Cross              = cross_section(DataSet, dims=(100,100), Interpolate=true, Start=(ustrip(Grid.min[1]),ustrip(Grid.max[2])), End=(ustrip(Grid.max[1]), ustrip(Grid.min[2])))

x_new = flatten_cross_section(Data_Cross)

# This flattened cross_section can be added to original Data_Cross by addfield()

Data_Cross = addfield(Data_Cross,"FlatCrossSection", x_new)
CartData
    size    : (100, 100, 1)
    x       ϵ [ 0.0 : 99.9]
    y       ϵ [ -10.0 : 20.0]
    z       ϵ [ -40.0 : 4.0]
    fields  : (:Depthdata, :FlatCrossSection)
  attributes: ["note"]

flatten_cross_section(V::GeoData)
This function takes a 3D cross section through a GeoData structure and computes the distance along the cross section for later 2D processing/plotting
```julia-repl
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, Start=(10,30),End=(20,40))
julia> x_profile        =   flatten_cross_section(Data_cross)
julia> Data_cross      =   addfield(Data_cross,"x_profile",x_profile)

```

Data = flip(Data::GeoData, dimension=3)

This flips the data in the structure in a certain dimension (default is z [3])

nProcX,nProcY,nProcZ, xc,yc,zc, nNodeX,nNodeY,nNodeZ = get_processor_partitioning(filename; is64bit=false)

Reads a LaMEM processor partitioning file, used to create marker files, and returns the parallel layout. By default this is done for a 32bit PETSc installation, which will fail if you actually use a 64bit version.

Data = getlonlatdepthmag_QuakeML(filename::String)

Extracts longitude, latitude, depth and magnitude from a QuakeML file that has been e.g. downloaded from ISC. The data is then returned in GeoData format.

inpoly(PolyX::Vector, PolyY::Vector, x::Number, y::Number, iSteps::Vector, jSteps::)

Checks if a point given by x and y is in or on (both cases return true) a polygon given by PolyX and PolyY, iSteps and jSteps provide the connectivity between the polygon edges. This function should be used through inpolygon!().

inpoly_fast(PolyX::Vector, PolyY::Vector, x::Number, y::Number, iSteps::Vector, jSteps::)

Faster version of inpoly() but will miss some points that are on the edge of the polygon.

inpolygon!(INSIDE::Matrix, PolyX::Vector, PolyY::Vector, X::Matrix, Y::Matrix; fast=false)

Checks if points given by matrices X and Y are in or on (both cases return true) a polygon given by PolyX and PolyY. Boolean fast will trigger faster version that may miss points that are exactly on the edge of the polygon. Speedup is a factor of 3.

inpolygon!(inside::Vector, PolyX::Vector, PolyY::Vector, x::Vector, y::Vector; fast=false)

Same as above but inside, X and Y and are vectors.

interpolate_data_fields_cross_section(V::CartData, X,Y,Z,Xcross)

Interpolates data fields along a cross-section defined by Xcross and Z

Surf_interp = interpolate_data_surface(V::GeoData, Surf::GeoData)

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

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)

Data_interp = interpolate_datafields(V::AbstractGeneralGrid, Lon, Lat, Depth)

Interpolates a data field V on a grid defined by Lon,Lat,Depth

Example

julia> x        =   0:2:10
julia> y        =   -5:5
julia> z        =   -10:2:2
julia> X,Y,Z    =   xyz_grid(x, y, z);
julia> Data     =   Z
julia> Data_set1=   CartData(X,Y,Z, (FakeData=Data,Data2=Data.+1.))
CartData
    size    : (6, 11, 7)
    x       ϵ [ 0.0 km : 10.0 km]
    y       ϵ [ -5.0 km : 5.0 km]
    z       ϵ [ -10.0 km : 2.0 km]
    fields  : (:FakeData, :Data2)
  attributes: ["note"]

julia> X,Y,Z    =   xyz_grid(0:4:10, -1:.1:1, -5:.1:1 );
julia> Data_set2= interpolate_datafields(Data_set1, X,Y,Z)

interpolate_datafields(V::UTMData, EW, NS, Depth)

Interpolates a data field V on a grid defined by UTM,Depth

interpolate_datafields_2D(V::CartData, X, Y)

Interpolates a data field V on a 2D CartData grid defined by X,Y. Typically used for horizontal surfaces

interpolate_datafields_2D(Original::CartData, New::CartData; Rotate=0.0, Translate=(0,0,0), Scale=(1.0,1.0,1.0))

Interpolates a data field Original on a 2D CartData grid New. Typically used for horizontal surfaces.

Note: Original should have orthogonal coordinates. If it has not, e.g., because it was rotated, you'll have to specify the angle Rotate that it was rotated by

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.

interpolate_datafields_2D(V::GeoData, Lon, Lat)

Interpolates a data field V on a 2D grid defined by Lon,Lat. Typically used for horizontal surfaces

interpolate_datafields_2D(V::UTMData, EW, NS)

Interpolates a data field V on a 2D grid defined by UTM. Typically used for horizontal surfaces

issurf = is_surface(surf::AbstractGeneralGrid)

Returns true if surf is a horizontal 3D surface.

inside = isinside_closed_STL(mesh::Mesh, Pt, eps=1e-3)

Determine whether a point Pt is inside a 3D closed triangular *.stl surface or not.

This implements the winding number method, following the python code: https://github.com/marmakoide/inside-3d-mesh

This again is described in the following paper by Alec Jacobson, Ladislav Kavan and Olga Sorkine-Hornung.

lithostatic_pressure!(Plithos::Array, Density::Array, dz::Number; g=9.81)

Computes lithostatic pressure from a 3D density array, assuming constant soacing dz in vertical direction. Optionally, the gravitational acceleration g can be specified.

load_GMG(filename::String, dir=pwd(); maxattempts=5)

Loads a GeoData/CartData/UTMData data set from jld2 file filename Note: the filename can also be a remote url, in which case we first download that file to a temporary directory before opening it. We make maxattempts attempts to download it before giving up.

Example 1 - Load local file

julia> data = load_GMG("test")
GeoData 
  size      : (4, 3, 3)
  lon       ϵ [ 1.0 : 10.0]
  lat       ϵ [ 11.0 : 19.0]
  depth     ϵ [ -20.0 : -10.0]
  fields    : (:DataFieldName,)
  attributes: ["note"]

Example 2 - remote download

julia> url  = "https://seafile.rlp.net/f/10f867e410bb4d95b3fe/?dl=1"
julia> load_GMG(url)
GeoData 
  size      : (149, 242, 1)
  lon       ϵ [ -24.875 : 35.375]
  lat       ϵ [ 34.375 : 71.375]
  depth     ϵ [ -11.76 : -34.7]
  fields    : (:MohoDepth,)
  attributes: ["author", "year"]

data::NamedTuple = load_GMG(data::GMG_Dataset)

Loads a dataset specified in GMG_Dataset data and returns it as a named tuple

Data = load_GMG(Datasets::Vector{GMG_Dataset})

This loads all the active datasets in Datasets, and returns a NamedTuple with Volume, Surface, Point, Screenshot and Topography data

Datasets = load_dataset_file(file_datasets::String)

This loads a CSV textfile that lists datasets, which is expected to have the following format:

• Name,Location,Type, [Active]
• AlpArray,./Seismicity/ALPARRAY/AlpArraySeis.jld2,Point, true
• Plomerova2022,https://seafile.rlp.net/f/abccb8d3302b4ef5af17/?dl=1,Volume

Note that the first line of the file is skipped.

Here, the meaning of the variables is:

• Name: The name of the dataset to be loaded
• Location: the location of the file (directory and filename) on your local machine, or an url where we can download the file from the web. The url is expected to start with "http".
• Type: type of the dataset (Volume, Surface, Point, Screenshot)
• Active: Do we want this file to be loaded or not? Optional parameter that defaults to true

Lon, Lat, Depth = lonlatdepth_grid(Lon::Any, Lat::Any, Depth:Any)

Creates 3D arrays of Lon, Lat, Depth from 1D vectors or numbers

Example 1: Create 3D grid

julia> Lon,Lat,Depth =  lonlatdepth_grid(10:20,30:40,(-10:-1)km);
julia> size(Lon)
(11, 11, 10)

Example 2: Create 2D lon/lat grid @ a given depth

julia> Lon,Lat,Depth =  lonlatdepth_grid(10:20,30:40,-50km);
julia> size(Lon)
(11, 11)

Example 3: Create 2D lon/depth grid @ a given lat

julia> Lon,Lat,Depth =  lonlatdepth_grid(10:20,30,(-10:-1)km);
julia> size(Lon)
(11, 11)

Example 4: Create 1D vertical line @ a given lon/lat point

julia> Lon,Lat,Depth =  lonlatdepth_grid(10,30,(-10:-1)km);
julia> size(Lon)
(10, )

make_paraview_collection(; dir=pwd(), pvd_name=nothing, files=nothing, file_extension = ".vts", time = nothing)

In case one has a list of *.vtk files, this routine creates a *.pvd file that can be opened in Paraview. This is useful if you previously saved vtk files but didnt save it as a collection in the code itself.

Optional options

• dir: directory where the *.vtk are stored.
• pvd_name: filename of the resulting *.pvd file without extension; if not specified, full_simulation is used.
• files: Vector of the *.vtk files without extension; if not specified, all *.vtk files in the directory are used.
• file_extension: file extension of the vtk files. Default is .vts but all vt* work.
• time: Vector of the timesteps; if not specified, pseudo time steps are assigned.

makevolctopo(Grid::LaMEM_grid; center::Array{Float64, 1}, height::Float64, radius::Float64, crater::Float64, base=0.0m, background=nothing)

Creates a generic volcano topography (cones and truncated cones)

Parameters

• Grid - LaMEM grid (created by readLaMEMinputfile)
• center - x- and -coordinates of center of volcano
• height - height of volcano
• radius - radius of volcano

Optional Parameters

• crater - this will create a truncated cone and the option defines the radius of the flat top
• base - this sets the flat topography around the volcano
• background - this allows loading in a topography and only adding the volcano on top (also allows stacking of several cones to get a volcano with different slopes)

Example

Cylinder with constant phase and temperature:

julia> Grid = read_LaMEM_inputfile("test_files/SaltModels.dat")
LaMEM Grid:
  nel         : (32, 32, 32)
  marker/cell : (3, 3, 3)
  markers     : (96, 96, 96)
  x           ϵ [-3.0 : 3.0]
  y           ϵ [-2.0 : 2.0]
  z           ϵ [-2.0 : 0.0]
julia> Topo = make_volc_topo(Grid, center=[0.0,0.0], height=0.4, radius=1.5, crater=0.5, base=0.1)
CartData
    size    : (33, 33, 1)
    x       ϵ [ -3.0 : 3.0]
    y       ϵ [ -2.0 : 2.0]
    z       ϵ [ 0.1 : 0.4]
    fields  : (:Topography,)
  attributes: ["note"]
julia> Topo = make_volc_topo(Grid, center=[0.0,0.0], height=0.8, radius=0.5, crater=0.0, base=0.4, background=Topo.fields.Topography)
CartData
    size    : (33, 33, 1)
    x       ϵ [ -3.0 : 3.0]
    y       ϵ [ -2.0 : 2.0]
    z       ϵ [ 0.1 : 0.8]
    fields  : (:Topography,)
  attributes: ["note"]
julia> write_paraview(Topo,"VolcanoTopo")           # Save topography to paraview
Saved file: VolcanoTopo.vts

meshgrid(vx,vy,vz)

Computes an (x,y,z)-grid from the vectors (vx,vy,vz). For more information, see the MATLAB documentation.

movie_from_images(; dir=pwd(), file=nothing, outfile=nothing, framerate=10, copy_to_current_dir=true, type=:mp4_default, collect=true)

The typical way to create animations with Paraview is to use the Save Animation option to save a series of *.png images.

This function combines these images to an *.mp4 movie.

Optional options

• dir: directory where the images are stored.
• file: filename of the image series without extension and numbers. Required if >1 image series is stored in the same directory. By default we reconstruct this name from the available files.
• outfile: filename of the resulting movie without extension; if not specified, file is used.
• framerate: number of frames/second.
• copy_to_current_dir: copies the final movie to the current directory if true (default); otherwise it will be stored in dir.
• type: type of movie that is created; possible options are:
  • :mp4_default: Default option that saves a well-compressed mp4 movie that works well for us on ipad and embedded in a powerpoint presentation.
  • :mov_hires: Higher-resolution quicktime movie (larger filesize & not compatible with windows)
• collect: suppresses output of FFMPEG if true (default).

pvd = movie_paraview(; name="Movie", pvd=pvd, Finalize::Bool=false, Initialize::Bool=true)

If you want to make a movie of your data set, you can use this routine to initialize and to finalize the movie-file. It will create a *.pvd file, which you can open in Paraview

Individual timesteps are added to the movie by passing pvd and the time of the timestep to the write_paraview routine.

Example

Usually this is used inside a *.jl script, as in this pseudo-example:

movie = movie_paraview(name="Movie", Initialize=true)
for itime=1:10
    name = "test"*string(itime)
    movie = write_paraview(Data, name, pvd=movie, time=itime)
end
movie_paraview(pvd=movie, Finalize=true)

ind = nearest_point_indices(X::Array,Y::Array,Z::Array, X_pt::Vector,Y_pt::Vector,Z_pt::Vector)

Returns the index of the nearest point in (X_pt,Y_pt,Z_pt) to (X,Y,Z) and returns the index

ind = nearest_point_indices(X::Array,Y::Array, X_pt::Vector, Y_pt::Vector)

Returns the index of the nearest point in (X_pt,Y_pt) to (X,Y) and returns the index

ind = nearest_point_indices(X::Array, X_pt::Vector)

Returns the index of the nearest point in (X_pt) to (X) and returns the index

parse_columns_CSV(data_file, num_columns)

This parses numbers from CSV file that is read in with CSV.File. That is useful in case the CSV files has tables that contain both strings (e.g., station names) and numbers (lat/lon/height) and you are only interested in the numbers

Example

This example assumes that the data starts at line 18, that the columns are separated by spaces, and that it contains at most 4 columns with data:

julia> using CSV
julia> data_file        =   CSV.File("FileName.txt",datarow=18,header=false,delim=' ')
julia> data = parse_columns_CSV(data_file, 4)

inside = point_in_tetrahedron(p::_T, a::_T, b::_T, c::_T, d::_T, tol=1e-10)

Determines if a point p is inside a tetrahedron specified by a,b,c,d or not

count = point_to_nearest_grid(pt_x,pt_y,pt_z, X,Y,Z; radius_factor=1)

This uses nearest neighbour interpolation to count how many points (given by pt_x,pt_y,pt_z coordinate vectors) are in the vicinity of 3D grid point specified by X,Y,Z 3D coordinate arrays, with regular spacing (Δx,Δy,Δz). The search radius is R=radius_factor*(Δx² + Δy² + Δz²)^(1/3)

Grid_counts = point_to_nearest_grid(pt_x,pt_y,pt_z, Grid::CartData; radius_factor=1)

Uses nearest neighbour interpolation to count how many points (given by pt_x,pt_y,pt_z coordinate vectors) are in the vicinity of 3D CartGrid specified by Grid. The search radius is R=radius_factor*(Δx² + Δy² + Δz²)^(1/3)

Grid_counts is Grid but with an additional field Count that has the number of hits

Grid_counts = point_to_nearest_grid(pt_x,pt_y,pt_z, Grid::GeoData; radius_factor=1)

Uses nearest neighbour interpolation to count how many points (given by pt_x,pt_y,pt_z coordinate vectors) are in the vicinity of 3D GeoData specified by Grid. The search radius is R=radius_factor*(Δx² + Δy² + Δz²)^(1/3)

Grid_counts is Grid but with an additional field Count that has the number of hits

Grid_counts = point_to_nearest_grid(Point::CartData, Grid::CartData; radius_factor=1)

Uses nearest neighbour interpolation to count how many points (given by Point) are in the vicinity of a 3D Grid. The search radius is R=radius_factor*(Δx² + Δy² + Δz²)^(1/3)

Point should have 1D coordinate vectors

Grid_counts is Grid but with an additional field Count that has the number of hits

Grid_counts = point_to_nearest_grid(Point::GeoData, Grid::GeoData; radius_factor=1)

Uses nearest neighbour interpolation to count how many points (given by Point) are in the vicinity of a 3D Grid. The search radius is R=radius_factor*(Δx² + Δy² + Δz²)^(1/3)

Point should have 1D coordinate vectors

Grid_counts is Grid but with an additional field Count that has the number of hits

d_cart = project_CartData(d_cart::CartData, d::GeoData, p::ProjectionPoint)

Projects all datafields from the GeoData struct d to the CartData struct d_cart, around the projection point p. d_cart must be an orthogonal cartesian grid (deformed doesn't work; use convert2CartData(d, proj), where proj is a projection point in that case).

Note:

• If d_cart and d are horizontal surfaces (3rd dimension has size==1), it also interpolates the depth coordinate.

d_cart = project_CartData(d_cart::CartData, d::GeoData, p::ProjectionPoint)

Projects all datafields from the GeoData struct d to the CartData struct d_cart, around the projection point p. d_cart must be an orthogonal cartesian grid (deformed doesn't work; use convert2CartData(d, proj), where proj is a projection point in that case).

Note:

• If d_cart and d are horizontal surfaces (3rd dimension has size==1), it also interpolates the depth coordinate.

d_cart = project_CartData(d_cart::CartData, d::UTMData, p::ProjectionPoint)

Projects all datafields from the UTMData struct d to the CartData struct d_cart, around the projection point p. d_cart must be an orthogonal cartesian grid (deformed doesn't work; use convert2CartData(d, proj), where proj is a projection point in that case).

# Note:    
- If `d_cart` and `d` are horizontal surfaces (3rd dimension has size==1), it also interpolates the depth coordinate.

data_cart = project_FEData_CartData(data_cart::CartData, data_fe::FEData)

Projects a FEData object with tetrahedrons (e.g., from Gmsh) to a Cartesian grid

Grid::LaMEM_grid = read_LaMEM_inputfile(file, args::Union{String,Nothing}=nothing)

Parses a LaMEM input file and stores grid information in the Grid structure. Optionally, you can pass LaMEM command-line arguments as well.

Example 1

julia> Grid = read_LaMEM_inputfile("SaltModels.dat")
LaMEM Grid:
nel         : (32, 32, 32)
marker/cell : (3, 3, 3)
markers     : (96, 96, 96)
x           ϵ [-3.0 : 3.0]
y           ϵ [-2.0 : 2.0]
z           ϵ [-2.0 : 0.0]

Example 2 (with command-line arguments)

julia> Grid = read_LaMEM_inputfile("SaltModels.dat", args="-nel_x 64 -coord_x -4,4")
LaMEM Grid:
  nel         : (64, 32, 32)
  marker/cell : (3, 3, 3)
  markers     : (192, 96, 96)
  x           ϵ [-4.0 : 4.0]
  y           ϵ [-2.0 : 2.0]
  z           ϵ [-2.0 : 0.0]

Data::ParaviewData = read_data_PVTR(fname, dir)

Reads a parallel, rectilinear, *.vts file with the name fname and located in dir and create a 3D Data struct from it.

Example

julia> Data = read_data_PVTR("Haaksbergen.pvtr", "./Timestep_00000005_3.35780500e-01/")
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)

coord, Data_3D_Arrays, Name_Vec = read_data_VTR(fname)

Reads a VTR (structured grid) VTK file fname and extracts the coordinates, data arrays and names of the data. In general, this only contains a piece of the data, and one should open a *.pvtr file to retrieve the full data

This reads the picked profiles from disk and returns a vector of ProfileData

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.

V = removefield(V::AbstractGeneralGrid,field_name::String)

Removes the field with name field_name from the GeoData or CartData dataset

Data_R = rotate_translate_scale(Data::Union{ParaviewData, CartData}; Rotate=0, Translate=(0,0,0), Scale=(1.0,1.0,1.0), Xc=(0.0,0.0))

Does an in-place rotation, translation and scaling of the Cartesian dataset Data.

Parameters

Note that we apply the transformations in exactly this order:

• Scale: scaling applied to the x,y,z coordinates of the data set
• Rotate: rotation around the x/y axis (around the center of the box)
• Translate: translation
• Xc: center of rotation

Example

julia> X,Y,Z   =   xyz_grid(10:20,30:40,-50:-10);
julia> Data_C  =   ParaviewData(X,Y,Z,(Depth=Z,))
ParaviewData
  size  : (11, 11, 41)
  x     ϵ [ 10.0 : 20.0]
  y     ϵ [ 30.0 : 40.0]
  z     ϵ [ -50.0 : -10.0]
  fields: (:Depth,)
julia> Data_R = rotate_translate_scale(Data_C, Rotate=30);
julia> Data_R
ParaviewData
  size  : (11, 11, 41)
  x     ϵ [ 8.169872981077807 : 21.83012701892219]
  y     ϵ [ 28.16987298107781 : 41.83012701892219]
  z     ϵ [ -50.0 : -10.0]
  fields: (:Depth,)

save_GMG(filename::String, data::Union{GeoData, CartDat, UTMData}; dir=pwd())

Saves the dataset data to a JLD2 file (name without extension) in the directory dir

Example

julia> Lon3D,Lat3D,Depth3D = lonlatdepth_grid(1.0:3:10.0, 11.0:4:20.0, (-20:5:-10)*km);
julia> Data_set    =   GeophysicalModelGenerator.GeoData(Lon3D,Lat3D,Depth3D,(DataFieldName=Depth3D,))   
julia> save_GMG("test",Data_set)

save_LaMEM_markers_parallel(Grid::CartData; PartitioningFile=empty, directory="./markers", verbose=true, is64bit=false)

Saves a LaMEM marker file from the CartData structure Grid. It must have a field called Phases, holding phase information (as integers) and optionally a field Temp with temperature info. It is possible to provide a LaMEM partitioning file PartitioningFile. If not, output is assumed to be for one processor. By default it is assumed that the partitioning file was generated on a 32bit PETSc installation. If Int64 was used instead, set the flag.

The size of Grid should be consistent with what is provided in the LaMEM input file. In practice, the size of the mesh can be retrieved from a LaMEM input file using read_LaMEM_inputfile.

Example

julia> Grid    = read_LaMEM_inputfile("LaMEM_input_file.dat")
julia> Phases  = zeros(Int32,size(Grid.X));
julia> Temp    = ones(Float64,size(Grid.X));
julia> Model3D = CartData(Grid, (Phases=Phases,Temp=Temp))
julia> save_LaMEM_markers_parallel(Model3D)
Writing LaMEM marker file -> ./markers/mdb.00000000.dat

If you want to create a LaMEM input file for multiple processors:

julia> save_LaMEM_markers_parallel(Model3D, PartitioningFile="ProcessorPartitioning_4cpu_1.2.2.bin")
Writing LaMEM marker file -> ./markers/mdb.00000000.dat
Writing LaMEM marker file -> ./markers/mdb.00000001.dat
Writing LaMEM marker file -> ./markers/mdb.00000002.dat
Writing LaMEM marker file -> ./markers/mdb.00000003.dat

save_LaMEM_topography(Topo::CartData, filename::String)

This writes a topography file Topo for use in LaMEM, which should have size (nx,ny,1) and contain the field :Topography

Data = screenshot_to_CartData(filename::String, Corner_LowerLeft, Corner_UpperRight; Corner_LowerRight=nothing, Corner_UpperLeft=nothing)

Does the same as screenshot_to_GeoData, but returns a CartData structure

screenshot_to_GeoData(filename::String, Corner_LowerLeft, Corner_UpperRight; Corner_LowerRight=nothing, Corner_UpperLeft=nothing, Cartesian=false, UTM=false, UTMzone, isnorth=true, fieldname::Symbol=:colors)

Take a screenshot of Georeferenced image either a lat/lon, x,y (if Cartesian=true) or in UTM coordinates (if UTM=true) at a given depth or along profile and converts it to a GeoData, CartData or UTMData struct, which can be saved to Paraview

The lower left and upper right coordinates of the image need to be specified in tuples of (lon,lat,depth) or (UTM_ew, UTM_ns, depth), where depth is negative inside the Earth (and in km).

The lower right and upper left corners can be specified optionally (to take non-orthogonal images into account). If they are not specified, the image is considered orthogonal and the corners are computed from the other two.

Note: if your data is in UTM coordinates you also need to provide the UTMzone and whether we are on the northern hemisphere or not (isnorth).

Data = screenshot_to_UTMData(filename::String, Corner_LowerLeft, Corner_UpperRight; Corner_LowerRight=nothing, Corner_UpperLeft=nothing, UTMzone::Int64=nothing, isnorth::Bool=true, fieldname=:colors)

Does the same as screenshot_to_GeoData, but returns for UTM data Note that you have to specify the UTMzone and isnorth

V_sub = subtract_horizontalmean(V::AbstractArray{T, 3}; Percentage=false)

Subtracts the horizontal average of the 3D data array V.

If Percentage=true, the result is given as percentage; otherwise absolute values are returned

V_sub = subtract_horizontalmean(V::AbstractArray{T, 2}; Percentage=false)

Subtracts the horizontal average of the 2D data array V.

If Percentage=true, the result is given as percentage; otherwise absolute values are returned

fe_swap = swap_yz_dims(fe_data::FEData)

This swaps the y and z dimensions of the FEData object, which is useful for pTatin as it uses y for what is z in GMG.

velocity_spherical_to_cartesian!(Data::GeoData, Velocity::Tuple)

In-place conversion of velocities in spherical velocities [Veast, Vnorth, Vup] to cartesian coordinates (for use in paraview).

NOTE: the magnitude of the vector will be the same, but the individual [Veast, Vnorth, Vup] components will not be retained correctly (as a different [x,y,z] coordinate system is used in paraview). Therefore, if you want to display or color that correctly in Paraview, you need to store these magnitudes as separate fields

votemap(DataSets::Vector{GeoData}, criteria::Vector{String}, dims=(50,50,50))

Creates a Vote map which shows consistent features in different 2D/3D tomographic datasets.

The way it works is:

• Find a common region between the different GeoData sets (overlapping lon/lat/depth regions)
• Interpolate the fields of all DataSets to common coordinates
• Filter data points in one model (e.g., areas with a velocity anomaly > 2 percent). Set everything that satisfies this criteria to 1 and everything else to 0.
• Sum the results of the different datasets

If a feature is consistent between different datasets, it will have larger values.

Example

We assume that we have 2 seismic velocity datasets Data_Zhao_Pwave and DataKoulakov_Alps:

julia> Data_Zhao_Pwave
GeoData
  size  : (121, 94, 101)
  lon   ϵ [ 0.0 : 18.0]
  lat   ϵ [ 38.0 : 51.95]
  depth ϵ [ -1001.0 km : -1.0 km]
  fields: (:dVp_Percentage,)
julia> DataKoulakov_Alps
  GeoData
    size  : (108, 81, 35)
    lon   ϵ [ 4.0 : 20.049999999999997]
    lat   ϵ [ 37.035928143712574 : 49.01197604790419]
    depth ϵ [ -700.0 km : -10.0 km]
    fields: (:dVp_percentage, :dVs_percentage)

You can create a votemap which combines the two data sets with:

julia> Data_VoteMap = votemap([Data_Zhao_Pwave,DataKoulakov_Alps],["dVp_Percentage>2.5","dVp_percentage>3.0"])
GeoData
  size  : (50, 50, 50)
  lon   ϵ [ 4.0 : 18.0]
  lat   ϵ [ 38.0 : 49.01197604790419]
  depth ϵ [ -700.0 km : -10.0 km]
  fields: (:votemap,)

You can also create a votemap of a single dataset:

julia> Data_VoteMap = votemap(Data_Zhao_Pwave,"dVp_Percentage>2.5", dims=(50,51,52))
GeoData
  size  : (50, 51, 52)
  lon   ϵ [ 0.0 : 18.0]
  lat   ϵ [ 38.0 : 51.95]
  depth ϵ [ -1001.0 km : -1.0 km]
  fields: (:votemap,)

voxel_grav(X::Array{Float64, 3}, Y::Array{Float64, 3}, Z::Array{Float64, 3}, RHO::Array{Float64, 3};
refMod="AVG", lengthUnit="m", rhoTol=1e-9, Topo=[], outName="Bouguer", printing=true)

Computes Bouguer anomalies and gradients

Required arguments:

• X,Y,Z: 3D matrices with the coordinates of the grid (X should vary in the first dimension, Y in the second, Z (vertical) in the third)
• RHO: 3D matrix with the density at each grid point [kg/m^3]

Optional arguments:

• refMod: 1D vector with the reference density for each depth. Alternatively, the strings "NE", "SE", "SW", "NW", "AVG" can be used. In that case, one of the corners of RHO is used as reference model.In case of "AVG" the reference model is the average of each depth slice.
• lengthUnit: The unit of the coordinates and topography file. Either "m" or "km"
• rhoTol: density differences smaller than rhoTol will be ignored [kg/m^3]
• Topo: 2D matrix with the topography of the surface (only relevant for the paraview output)
• outName: name of the paraview output (do not include file type)
• printing: activate printing of additional information [true or false]

write_FEmesh_fields(data::FEData)

Writes cell and vertex fields to disk to be read by pTatin

write_FEmesh_mesh(vdata::FEData; out_file="md.bin", connectivity_zero_based=true)

Writes a binary file with the mesh information for pTatin

write_pTatin_mesh(fe_mesh::FEData; out_file="md.bin", connectivity_zero_based=true)

Write a binary file with the mesh information for pTatin

write_pTatin_mesh(q1_mesh::Q1Data; out_file="md.bin", connectivity_zero_based=true)

Write a binary file with the mesh information for pTatin

pvd = write_paraview(DataSet::ParaviewData, filename="test"; PointsData=false, pvd=nothing, time=nothing, directory=nothing, verbose=true)

Writes a structure with Geodata to a paraview (or VTK) file. If you have unstructured points (e.g., earthquake data), set PointsData=true. In case you want to create a movie in Paraview, and this is a timestep of that movie you also have to pass time and pvd

Example 1: Write a 3D volume

julia> Lon,Lat,Depth   =   lonlatdepth_grid(10:20,30:40,(-300:25:0)km);
julia> Data_set        =   GeoData(Lat,Lon,Depth,(Depthdata=Depth,LonData=Lon))  
julia> write_paraview(Data_set, "test_depth3D")

Example 2: Horizontal slice @ given depth

julia> Lon,Lat,Depth  =   lonlatdepth_grid(10:20,30:40,10km);
julia> Data_set       =   GeoData(Lat,Lon,Depth,(Topography=Depth,))  
julia> write_paraview(Data_set, "test")

Example 3: Case with topography

julia> Lon,Lat,Depth    =   lonlatdepth_grid(10:20,30:40,10km);
julia> Depth[2:4,2:4,1] .=  25km     
julia> Data_set         =   GeoData(Lat,Lon,Depth,(Topography=Depth,))  
julia> write_paraview(Data_set, "test2")

Example 4: Profile

julia> Lon,Lat,Depth  =   lonlatdepth_grid(10:20,35,(-300:25:0)km);
julia> Data_set       =   GeoData(Lat,Lon,Depth,(DataSet=Depth,Depth=Depth))  
julia> write_paraview(Data_set, "test")

Example 5: Velocity vectors

julia> Lon,Lat,Depth  =   lonlatdepth_grid(10:20,30:40,10km);
julia> Ve, Vn, Vz     =   ones(size(Depth)), ones(size(Depth))*0.5, zeros(size(Depth));
julia> Data_set       =   GeoData(Lat,Lon,Depth,(DataSet=Depth, Velocity=(Ve,Vn,Vz)))
GeoData 
  size  : (11, 11, 1)
  lon   ϵ [ 30.0 - 40.0]
  lat   ϵ [ 10.0 - 20.0]
  depth ϵ [ 10.0 km - 10.0 km]
  fields: (:DataSet, :Velocity)  
julia> write_paraview(Data_set, "test_Velocity")

Example 6: Unconnected points (e.g., earthquake locations)

Note that these points should be 1D vectors.

julia> Lon,Lat,Depth  =   lonlatdepth_grid(10:5:20,35:2:40,(-300:50:0)km);
julia> Lon=Lon[:]; Lat=Lat[:]; Depth=Depth[:];
julia> Data_set       =   GeoData(Lat,Lon,Depth,(DataSet=Depth[:],Depth=Depth*10));  
julia> write_paraview(Data_set, "test_Points", PointsData=true)

write_paraview(DataSet::FEData, filename="test"; directory=nothing, pvd=nothing, time=nothing, verbose=true)

Writes a FEData dataset (general finite element) to disk, which has cell and vertex field

write_paraview(DataSet::Q1Data, filename="test"; directory=nothing, pvd=nothing, time=nothing, verbose=true)

Writes a Q1Data dataset to disk, which has cell and vertex field

write_paraview(DataSet::CartData, filename::Any; PointsData=false, pvd=nothing, time=nothing, directory=nothing, verbose=true)

Writes a CartData structure to paraview.

write_paraview(DataSet::UTMData, filename::Any; PointsData=false, pvd=nothing, time=nothing, directory=nothing, verbose=true)

Writes a UTMData structure to paraview. Note that this data is not transformed into an Earth-like framework, but remains cartesian instead.

X,Y,Z = xyz_grid(X_vec::Any, Y_vec::Any, Z_vec::Any)

Creates a X,Y,Z grid. It works just as lonlatdepth_grid apart from the better suited name.

Example 1: Create 3D grid

julia> X,Y,Z =  xyz_grid(10:20,30:40,(-10:-1)km);
julia> size(X)
(11, 11, 10)

See lonlatdepth_grid for more examples.