Data structures
The main data structure used in GeophysicalModelGenerator.jl is GeoData, which contains info about the longitude,latitude, and depth of a data set, as well as several data sets itself. We also provide a UTMData, which is essentially the same but with UTM coordinates, and a CartData structure, which has Cartesian coordinates in kilometers (as used in many geodynamic codes). If one wishes to transfer GeoData to CartData, one needs to provide a ProjectionPoint. For plotting, we transfer this into the ParaviewData structure, which has Cartesian coordinates centered around the center of the Earth. We employ the wgs84 reference ellipsoid as provided by the Geodesy.jl package to perform this transformation.
GeophysicalModelGenerator.GeoData — TypeGeoData(lon::Any, lat:Any, depth::GeoUnit, fields::NamedTuple)Data structure that holds one or several fields with longitude, latitude and depth information.
depthcan have units of meter, kilometer or be unitless; it will be converted to km.fieldsshould ideally be a NamedTuple which allows you to specify the names of each of the fields.- In case you only pass one array we will convert it to a NamedTuple with default name.
- A single field should be added as
(DataFieldName=Data,)(don't forget the comma at the end). - Multiple fields can be added as well.
lon,lat,depthshould all have the same size as each of thefields. - In case you want to display a vector field in paraview, add it as a tuple:
(Velocity=(Veast,Vnorth,Vup), Veast=Veast, Vnorth=Vnorth, Vup=Vup); we automatically apply a vector transformation when transforming this to aParaviewDatastructure from which we generate Paraview output. As this changes the magnitude of the arrows, you will no longer see the[Veast,Vnorth,Vup]components in Paraview which is why it is a good ideas to store them as separate Fields. - Yet, there is one exception: if the name of the 3-component field is
colors, we do not apply this vector transformation as this field is regarded to contain RGB colors. Lat,Lon,Depthshould have the same size as theDataarray. The ordering of the arrays is important. If they are 3D arrays, as in the example below, we assume that the first dimension corresponds tolon, second dimension tolatand third dimension todepth(which should be in km). See below for an example.
Example
julia> Lat = 1.0:3:10.0;
julia> Lon = 11.0:4:20.0;
julia> Depth = (-20:5:-10)*km;
julia> Lon3D,Lat3D,Depth3D = lonlatdepth_grid(Lon, Lat, Depth);
julia> Lon3D
3×4×3 Array{Float64, 3}:
[:, :, 1] =
11.0 11.0 11.0 11.0
15.0 15.0 15.0 15.0
19.0 19.0 19.0 19.0
[:, :, 2] =
11.0 11.0 11.0 11.0
15.0 15.0 15.0 15.0
19.0 19.0 19.0 19.0
[:, :, 3] =
11.0 11.0 11.0 11.0
15.0 15.0 15.0 15.0
19.0 19.0 19.0 19.0
julia> Lat3D
3×4×3 Array{Float64, 3}:
[:, :, 1] =
1.0 4.0 7.0 10.0
1.0 4.0 7.0 10.0
1.0 4.0 7.0 10.0
[:, :, 2] =
1.0 4.0 7.0 10.0
1.0 4.0 7.0 10.0
1.0 4.0 7.0 10.0
[:, :, 3] =
1.0 4.0 7.0 10.0
1.0 4.0 7.0 10.0
1.0 4.0 7.0 10.0
julia> Depth3D
3×4×3 Array{Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(km,), 𝐋, nothing}}, 3}:
[:, :, 1] =
-20.0 km -20.0 km -20.0 km -20.0 km
-20.0 km -20.0 km -20.0 km -20.0 km
-20.0 km -20.0 km -20.0 km -20.0 km
[:, :, 2] =
-15.0 km -15.0 km -15.0 km -15.0 km
-15.0 km -15.0 km -15.0 km -15.0 km
-15.0 km -15.0 km -15.0 km -15.0 km
[:, :, 3] =
-10.0 km -10.0 km -10.0 km -10.0 km
-10.0 km -10.0 km -10.0 km -10.0 km
-10.0 km -10.0 km -10.0 km -10.0 km
julia> Data = zeros(size(Lon3D));
julia> Data_set = GeophysicalModelGenerator.GeoData(Lon3D,Lat3D,Depth3D,(DataFieldName=Data,))
GeoData
size : (3, 4, 3)
lon ϵ [ 11.0 : 19.0]
lat ϵ [ 1.0 : 10.0]
depth ϵ [ -20.0 km : -10.0 km]
fields : (:DataFieldName,)
attributes: ["note"]GeophysicalModelGenerator.CartData — TypeCartData(x::Any, y::Any, z::GeoUnit, fields::NamedTuple)Data structure that holds one or several fields with with Cartesian x/y/z coordinates. Distances are in kilometers
x,y,zcan have units of meters, kilometer or be unitless; they will be converted to kilometersfieldsshould ideally be a NamedTuple which allows you to specify the names of each of the fields.- In case you only pass one array we will convert it to a NamedTuple with default name.
- A single field should be added as
(DataFieldName=Data,)(don't forget the comma at the end). - Multiple fields can be added as well.
- In case you want to display a vector field in paraview, add it as a tuple:
(Velocity=(Vx,Vnorth,Vup), Veast=Veast, Vnorth=Vnorth, Vup=Vup); we automatically apply a vector transformation when transforming this to aParaviewDatastructure from which we generate Paraview output. As this changes the magnitude of the arrows, you will no longer see the[Veast,Vnorth,Vup]components in Paraview which is why it is a good ideas to store them as separate Fields. - Yet, there is one exception: if the name of the 3-component field is
colors, we do not apply this vector transformation as this field is regarded to contain RGB colors. x,y,zshould have the same size as theDataarray. The ordering of the arrays is important. If they are 3D arrays, as in the example below, we assume that the first dimension corresponds tox, second dimension toyand third dimension toz(which should be in km). See below for an example.
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_set = 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"]CartData is particularly useful in combination with cartesian geodynamic codes, such as LaMEM, which require cartesian grids. You can directly save your data to Paraview with
julia> write_paraview(Data_set, "Data_set")
1-element Vector{String}:
"Data_set.vts"If you wish, you can convert this to UTMData (which will simply convert the )
julia> Data_set1 = convert(GeoData, Data_set)
GeoData
size : (116, 96, 25)
lon ϵ [ 14.075969111533457 : 14.213417764154963]
lat ϵ [ 40.77452227533946 : 40.86110443583479]
depth ϵ [ -5.4 km : 0.6 km]
fields: (:FakeData, :Data2)which would allow visualizing this in paraview in the usual manner:
GeophysicalModelGenerator.UTMData — TypeUTMData(EW::Any, NS:Any, depth::GeoUnit, UTMZone::Int, NorthernHemisphere=true, fields::NamedTuple)Data structure that holds one or several fields with UTM coordinates (east-west), (north-south) and depth information.
depthcan have units of meters, kilometer or be unitless; it will be converted to meters (as UTMZ is usually in meters)fieldsshould ideally be a NamedTuple which allows you to specify the names of each of the fields.- In case you only pass one array we will convert it to a NamedTuple with default name.
- A single field should be added as
(DataFieldName=Data,)(don't forget the comma at the end). - Multiple fields can be added as well.
- In case you want to display a vector field in paraview, add it as a tuple:
(Velocity=(Veast,Vnorth,Vup), Veast=Veast, Vnorth=Vnorth, Vup=Vup); we automatically apply a vector transformation when transforming this to aParaviewDatastructure from which we generate Paraview output. As this changes the magnitude of the arrows, you will no longer see the[Veast,Vnorth,Vup]components in Paraview which is why it is a good ideas to store them as separate Fields. - Yet, there is one exception: if the name of the 3-component field is
colors, we do not apply this vector transformation as this field is regarded to contain RGB colors. Lat,Lon,Depthshould have the same size as theDataarray. The ordering of the arrays is important. If they are 3D arrays, as in the example below, we assume that the first dimension corresponds tolon, second dimension tolatand third dimension todepth(which should be in km). See below for an example.
Example
julia> ew = 422123.0:100:433623.0
julia> ns = 4.514137e6:100:4.523637e6
julia> depth = -5400:250:600
julia> EW,NS,Depth = xyz_grid(ew, ns, depth);
julia> Data = ustrip.(Depth);
julia> Data_set = UTMData(EW,NS,Depth,33, true, (FakeData=Data,Data2=Data.+1.))
UTMData
UTM zone : 33-33 North
size : (116, 96, 25)
EW ϵ [ 422123.0 : 433623.0]
NS ϵ [ 4.514137e6 : 4.523637e6]
depth ϵ [ -5400.0 m : 600.0 m]
fields : (:FakeData, :Data2)
attributes: ["note"]If you wish, you can convert this from UTMData to GeoData with
julia> Data_set1 = convert(GeoData, Data_set)
GeoData
size : (116, 96, 25)
lon ϵ [ 14.075969111533457 : 14.213417764154963]
lat ϵ [ 40.77452227533946 : 40.86110443583479]
depth ϵ [ -5.4 km : 0.6 km]
fields : (:FakeData, :Data2)
attributes: ["note"]which would allow visualizing this in paraview in the usual manner:
julia> write_paraview(Data_set1, "Data_set1")
1-element Vector{String}:
"Data_set1.vts"GeophysicalModelGenerator.Q1Data — TypeHolds a Q1 Finite Element Data set with vertex and cell data. The specified coordinates are the ones of the vertices.
GeophysicalModelGenerator.FEData — TypeFEData{dim, points_per_cell}Structure that holds Finite Element info with vertex and cell data. Works in 2D/3D for arbitrary elements
Parameters
verticeswith the points on the mesh (dimxNpoints)connectivitywith the connectivity of the mesh (points_per_cellxNcells)fieldswith the fields on the verticescellfieldswith the fields of the cells
GeophysicalModelGenerator.ParaviewData — TypeParaviewData(x::GeoUnit, y::GeoUnit, z::GeoUnit, values::NamedTuple)Cartesian data in x/y/z coordinates to be used with Paraview. This is usually generated automatically from the GeoData structure, but you can also invoke do this manually:
julia> Data_set = GeophysicalModelGenerator.GeoData(1.0:10.0,11.0:20.0,(-20:-11)*km,(DataFieldName=(-20:-11),))
julia> Data_cart = convert(ParaviewData, Data_set)GeophysicalModelGenerator.lonlatdepth_grid — FunctionLon, 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.xyz_grid — FunctionX,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.
GeophysicalModelGenerator.ProjectionPoint — Typestruct ProjectionPoint
Lon :: Float64
Lat :: Float64
EW :: Float64
NS :: Float64
zone :: Integer
isnorth :: Bool
endStructure that holds the coordinates of a point that is used to project a data set from Lon/Lat to a Cartesian grid and vice-versa.