Tools

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

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

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

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

3D Example with no interpolation:

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

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

3D Example with interpolation:

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

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

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

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

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

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

3D Example with no interpolation:

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

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

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

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

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

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

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

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

Example

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

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

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

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

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

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

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

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

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

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

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

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)

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

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

Point should have 1D coordinate vectors

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

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

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

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

Grid_counts = point_to_nearest_grid(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

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

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)

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.

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