List of all functions
Here an overview of all functions:
GeophysicalModelGenerator.BUNDLED_SERVER_INDEX Constant
The copy of the server's index that ships with the package, for when the server cannot be asked for the current one. It changes rarely (new sets are added now and then; the existing lines stay), so a stale copy is far better than none: without it a tiled resolution would be looked for as a single global file, which does not exist.
sourceGeophysicalModelGenerator.GMT_SERVERS Constant
The GMT data server mirrors. The first that answers is used; the others are tried in turn if a download fails, so one server being offline does not fail the import.
sourceGeophysicalModelGenerator.INDEX_TIMEOUT Constant
How long to wait for the server's index file, which is a hundred kilobytes: if that does not arrive within a few seconds the mirror is not going to deliver the tiles either.
sourceGeophysicalModelGenerator.SERVER_RETRY_AFTER Constant
Once every mirror has failed to deliver a file, how long (seconds) to leave the server alone before asking it again. Without this, every later call in the session – and in every other test worker, hence the marker file rather than a variable – waits out the full retry budget too, and a run that should fail in a minute takes hours.
sourceGeophysicalModelGenerator.TILE_TIMEOUT Constant
How long to wait for one tile before giving up on that mirror and trying the next. The transfers are a few megabytes, so a minute is generous; without a limit a throttled or stalled connection blocks indefinitely.
sourceGeophysicalModelGenerator.CartData Method
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.
GeophysicalModelGenerator.CartData Method
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.
GeophysicalModelGenerator.CartData Method
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"]GeophysicalModelGenerator.CartGrid Type
Structure that holds data for an orthogonal cartesian grid, which can be described with 1D vectorsGeophysicalModelGenerator.GMG_Dataset Type
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?
GeophysicalModelGenerator.GeoData Method
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,)GeophysicalModelGenerator.LaMEMPartitioningInfo Type
Structure that holds information about the LaMEM partitioning
sourceGeophysicalModelGenerator.ParaviewData Method
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.
GeophysicalModelGenerator.ParticlesDistribution Type
Structure that holds information about the LaMEM particles distribution for partitioning
sourceGeophysicalModelGenerator.ProfileData Type
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}
TopoData :: Union{Nothing, NamedTuple}
PointData :: Union{Nothing, NamedTuple}
ScreenshotData :: Union{Nothing, NamedTuple}
end
Structure to store cross section dataGeophysicalModelGenerator.ProjectionPoint Method
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
sourceGeophysicalModelGenerator.ProjectionPoint Method
ProjectionPoint(; Lat=49.9929, Lon=8.2473)Defines a projection point used for map projections, by specifying latitude and longitude
sourceGeophysicalModelGenerator.Q1Data Method
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"]GeophysicalModelGenerator.Trench Type
Trench structureStructure that defines the geometry of the trench and the slab.
Parameters
Start- Start of the trench (x,y) coordinatesEnd- End of the trench (x,y) coordinatesn_seg- The number of segment through which the slab is discretize along the dipLength- The length of the slabThickness- The thickness of the slabLb- 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 ofl(∈ [0,Length]).lis the actual distance along the slab length from the trench. In case: -:Linearmath θ(l) = ((θ_max - 0.0)/(Lb-0))*l; -:Ribemath θ(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
GeophysicalModelGenerator.Create1D_grid_vector Method
Returns 1D coordinate vectors of grid points and of marker locations for a regular spacing
sourceGeophysicalModelGenerator.Create1D_grid_vector Method
Returns 1D coordinate vectors of grid points and of marker locations for a regular spacing
sourceGeophysicalModelGenerator.InterpolateDataFields2D_vecs Method
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
GeophysicalModelGenerator.ParseValue_CommandLineArgs Method
This parses a LaMEM command line argument string and checks if the keyword exists thereGeophysicalModelGenerator.ParseValue_LaMEM_InputFile Method
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")GeophysicalModelGenerator.PetscBinaryWrite_Vec Method
PetscBinaryWrite_Vec(filename, A)Writes a vector A to disk, such that it can be read with PetscBinaryRead (which assumes a Big Endian type)
GeophysicalModelGenerator.Rot3D Method
xrot, yrot, zrot = Rot3D(X::Number,Y::Number,Z::Number, cosStrikeAngle, sindStrikeAngle, cosDipAngle, sinDipAngle)
Perform rotation for a point in 3D spaceGeophysicalModelGenerator.above_surface Method
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
GeophysicalModelGenerator.above_surface Method
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.
GeophysicalModelGenerator.above_surface Method
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
GeophysicalModelGenerator.above_surface Method
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
GeophysicalModelGenerator.above_surface Method
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
GeophysicalModelGenerator.add_box! Method
add_box!(Phase, Temp, Grid::AbstractGeneralGrid,
bounds::AbstractVector{<:AbstractVector{<:Real}};
Origin = nothing, StrikeAngle = 0, DipAngle = 0,
phase = ConstantPhase(1),
T = nothing,
segments = nothing,
cell = false )Add box function but getting bounds as a vector bounds in a way of [[xmin,xmax],[ymin,ymax],[zmin,zmax]]. If bounds is empty return nothing
sourceGeophysicalModelGenerator.add_box! Method
add_box!(Phase, Temp, Grid::AbstractGeneralGrid; xlim::Tuple = (20,100), [ylim::Tuple = (1,10)], zlim::Tuple = (10,80),
Origin=nothing, StrikeAngle=0, DipAngle=0,
phase = ConstantPhase(1),
T=nothing,
segments=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 inGMG)xlim- left/right coordinates of boxylim- front/back coordinates of box [optional; if not specified we use the whole box]zlim- bottom/top coordinates of boxOrigin- the origin, used to rotate the box around. Default is the left-front-top cornerStrikeAngle- strike angle of slabDipAngle- dip angle of slabphase- specifies the phase of the box. SeeConstantPhase(),LithosphericPhases()T- specifies the temperature of the box. SeeConstantTemp(),LinearTemp(),HalfspaceCoolingTemp(),SpreadingRateTemp(),LithosphericTemp()segments- optional parameter to define multiple ridge segments within the boxcell- if true,PhaseandTempare 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"Example 3) Box with ridge thermal structure ```julia-repl julia> Grid = CartData(xyz_grid(-1000:10:1000, -1000:10:1000, -660:5:0)) julia> Phases = fill(2, size(Grid)); julia> Temp = fill(1350.0, size(Grid)); julia> segments = [((-500.0, -1000.0), (-500.0, 0.0)), ((-250.0, 0.0), (-250.0, 200.0)), ((-750.0, 200.0), (-750.0, 1000.0))]; julia> lith = LithosphericPhases(Layers=[15 55], Phases=[1 2], Tlab=1250); julia> add_box!(Phases, Temp, Grid; xlim=(-1000.0, 0.0), ylim=(-500.0, 500.0), zlim=(-80.0, 0.0), phase=lith, T=SpreadingRateTemp(SpreadingVel=3), segments=segments) julia> Grid = addfield(Grid, (; Phases, Temp)); # Add to Cartesian model julia> write_paraview(Grid, "Ridge_Thermal_Structure") # Save model to Paraview 1-element Vector{String}: "Ridge_Thermal_Structure.vts"
sourceGeophysicalModelGenerator.add_cylinder! Method
add_cylinder!(Phase, Temp, Grid::AbstractGeneralGrid; base::Tuple = (-1,-1,-1.5), cap::Tuple = (-1,-1,-0.5), radius::Number,
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 withread_LaMEM_inputfile)base- center coordinate of bottom of cylindercap- center coordinate of top of cylinderradius- radius of the cylinderphase- specifies the phase of the box. SeeConstantPhase(),LithosphericPhases()T- specifies the temperature of the box. SeeConstantTemp(),LinearTemp(),HalfspaceCoolingTemp(),SpreadingRateTemp()cell- if true,PhaseandTempare 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"GeophysicalModelGenerator.add_ellipsoid! Method
add_ellipsoid!(Phase, Temp, Grid::AbstractGeneralGrid; cen::Tuple = (-1,-1,-1), axes::Tuple = (0.2,0.1,0.5),
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 read_LaMEM_inputfile)cen- center coordinates of sphereaxes- semi-axes of ellipsoid in X,Y,ZOrigin- the origin, used to rotate the box around. Default is the left-front-top cornerStrikeAngle- strike angle of slabDipAngle- dip angle of slabphase- specifies the phase of the box. SeeConstantPhase(),LithosphericPhases()T- specifies the temperature of the box. SeeConstantTemp(),LinearTemp(),HalfspaceCoolingTemp(),SpreadingRateTemp()cell- if true,PhaseandTempare 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"GeophysicalModelGenerator.add_fault! Method
add_fault!(Phase, Temp, Grid::AbstractGeneralGrid;
Start=(20,100), End=(10,80),
Fault_thickness=10.0,
Depth_extent=nothing,
DipAngle=0e0,
phase=ConstantPhase(1),
T=nothing,
cell=false)Adds a fault to the given 3D grid by modifying the Phase and Temp arrays. For a 2D grid, use add_box instead.
Arguments
Phase: Phase arrayTemp: Temp arrayGrid: The grid on which the fault is to be added.Start: Tuple representing the starting coordinates of the fault (X, Y).End: Tuple representing the ending coordinates of the fault (X, Y).Fault_thickness: Thickness of the fault.Depth_extent: Depth extent of the fault. Ifnothing, the fault extends through the entire domain.DipAngle: Dip angle of the fault.phase: Phase to be assigned to the fault.T: Temperature to be assigned to the fault. Ifnothing, the temperature is not modified.
Example
add_fault!(Phase, Temp, Grid;
Start=(20,100), End=(10,80),
Fault_thickness=10.0,
Depth_extent=(-25.0, 0.0),
DipAngle=-10.0,
phase=ConstantPhase(1)
)GeophysicalModelGenerator.add_layer! Method
add_layer!(Phase, Temp, Grid::AbstractGeneralGrid; xlim::Tuple = (1,100), [ylim::Tuple = (0,20)], zlim::Tuple = (0,-100),
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 read_LaMEM_inputfile, but can also be other grid types)xlim- left/right coordinates of boxylim- front/back coordinates of boxzlim- bottom/top coordinates of boxphase- specifies the phase of the box. SeeConstantPhase(),LithosphericPhases()T- specifies the temperature of the box. SeeConstantTemp(),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"GeophysicalModelGenerator.add_plate! Method
add_plate!(Phase, Temp, Grid::AbstractGeneralGrid,
bounds::AbstractVector{<:AbstractVector{<:Real}};
phase = ConstantPhase(1),
T = nothing,
segments = nothing,
cell = false )Add plate function but getting bounds as a vector bounds in a way of [[xmin,xmax],[ymin,ymax],[zmin,zmax]]. If bounds is empty return nothing
sourceGeophysicalModelGenerator.add_polygon! Method
add_polygon!(Phase, Temp, Grid::AbstractGeneralGrid,
bounds::AbstractVector{<:AbstractVector{<:Real}};
phase = ConstantPhase(1),
T = nothing,
cell = false )Add polygon function but getting bounds as a vector bounds in a way of [[xmin,xmax],[ymin,ymax],[zmin,zmax]]. If bounds is empty return nothing
sourceGeophysicalModelGenerator.add_polygon! Method
add_polygon!(Phase, Temp, Grid::AbstractGeneralGrid; xlim=(), ylim::Tuple = (0.0,0.8), zlim=(), 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 read_LaMEM_inputfile)xlim-x-coordinate of the polygon points, same ordering as zlim, number of points unlimitedylim-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 unlimitedphase- specifies the phase of the box. SeeConstantPhase()T- specifies the temperature of the box. SeeConstantTemp(),LinearTemp(),HalfspaceCoolingTemp(),SpreadingRateTemp()cell- if true,PhaseandTempare 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, 1.6, 2.0),ylim=(0,0.8), zlim=(0,-1,-2,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"GeophysicalModelGenerator.add_slab! Method
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 inGMG)trench- Trench structurephase- specifies the phase of the box. SeeConstantPhase(),LithosphericPhases()T- specifies the temperature of the box. SeeConstantTemp(),LinearTemp(),HalfspaceCoolingTemp(),SpreadingRateTemp(),LithosphericTemp()cell- if true,PhaseandTempare 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)GeophysicalModelGenerator.add_sphere! Method
add_sphere!(Phase, Temp, Grid::AbstractGeneralGrid; cen::Tuple = (0,0,-1), radius::Number,
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 read_LaMEM_inputfile)cen- center coordinates of sphereradius- radius of spherephase- specifies the phase of the box. SeeConstantPhase(),LithosphericPhases()T- specifies the temperature of the box. SeeConstantTemp(),LinearTemp(),HalfspaceCoolingTemp(),SpreadingRateTemp()cell- if true,PhaseandTempare 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"GeophysicalModelGenerator.add_stripes! Method
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 read_LaMEM_inputfile, 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 ZstripeWidth- width of the stripestripeSpacing- space between two stripesOrigin- the origin, used to rotate the box around. Default is the left-front-top cornerStrikeAngle- strike angleDipAngle- dip anglephase- specifies the phase we want to apply stripes tostripePhase- specifies the stripe phasecell- if true,PhaseandTempare 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"GeophysicalModelGenerator.add_volcano! Method
add_volcano!( Phases, Temp, Grid::CartData; volcanic_phase, center, height, radius, crater, base, background, T, )
Adds a volcano topography (cones and truncated cones)
Parameters
Phases - Phase array (consistent with Grid)
Temp - Temperature array (consistent with Grid)
Grid - CartData
Optional Parameters
volcanic_phase - phase number of the volcano,
center - x- and -coordinates of center of volcano
height - height of volcano
radius - radius of volcano
T - temperature structure of the volcano
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)
GeophysicalModelGenerator.addfield Method
V = addfield(V::AbstractGeneralGrid,field_name::String,data::Any)Add Fields Data to GeoData or CartData
sourceGeophysicalModelGenerator.addfield Method
V = addfield(V::CartData,new_fields::NamedTuple)Add new_fields fields to a CartData dataset
GeophysicalModelGenerator.addfield Method
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
GeophysicalModelGenerator.addfield Method
V = addfield(V::GeoData,new_fields::NamedTuple)Add new_fields fields to a GeoData dataset
GeophysicalModelGenerator.addfield Method
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
GeophysicalModelGenerator.average_q1 Method
out = average_q1(d::Array)3D linear averaging of a 3D array
sourceGeophysicalModelGenerator.below_surface Function
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
GeophysicalModelGenerator.below_surface Method
Below = below_surface(Grid::CartGrid, DataSurface_Cart::CartData)
Determines if points described by the `Grid` CartGrid structure are above the Cartesian surface `DataSurface_Cart`GeophysicalModelGenerator.below_surface Method
Below = below_surface(Data::GeoData, DataSurface::GeoData)Determines if points within the 3D Data structure are below the GeoData surface DataSurface
GeophysicalModelGenerator.below_surface Method
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
GeophysicalModelGenerator.below_surface Method
Below = below_surface(Data_Cart::ParaviewData, DataSurface_Cart::ParaviewData)Determines if points within the 3D Data_Cart structure are below the Cartesian surface DataSurface_Cart
sourceGeophysicalModelGenerator.bilinear Method
bilinear(xs, ys, A, x, y)Sample A at (x, y) by interpolating between the four surrounding points, which is what GMT does when it resamples a coarser set into a finer grid. Falls back to the nearest edge value outside the grid.
GeophysicalModelGenerator.cell_area Method
area_m2 = cell_area(Topo::GeoData)Returns the cell area for a Topographic dataset in m² (required for upstream area calculation)
sourceGeophysicalModelGenerator.check_multigrid_compatibility Function
check_multigrid_compatibility(nx::Int, ny::Int, nz::Int, total_ranks::Int)Check if the model resolution and MPI ranks are optimal for multigrid methods.
sourceGeophysicalModelGenerator.combine_vol_data Method
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.
GeophysicalModelGenerator.compute_bending_angle Method
θ = 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]
GeophysicalModelGenerator.compute_phase Method
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 read_LaMEM_inputfile)
GeophysicalModelGenerator.compute_slab_surface Method
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-Trenchstructure 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.
GeophysicalModelGenerator.compute_thermal_structure Method
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) -> the weight of F1 increase from the origin, while F2 decreases.
This function has been conceived for averaging the solution of McKenzie and half space cooling models, but it 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
then modify F.
GeophysicalModelGenerator.compute_thermal_structure Method
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 arrayX: X ArrayY: Y ArrayZ: Z ArrayPhase: Phase arrays:McKenzie_subducting_slab
GeophysicalModelGenerator.convert2CartData Method
convert2CartData(d::GeoData, proj::ProjectionPoint)Converts a GeoData structure to a CartData structure, which essentially transfers the dimensions to km
GeophysicalModelGenerator.convert2CartData Method
convert2CartData(d::UTMData, proj::ProjectionPoint)Converts a UTMData structure to a CartData structure, which essentially transfers the dimensions to km
GeophysicalModelGenerator.convert2FEData Method
fe_data::FEData = convert2FEData(d::Q1Data)Creates a Q1 FEM mesh from the Q1Data data which holds the vertex coordinates and cell/vertex fields
GeophysicalModelGenerator.convert2UTMzone Method
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
GeophysicalModelGenerator.convert2UTMzone Method
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.
GeophysicalModelGenerator.coordinate_grids Method
X,Y,Z = coordinate_grids(Data::CartData; cell=false)Returns 3D coordinate arrays
sourceGeophysicalModelGenerator.coordinate_grids Method
X,Y,Z = coordinate_grids(Data::CartGrid; cell=false)Returns 3D coordinate arrays
sourceGeophysicalModelGenerator.coordinate_grids Method
LON,LAT,Z = coordinate_grids(Data::GeoData; cell=false)Returns 3D coordinate arrays
sourceGeophysicalModelGenerator.coordinate_grids Method
X,Y,Z = coordinate_grids(Data::LaMEM_grid; cell=false)Returns 3D coordinate arrays
sourceGeophysicalModelGenerator.coordinate_grids Method
X,Y,Z = coordinate_grids(Data::ParaviewData; cell=false)Returns 3D coordinate arrays
sourceGeophysicalModelGenerator.coordinate_grids Method
X,Y,Z = coordinate_grids(Data::Q1Data; cell=false)Returns 3D coordinate arrays
sourceGeophysicalModelGenerator.coordinate_grids Method
EW,NS,Z = coordinate_grids(Data::UTMData; cell=false)Returns 3D coordinate arrays
sourceGeophysicalModelGenerator.countmap Method
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
sourceGeophysicalModelGenerator.create_CartGrid Method
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)GeophysicalModelGenerator.create_partitioning_file Method
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).
GeophysicalModelGenerator.create_profile_volume! Method
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.
GeophysicalModelGenerator.create_profile_volume! Method
create_profile_volume!(Profile::ProfileData, VolData; DimsVolCross::NTuple=(100,100), Depth_extent=nothing)Fallback that normalizes VolData: nothing and empty containers (NamedTuple(), ()) are treated as "no volume data" and leave Profile unchanged. Any other unsupported type throws an ArgumentError.
GeophysicalModelGenerator.create_profile_volume! Method
create_profile_volume!(Profile::ProfileData, VolData::NamedTuple; 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. This function allows to pass the data as NamedTuples instead of a GeoData object.
GeophysicalModelGenerator.crop_bounds Method
crop_bounds(uncropped_bounds, proc_bounds, x, y, z)
Crop boundaries from the whole model to only the extent of the current processor
uncropped_bounds: 3 element vector with maximum and minimum values[[x_min, x_max],[y_min, y_max],[z_min, z_max]] for the whole model
proc_bounds: 3 element vector with maximum and minimum values[[x_min, x_max],[y_min, y_max],[z_min, z_max]] for the current processor
x,y,z: 3 element vector with maximum and minimum values[[x_min, x_max],[y_min, y_max],[z_min, z_max]] for the current processor
Returns a 3 element vector with maximum and minimum values[[x_min, x_max],[y_min, y_max],[z_min, z_max]] for the current processorGeophysicalModelGenerator.cross_section Method
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_levelis specifiedThey can also be vertical, either by specifying
Lon_levelorLat_level(for a fixed lon/lat), or by defining bothStart=(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:
Volume data: data will be interpolated or directly extracted from the data set.
Surface data: surface data will be interpolated or directly extracted from the data set
Point data: data will be projected to the chosen profile. Only data within a chosen distance (default is 50 km) will be used
Interpolateindicates 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 indimsNOTE: 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)GeophysicalModelGenerator.cross_section_points Method
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
sourceGeophysicalModelGenerator.cross_section_surface Method
cross_section_surface(Surface::GeoData; dims=(100,), Interpolate=false, Depth_level=nothing; Lat_level=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_levelis specifiedThey can also be vertical, either by specifying
Lon_levelorLat_level(for a fixed lon/lat), or by defining bothStart=(lon,lat)&End=(lon,lat)points. Start and End points will be in km!IMPORTANT: The surface to be extracted has to be given as a gridded GeoData object. It may also contain NaNs where it is not defined. Any points lying outside of the defined surface will be considered NaN.
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"]GeophysicalModelGenerator.cross_section_volume Method
cross_section_volume(Volume::AbstractGeneralGrid; dims=(100,100), Interpolate=false, Depth_level=nothing; Lat_level=nothing; Lon_level=nothing; Start=nothing, End=nothing, Depth_extent=nothing )
Creates a cross-section through a volumetric (3D) GeoData object.
Cross-sections can be horizontal (map view at a given depth), if
Depth_levelis specifiedThey can also be vertical, either by specifying
Lon_levelorLat_level(for a fixed lon/lat), or by defining bothStart=(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 Depth_extent=(depth_min,depth_max)Interpolateindicates 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 indimsDepth_extentis 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)GeophysicalModelGenerator.dataset_scale Method
dataset_scale(dataset, res, reg)The scale factor and offset that turn the integers in a tile into metres. The server publishes these in gmt_data_server.txt, one line per set: earth_relief_15s_p is stored in half metres, for instance, while the SRTM sets are stored directly in metres. Reading them rather than assuming keeps this correct if a set is ever re-encoded.
Falls back to (1, 0) if the index cannot be fetched, which is right for every set that is not scaled.
GeophysicalModelGenerator.decode_jp2 Method
decode_jp2(file, n)Decode a JPEG2000 tile into an n x n matrix of Int16 elevations, using the decoder that comes with OpenJpeg_jll. The file stores its rows from north to south.
GeophysicalModelGenerator.decompose_mpi_ranks Method
decompose_mpi_ranks(total_ranks::Int, nx::Int, ny::Int, nz::Int) -> Tuple{Int,Int,Int}Decompose total number of MPI ranks into a 3D processor grid (px, py, pz), optimizing for cell aspect ratio closest to 1.0.
sourceGeophysicalModelGenerator.distance_to_linesegment Method
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.
GeophysicalModelGenerator.download_data Function
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"GeophysicalModelGenerator.download_tile Method
download_tile(url, dest)Fetch one file, trying each mirror in turn. Returns true if it was downloaded, false if the file is not on the server – which is not an error: the 3 arcsecond set covers land only, so a tile that is entirely ocean simply does not exist.
GeophysicalModelGenerator.drape_on_topo Method
drape_on_topo(Topo::CartData, Data::CartData)Drapes Cartesian Data on topography
sourceGeophysicalModelGenerator.drape_on_topo Method
Topo = drape_on_topo(Topo::GeoData, Data::GeoData)This drapes fields of a data set Data on the topography Topo
GeophysicalModelGenerator.extract_ProfileData! Function
extract_ProfileData!(Profile::ProfileData, VolData=nothing, SurfData=nothing, PointData=nothing; DimsVolCross=(100,100), Depth_extent=nothing, DimsSurfCross=(100,), section_width=50km, ScreenshotData=nothing, TopoData=nothing)Extracts data along a vertical or horizontal profile. VolData can be a GeoData/AbstractGeneralGrid or a NamedTuple of those; the other data sets are NamedTuples. Any data argument can be nothing or empty (NamedTuple() or ()), in which case it is skipped.
GeophysicalModelGenerator.extract_ProfileData Method
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
sourceGeophysicalModelGenerator.extract_subvolume Method
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_levelandDepth_levelshould be tuples that indicate(minimum_value, maximum_value)along the respective direction. If not specified we use the full range.By default,
Interpolate=falseand 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=truewe interpolate the data onto a new grid that has dimensionsdims. 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 Lon_level/Lat_level/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}GeophysicalModelGenerator.extract_subvolume Method
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_levelandDepth_levelshould be tuples that indicate(minimum_value, maximum_value)along the respective direction. If not specified we use the full range.By default,
Interpolate=falseand 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=truewe interpolate the data onto a new grid that has dimensionsdims. 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 Lon_level/Lat_level/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)GeophysicalModelGenerator.fill_gaps! Method
fill_gaps!(Z, lon, lat, filler)Replace the cells of Z that carry no data – the flat zeros that SRTM writes for the sea, and whatever a missing tile left behind – with values interpolated from the coarser filler set that the server names for the purpose.
GeophysicalModelGenerator.filler_dataset Method
filler_dataset(dataset, res, reg)Some sets cover only part of the globe and are published with a coarser set to fill the rest: the 3 arcsecond relief is SRTM, which is land only, and the server names earth_relief_15s_p as what belongs in the sea. GMT resamples that filler into the gaps, and so do we – otherwise the ocean comes out flat at zero.
Returns nothing when a set needs no filler.
GeophysicalModelGenerator.find_slab_distance! Method
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.
GeophysicalModelGenerator.fit_surface_to_points Method
surf_new = fit_surface_to_points(surf::CartData, X_pt::Vector, Y_pt::Vector, Z_pt::Vector)This fits the z values of the surface surf to the z value of the closest-by-points in (X_pt, Y_pt, Z_pt)
GeophysicalModelGenerator.fit_surface_to_points Method
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)
GeophysicalModelGenerator.flatten_cross_section Method
flatten_cross_section(V::CartData)Takes a diagonal 3D cross_section and flattens it to be converted to a 2D Grid by create_CartGrid
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"]GeophysicalModelGenerator.flatten_cross_section Method
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)
```GeophysicalModelGenerator.flatten_index_dimensions Method
ind2D = flatten_index_dimensions(Phase, ind_vec::Vector{CartesianIndex{3}})This converts the indices to purely 2D indices if the array phase is 2D
GeophysicalModelGenerator.flatten_index_dimensions Method
ind2D = flatten_index_dimensions(Phase, ind_vec::Vector{CartesianIndex{3}})This converts the indices to purely 2D indices if the array phase is 2D
GeophysicalModelGenerator.flip Function
Data = flip(Data::GeoData, dimension=3)This flips the data in the structure in a certain dimension (default is z [3])
sourceGeophysicalModelGenerator.generate_valid_rank_sequence Function
generate_valid_rank_sequence(max_value::Int = 196608, multipliers::Vector{Int} = [3])Generate sequence of valid MPI ranks up to max_value. Sequence includes powers of 2 and their multiples by specified multipliers.
sourceGeophysicalModelGenerator.get_particles_distribution Method
p_dist = get_particles_distribution(Grid, P)Get the distribution of particles in the grid Grid: LaMEM_grid P: LaMEMPartitioningInfo Returns a LaMEMPartitioningInfo object with the distribution of particles in the grid
sourceGeophysicalModelGenerator.get_proc_bound Method
proc_bounds = get_proc_bound(Grid, p_dist, proc_num) Get the bounds of the current processor in x, y, z direction Grid: LaMEM_grid p_dist: LaMEMPartitioningInfo proc_num: processor number Returns a 3 element vector with maximum and minimum values[[x_min, x_max],[y_min, y_max],[z_min, z_max]] for current processor proc_num
Example for a model with 8 MPI ranks for current processor number 2
And gives us coordinates for the current processor
Grid_example = LaMEM_grid(
3, 3, 3,
12, 12, 12,
4, 4, 4,
10, 5, 5,
[10.0, 18.0], [20,28], [30,38],
[],[],[],
[],[],[],
[],[],[],
[],[],[]
)
P_example = setup_model_domain(Grid_example.coord_x, Grid_example.coord_y, Grid_example.coord_z, Grid_example.nel_x, Grid_example.nel_x, Grid_example.nel_x, 8)
p_dist_example = get_ParticlesDistribution(Grid_example,P_example)
proc_bounds = get_proc_bound(Grid_example,p_dist_example,2)GeophysicalModelGenerator.get_proc_grid Method
Grid = get_proc_grid(Grid_info, p_dist, proc_bounds, proc_num, RandomNoise)Get the local grid for the current processor Grid_info: LaMEM_grid p_dist: LaMEMPartitioningInfo proc_bounds: bounds of the current processor proc_num: processor number RandomNoise: add random noise to the grid (0/1) Returns a LaMEM_grid object with the local grid for the current processor
sourceGeophysicalModelGenerator.get_processor_partitioning Method
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.
sourceGeophysicalModelGenerator.getlonlatdepthmag_QuakeML Method
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.
sourceGeophysicalModelGenerator.grid_from_single_file Method
grid_from_single_file(limits, dataset, res, reg)The coarse resolutions are published as one NetCDF file for the whole globe, which is small enough to download and cut down here. Reading it needs NCDatasets, which is a weak dependency – the tiled resolutions, which are the interesting ones, need nothing extra.
GeophysicalModelGenerator.grid_from_tiles Method
grid_from_tiles(limits, dataset, res, reg)Assemble the tiles that cover limits and cut out the requested region. Tiles that are not on the server are taken to be at sea level, which is what the land-only 3 arcsecond set needs.
GeophysicalModelGenerator.import_topo Method
import_topo(limits; dataset="earth_relief", res="01m", reg="p")Download topography (and bathymetry) for the region limits = [lonmin, lonmax, latmin, latmax] from the GMT data server, and return it as GeoData whose Topography field is in km.
This does not require GMT: the grids are fetched over HTTP and decoded here.
Datasets
"earth_relief"(the default) blends land topography with ocean bathymetry"earth_gebco"is GEBCO bathymetry"earth_synbath"is GEBCO with synthetic bathymetry where the sea floor is unsurveyed
Resolutions
"01d", "30m", "20m", "15m", "10m", "06m" come as one file for the globe; "05m" through "01m", "30s", "15s" and "03s" are tiled.
Note that "03s" is SRTM, which covers land only – the sea is flat zero there, and the finest resolution that carries real bathymetry is "15s".
Example
julia> Topo = import_topo([-18.7, -17.1, 28.0, 29.2], res="03s")
julia> Topo = import_topo(lon=[4,20], lat=[37,49])Tiles are cached, so importing the same region again does not download them a second time.
sourceGeophysicalModelGenerator.inpoly Method
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!().
sourceGeophysicalModelGenerator.inpoly_fast Method
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.
sourceGeophysicalModelGenerator.inpolygon! Method
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.
GeophysicalModelGenerator.inpolygon! Method
inpolygon!(inside::Vector, PolyX::Vector, PolyY::Vector, x::Vector, y::Vector; fast=false)Same as above but inside, X and Y and are vectors.
GeophysicalModelGenerator.insert_cpu_lines_after_nelz Method
insert_cpu_lines_after_nelz(filename::String, P) Inserts or updates CPU lines in a LaMEM input file after the nel_z line. If CPU lines already exist, they are replaced with the new values from P. If they do not exist, new lines are inserted. Usage:
- Generate LaMEM model and save it to a setup file, and update to correct target resolution after saving.
Target resolution is 256x256x256 and we want to run on 128 processors n_ranks
nx = 256; ny = 256; nz = 256; n_ranks = 128 ParamFile_name = "LaMEM_setup.txt" model = Model( LaMEM.Grid(x=[-80.,80.], y=[-95.,100.], z=[-60.0,10.0] , nel=(8,8,8)),Output(param_file_name = ParamFile_name,)) write_LaMEM_inputFile(model,ParamFile_name). update_nel_xyz_in_file(ParamFile_name, nx, ny, nz)
- Create a partitioning info object
Pwith the number of processes.
Grid_info = get_LaMEM_grid_info(ParamFile_name) P = setup_model_domain(Grid_info.coord_x, Grid_info.coord_y, Grid_info.coord_z, nx, ny, nz, n_ranks)
- Then we can use this function to insert or update CPU lines in the setup file.
insert_cpu_lines_after_nelz(ParamFile_name, P)
Use-case:This allows to run correctly LaMEM simulation with the correct number of processes on each direction if partitioning file is generated by julia routine: part_filename = write_processor_partitioning_LaMEM(P; is64bit=true) # always use is64bit=true insert_cpu_lines_after_nelz(ParamFile_name, P) Grid_o = read_LaMEM_inputfile(ParamFile_name); Phase_o = zeros(UInt8, size(Grid_o.X)); Temp_o = zeros(size(Phase_o)); Model3D_o = CartData(Grid_o, (Phases=Phase_o,Temp=Temp_o)) save_LaMEM_markers_parallel(Model3D_o, PartitioningFile=part_filename,verbose = true, is64bit = true) run_lamem(ParamFile_name,n_ranks)
sourceGeophysicalModelGenerator.interpolate_data_fields_cross_section Method
interpolate_data_fields_cross_section(V::CartData, X,Y,Z,Xcross)Interpolates data fields along a cross-section defined by Xcross and Z
GeophysicalModelGenerator.interpolate_data_surface Method
Surf_interp = interpolate_data_surface(V::GeoData, Surf::GeoData)Interpolates a 3D data set V on a surface defined by Surf
GeophysicalModelGenerator.interpolate_data_surface Method
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)GeophysicalModelGenerator.interpolate_datafields Method
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)GeophysicalModelGenerator.interpolate_datafields Method
interpolate_datafields(V::UTMData, EW, NS, Depth)Interpolates a data field V on a grid defined by UTM,Depth
GeophysicalModelGenerator.interpolate_datafields_2D Method
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
GeophysicalModelGenerator.interpolate_datafields_2D Method
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
GeophysicalModelGenerator.interpolate_datafields_2D Method
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.
GeophysicalModelGenerator.interpolate_datafields_2D Method
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
GeophysicalModelGenerator.interpolate_datafields_2D Method
interpolate_datafields_2D(Original::GeoData, New::GeoData; Rotate=0.0, Translate=(0,0,0), Scale=(1.0,1.0,1.0))Interpolates a data field Original on a 2D GeoData 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
GeophysicalModelGenerator.interpolate_datafields_2D Method
interpolate_datafields_2D(V::UTMData, EW, NS)Interpolates a data field V on a 2D grid defined by UTM. Typically used for horizontal surfaces
GeophysicalModelGenerator.is_surface Method
issurf = is_surface(surf::AbstractGeneralGrid)Returns true if surf is a horizontal 3D surface.
GeophysicalModelGenerator.isinside_closed_STL Function
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.
sourceGeophysicalModelGenerator.lithostatic_pressure! Method
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.
GeophysicalModelGenerator.load_GMG Function
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"]GeophysicalModelGenerator.load_GMG Method
data::NamedTuple = load_GMG(data::GMG_Dataset)Loads a dataset specified in GMG_Dataset data and returns it as a named tuple
GeophysicalModelGenerator.load_GMG Method
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
GeophysicalModelGenerator.load_dataset_file Method
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 loadedLocation: 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 totrue
GeophysicalModelGenerator.lonlatdepth_grid Method
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, )GeophysicalModelGenerator.make_paraview_collection Method
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*.vtkare stored.pvd_name: filename of the resulting*.pvdfile without extension; if not specified,full_simulationis used.files: Vector of the*.vtkfiles without extension; if not specified, all*.vtkfiles in the directory are used.file_extension: file extension of the vtk files. Default is.vtsbut allvt*work.time: Vector of the timesteps; if not specified, pseudo time steps are assigned.
GeophysicalModelGenerator.make_volc_topo Method
make_volc_topo(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 read_LaMEM_inputfile)
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.vtsGeophysicalModelGenerator.meshgrid Method
meshgrid(vx,vy,vz)Computes an (x,y,z)-grid from the vectors (vx,vy,vz). For more information, see the MATLAB documentation.
sourceGeophysicalModelGenerator.movie_from_images Method
movie_from_images(; dir=pwd(), file=nothing, outfile=nothing, framerate=10, copy_to_current_dir=true, type=:mp4_default, gif=false, 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 (default) or a *.gif animation.
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,fileis used.framerate: number of frames/second.copy_to_current_dir: copies the final movie to the current directory iftrue(default); otherwise it will be stored indir.type: type of movie that is created; possible options are::mp4_default: Default option that saves a well-compressedmp4movie 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):gif: Animated GIF (same as passinggif=true).
gif: iftrue, creates a*.gifanimation in addition to themp4movie. Iftype=:gif, only a gif is created.collect: suppresses output ofFFMPEGiftrue(default).
GeophysicalModelGenerator.movie_paraview Method
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)GeophysicalModelGenerator.native_registration Method
native_registration(dataset, res, reg)Which registration this set is actually published in. Most come in both, and then the caller's reg is used; the SRTM-derived sets are gridline only and 15s is pixel only, so asking for the other one would 404 on every tile.
GeophysicalModelGenerator.nearest_point_indices Method
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
GeophysicalModelGenerator.nearest_point_indices Method
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
GeophysicalModelGenerator.nearest_point_indices Method
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
GeophysicalModelGenerator.parse_columns_CSV Method
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)GeophysicalModelGenerator.parse_topo_file Method
parse_topo_file(file)Take apart the "@earth_relief_01m" spelling into the dataset, resolution and registration that name it on the server.
GeophysicalModelGenerator.point_in_tetrahedron Method
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
GeophysicalModelGenerator.point_to_nearest_grid Method
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)
GeophysicalModelGenerator.point_to_nearest_grid Method
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
GeophysicalModelGenerator.point_to_nearest_grid Method
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
GeophysicalModelGenerator.point_to_nearest_grid Method
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
GeophysicalModelGenerator.point_to_nearest_grid Method
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
GeophysicalModelGenerator.project_CartData Method
d_cart = project_CartData(d_cart::CartData, d_cart_data0::CartData)Projects all datafields from the CartData struct d_cart_data0 to the CartData struct d_cart. d_cart_data0 must be an orthogonal cartesian grid.
Note:
- If
d_cartis a horizontal surface (3rd dimension has size==1), the depth coordinate ofd_cart_data0is interpolated as well.
GeophysicalModelGenerator.project_CartData Method
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_cartanddare horizontal surfaces (3rd dimension has size==1), it also interpolates the depth coordinate.
GeophysicalModelGenerator.project_CartData Method
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.GeophysicalModelGenerator.project_FEData_CartData Method
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
sourceGeophysicalModelGenerator.read_LaMEM_inputfile Method
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]GeophysicalModelGenerator.read_data_PVTR Method
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)GeophysicalModelGenerator.read_data_VTR Method
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
GeophysicalModelGenerator.read_picked_profiles Method
This reads the picked profiles from disk and returns a vector of ProfileData
sourceGeophysicalModelGenerator.read_topo_netcdf Function
read_topo_netcdf(file)Read a global .grd grid, returning (lon, lat, z). Filled in by the NCDatasets extension; the tiled resolutions do not go through here.
GeophysicalModelGenerator.remove_NaN_surface! Method
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.
GeophysicalModelGenerator.removefield Method
V = removefield(V::AbstractGeneralGrid,field_name::String)Removes the field with name field_name from the GeoData or CartData dataset
GeophysicalModelGenerator.removefield Method
V = removefield(V::AbstractGeneralGrid,field_name::Symbol)Removes the field with name field_name from the GeoData or CartData dataset
GeophysicalModelGenerator.removefield Method
V = removefield(V::AbstractGeneralGrid,field_name::NTuple{N,Symbol})Removes the fields in the tuple field_name from the GeoData or CartData dataset
GeophysicalModelGenerator.reset_topo_server Method
reset_topo_server()Forget that the GMT data server was unreachable, so that the next import_topo asks it again right away rather than after SERVER_RETRY_AFTER seconds.
GeophysicalModelGenerator.rotate_translate_scale Method
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 thex,y,zcoordinates of the data setRotate: rotation around thex/yaxis (around the center of the box)Translate: translationXc: 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,)GeophysicalModelGenerator.save_GMG Method
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)GeophysicalModelGenerator.save_LaMEM_markers_parallel Method
save_LaMEM_markers_parallel(Grid::CartData; PartitioningFile=empty, directory="./markers", verbose=true, is64bit=false, add_APS=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. if add_APS is true it will add a new field to all particles (needed in newer versions of LaMEM, but not supported in LaMEM 2.2.0 or before)
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.datIf 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.datGeophysicalModelGenerator.save_LaMEM_topography Method
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
GeophysicalModelGenerator.screenshot_to_CartData Method
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
GeophysicalModelGenerator.screenshot_to_GeoData Method
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).
GeophysicalModelGenerator.screenshot_to_UTMData Method
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
GeophysicalModelGenerator.server_down_for Method
server_down_for()How many seconds are left of leaving the server alone after every mirror failed to deliver a file, or 0.0 if it may be asked. The marker is a file in the cache rather than a variable so that separate processes – the parallel test workers – share it.
GeophysicalModelGenerator.server_entry Method
server_entry(dataset, res, reg)The line that gmt_data_server.txt holds for this set, as a named tuple. Everything that varies from set to set – how large a tile is, how the integers in it scale to metres, which coarser set fills the gaps – is published there, so it is read rather than assumed.
GeophysicalModelGenerator.server_index_file Method
server_index_file()The server's own index of its datasets: downloaded once into the cache and read from there after that. If it cannot be downloaded, the copy bundled with the package is used instead. The download is attempted once per session, not once per call: the index is consulted several times per import, and with the server down each attempt would cost a full round of the mirrors.
sourceGeophysicalModelGenerator.setup_model_domain Method
setup_model_domain(coord_x::Vector{Float64},
coord_y::Vector{Float64},
coord_z::Vector{Float64},
nx::Int, ny::Int, nz::Int,
n_ranks::Int) -> ModelDomainSetup model domain decomposition using domain boundaries and resolution.
Parameters:
coord_x, coord_y, coord_z: 2-element vectors specifying [min, max] for each direction
nx, ny, nz: Number of cells in each direction
n_ranks: Total number of MPI ranks
Returns:
- ModelDomain struct containing all domain decomposition information
GeophysicalModelGenerator.spacing Method
dlon, dlat = spacing(lon,lat)Computes the spacing with central differences
sourceGeophysicalModelGenerator.subset_indices Method
subset_indices(coords, lo, hi, px)The samples that make up lo..hi, following what GMT returns. With gridline registration the samples are the cell corners and GMT widens the request out to the nodes that enclose it; with pixel registration they are the cell centres and a cell belongs to the region when its centre does.
GeophysicalModelGenerator.subtract_horizontalmean Method
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
GeophysicalModelGenerator.subtract_horizontalmean Method
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
GeophysicalModelGenerator.swap_yz_dims Method
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.
GeophysicalModelGenerator.tile_directory Method
tile_directory(dataset, res, reg)The directory on the server that holds the tiles of this set.
sourceGeophysicalModelGenerator.tile_filename Method
tile_filename(lat, lon, dataset, res, reg)The name of the tile whose south-west corner is at (lat, lon).
The 3 arcsecond set is the odd one out: it is SRTM repackaged, and is named after that rather than after the dataset it is served under.
sourceGeophysicalModelGenerator.tile_size_deg Method
tile_size_deg(res)The width and height, in degrees, of one tile at resolution res, or nothing if that resolution is published as a single global file rather than as tiles.
GeophysicalModelGenerator.tile_size_px Method
tile_size_px(res, deg, reg)How many pixels a tile has along each side: the tile width in degrees divided by the spacing that res names (m = arcminutes, s = arcseconds), plus one for gridline registration.
GeophysicalModelGenerator.topo_cache_dir Method
topo_cache_dir()Where downloaded tiles are kept, so that repeated imports of the same region do not download them again. Uses a scratch space, which is removed when the package is.
sourceGeophysicalModelGenerator.update_nel_xyz_in_file Method
update_nel_xyz_in_file(ParamFile_name, nx, ny, nz) Generates or updates CPU lines in a LaMEM input file after the nel_z line. Usage:
- Generate LaMEM model
Target resolution is 256x256x256.
nx = 256; ny = 256; nz = 256 ParamFile_name = "LaMEM_setup.txt" model = Model( LaMEM.Grid(x=[-80.,80.], y=[-95.,100.], z=[-60.0,10.0] , nel=(8,8,8)),Output(param_file_name = ParamFile_name,)) Notice that here nel is a very small number. Thick hacky way prevents us to generate huge grid at this stage.
- Once populated model with information, we save model as setup file in text.
write_LaMEM_inputFile(model,ParamFile_name).
- Then we can use this function to update back amount of elements of each direction to the target one, so setup file is correct.
update_nel_xyz_in_file(ParamFile_name, nx, ny, nz)
sourceGeophysicalModelGenerator.velocity_spherical_to_cartesian! Method
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
GeophysicalModelGenerator.visualise Function
visualiseInteractive widget that allows you to explore a 3D data set DataSet in an interactive manner. It requires you to load GLMakie.
GeophysicalModelGenerator.votemap Method
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,)GeophysicalModelGenerator.voxel_grav Method
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 ofRHOis 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 thanrhoTolwill 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 [trueorfalse]
GeophysicalModelGenerator.write_FEmesh_fields Method
write_FEmesh_fields(data::FEData)Writes cell and vertex fields to disk to be read by pTatin
sourceGeophysicalModelGenerator.write_FEmesh_mesh Method
write_FEmesh_mesh(vdata::FEData; out_file="md.bin", connectivity_zero_based=true)Writes a binary file with the mesh information for pTatin
sourceGeophysicalModelGenerator.write_pTatin_mesh Method
write_pTatin_mesh(fe_mesh::FEData; out_file="md.bin", connectivity_zero_based=true)Write a binary file with the mesh information for pTatin
sourceGeophysicalModelGenerator.write_pTatin_mesh Method
write_pTatin_mesh(q1_mesh::Q1Data; out_file="md.bin", connectivity_zero_based=true)Write a binary file with the mesh information for pTatin
sourceGeophysicalModelGenerator.write_paraview Function
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
GeophysicalModelGenerator.write_paraview Function
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)GeophysicalModelGenerator.write_paraview Function
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
GeophysicalModelGenerator.write_paraview Method
write_paraview(DataSet::CartData, filename::Any; PointsData=false, pvd=nothing, time=nothing, directory=nothing, verbose=true)Writes a CartData structure to paraview.
GeophysicalModelGenerator.write_paraview Method
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.
GeophysicalModelGenerator.write_processor_partitioning_LaMEM Method
write_processor_partitioning_LaMEM(P::LaMEMPartitioningInfo; is64bit::Bool = false)
Writes the processor partitioning information P to a binary file in the format used by LaMEM. The file is named ProcessorPartitioning_(P.nProcX*P.nProcY*P.nProcZ)cpu_(P.nProcX).(P.nProcY).(P.nProcZ).bin and returned as string. The coordinates are written as Float64, and the processor counts and node counts as Int64 or Int32, depending on the is64bit flag.
GeophysicalModelGenerator.xyz_grid Method
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.
WhereTheWaterFlows.waterflows Function
Topo_water, sinks, pits, bnds = waterflows(Topo::GeoData;
flowdir_fn=WhereTheWaterFlows.d8dir_feature, feedback_fn=nothing, drain_pits=true, bnd_as_sink=true,
rainfall = nothing,
minsize=300)Takes a GMG GeoData object of a topographic map and routes water through the grid. Optionally, you can specify rainfall in which case we accumulate the rain as specified in this 2D array instead of the cellarea. This allows you to, for example, sum, up water if you have variable rainfall in the area. The other options are as in the waterflows function of the package WhereTheWaterFlows.
Example
# Download some topographic data
julia> Topo = import_topo([6.5,7.3,50.2,50.6], file="@earth_relief_03s");
# Flow the water through the area:
julia> Topo_water, sinks, pits, bnds = waterflows(Topo)
julia> Topo_water
GeoData
size : (961, 481, 1)
lon ϵ [ 6.5 : 7.3]
lat ϵ [ 50.2 : 50.59999999999999]
depth ϵ [ 0.045 : 0.724]
fields : (:Topography, :area, :slen, :dir, :nout, :nin, :c)