Tools

We have a number of functions with which we can extract sub-data from a 2D or 3D GeoData structure.

Missing docstring.

Missing docstring for cross_section. Check Documenter's build log for details.

GeophysicalModelGenerator.extract_subvolume Function

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

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

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

• Lon_level,Lat_level and Depth_level should be tuples that indicate (minimum_value, maximum_value) along the respective direction. If not specified we use the full range.

• By default, Interpolate=false and we find the closest indices within the data set (so your new data set will not go exactly from minimum to maximum).

• Alternatively, if Interpolate=true we interpolate the data onto a new grid that has dimensions dims. This can be useful to compare data sets that are originally given in different resolutions.

3D Example with no interpolation:

julia> Lon,Lat,Depth   =   lonlatdepth_grid(10:20,30:40,(-300:25:0)km);
julia> Data            =   Depth*2;                # some data
julia> Vx,Vy,Vz        =   ustrip(Data*3),ustrip(Data*4),ustrip(Data*5);
julia> Data_set3D      =   GeoData(Lon,Lat,Depth,(Depthdata=Data,LonData=Lon, Velocity=(Vx,Vy,Vz)))
GeoData
  size  : (11, 11, 13)
  lon   ϵ [ 10.0 : 20.0]
  lat   ϵ [ 30.0 : 40.0]
  depth ϵ [ -300.0 km : 0.0 km]
  fields: (:Depthdata, :LonData, :Velocity)
julia> Data_extracted = extract_subvolume(Data_set3D,Lon_level=(10,12),Lat_level=(35,40))
GeoData
  size  : (3, 6, 13)
  lon   ϵ [ 10.0 : 12.0]
  lat   ϵ [ 35.0 : 40.0]
  depth ϵ [ -300.0 km : 0.0 km]
  fields: (:Depthdata, :LonData, :Velocity)

By default it extracts the data points closest to the area defined by 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)

source

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

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

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

• Lon_level,Lat_level and Depth_level should be tuples that indicate (minimum_value, maximum_value) along the respective direction. If not specified we use the full range.

• By default, Interpolate=false and we find the closest indices within the data set (so your new data set will not go exactly from minimum to maximum).

• Alternatively, if Interpolate=true we interpolate the data onto a new grid that has dimensions dims. This can be useful to compare data sets that are originally given in different resolutions.

3D Example with no interpolation:

julia> Lon,Lat,Depth   =   lonlatdepth_grid(10:20,30:40,(-300:25:0)km);
julia> Data            =   Depth*2;                # some data
julia> Vx,Vy,Vz        =   ustrip(Data*3),ustrip(Data*4),ustrip(Data*5);
julia> Data_set3D      =   GeoData(Lon,Lat,Depth,(Depthdata=Data,LonData=Lon, Velocity=(Vx,Vy,Vz)))
GeoData
  size  : (11, 11, 13)
  lon   ϵ [ 10.0 : 20.0]
  lat   ϵ [ 30.0 : 40.0]
  depth ϵ [ -300.0 km : 0.0 km]
  fields: (:Depthdata, :LonData, :Velocity)
julia> Data_extracted = extract_subvolume(Data_set3D,Lon_level=(10,12),Lat_level=(35,40))
GeoData
  size  : (3, 6, 13)
  lon   ϵ [ 10.0 : 12.0]
  lat   ϵ [ 35.0 : 40.0]
  depth ϵ [ -300.0 km : 0.0 km]
  fields: (:Depthdata, :LonData, :Velocity)

By default it extracts the data points closest to the area defined by 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}

source

GeophysicalModelGenerator.interpolate_datafields Function

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)

source

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

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

source

GeophysicalModelGenerator.votemap Function

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,)

source

GeophysicalModelGenerator.subtract_horizontalmean Function

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

source

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

source

Missing docstring.

Missing docstring for above_surface. Check Documenter's build log for details.

Missing docstring.

Missing docstring for below_surface. Check Documenter's build log for details.

Missing docstring.

Missing docstring for interpolate_data_surface. Check Documenter's build log for details.

GeophysicalModelGenerator.parse_columns_CSV Function

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)

source

GeophysicalModelGenerator.rotate_translate_scale Function

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

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

Parameters

Note that we apply the transformations in exactly this order:

• Scale: scaling applied to the x,y,z coordinates of the data set

• Rotate: rotation around the x/y axis (around the center of the box)

• Translate: translation

• Xc: center of rotation

Example

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

source

GeophysicalModelGenerator.point_to_nearest_grid Function

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

source

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

source

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

source

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

source

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)

source

Missing docstring.

Missing docstring for convert2UTMzone. Check Documenter's build log for details.

Missing docstring.

Missing docstring for convert2CartData. Check Documenter's build log for details.

Missing docstring.

Missing docstring for project_CartData. Check Documenter's build log for details.

Missing docstring.

Missing docstring for drape_on_topo. Check Documenter's build log for details.

GeophysicalModelGenerator.lithostatic_pressure! Function

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.

source

GeophysicalModelGenerator.countmap Function

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

source