Tools

CrossSection(Volume::GeoData; dims=(100,100), Interpolate=false, Depth_level=nothing; Lat_level=nothing; Lon_level=nothing; Start=nothing, End=nothing )

Creates a cross-section through a volumetric (3D) GeoData object.

• Cross-sections can be horizontal (map view at a given depth), if Depth_level is specified
• They can also be vertical, either by specifying Lon_level or Lat_level (for a fixed lon/lat), or by defining both Start=(lon,lat) & End=(lon,lat) points.
• Interpolate indicates whether we want to simply extract the data from the 3D volume (default) or whether we want to linearly interpolate it on a new grid, which has dimensions as specified in dims

Example:

julia> Lon,Lat,Depth   =   LonLatDepthGrid(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      =   CrossSection(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)

ExtractSubvolume(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   =   LonLatDepthGrid(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 = ExtractSubvolume(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 = ExtractSubvolume(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)

InterpolateDataFields(V::GeoData, Lon, Lat, Depth)

Interpolates a data field V on a grid defined by Lon,Lat,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 = SubtractHorizontalMean(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 = SubtractHorizontalMean(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

AboveSurface(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   =   LonLatDepthGrid(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   =   LonLatDepthGrid(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       =   AboveSurface(Data_set3D, Data_Moho); 

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

Above = AboveSurface(Data_Cart::CartData, DataSurface_Cart::CartData; above=true)

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

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

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

Below = BelowSurface(Data_Cart::CartData, DataSurface_Cart::CartData)

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

Surf_interp = InterpolateDataOnSurface(V::CartData, Surf::CartData)

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

Example

julia> Data
CartData 
  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
CartData 
  size  : (96, 96, 1)
  x     ϵ [ -2.9671875 : 3.2671875]
  y     ϵ [ -1.9791666666666667 : 1.9791666666666667]
  z     ϵ [ -1.5353766679763794 : -0.69925457239151]
  fields: (:Depth,)
julia> Surf_interp = InterpolateDataOnSurface(Data, surf)
  CartData 
    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)

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

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

ParseColumns_CSV_File(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 intested in the numbers

Example

This example assumes that the data starts at line 18, that the colums 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 = ParseColumns_CSV_File(data_file, 4)

RotateTranslateScale!(Data::CartData; Rotate=0, Translate=(0,0,0), Scale=(1.0,1.0,1.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

Example

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