/
Tutorial_LaPalma.jl
94 lines (72 loc) · 4.62 KB
/
Tutorial_LaPalma.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
# # Create a Cartesian Model Setup for La Palma
#
# ## Aim
# In this tutorial, your will learn how to use real data to create a geodynamic model setup in `CartData`.
# We will use the data of the Cumbre Viejo eruption in La Palma, which is a volcanic island that erupted from mid september 2021 - december 2021.
# The seimsicity from that time-period is used as an inspiration to set the locations of magma intrusions
#
# ## 1. Load data
# We will use two types of data to create the model
# 1) Topography
# 2) Earthquake locations
# We start with loading the required packages, which includes `GMT` to download topography (an optional dependency for `GeophysicalModelGenerator`)
using GeophysicalModelGenerator, GMT, DelimitedFiles
# We will use GMT to download the topography with:
Topo = import_topo(lon = [-18.2, -17.5], lat=[28.4, 29.0], file="@earth_relief_15s.grd")
# Next, lets load the seismicity. The earthquake data is available on [https://www.ign.es/web/ign/portal/vlc-catalogo](https://www.ign.es/web/ign/portal/vlc-catalogo).
# We have filtered them and prepared a file with earthquake locations up to early November 2021 (from january 2021).
# Download that:
download_data("https://zenodo.org/records/10738510/files/EQ_events_all_info5_LaPalma_2021.dat","EQ_events_all_info5_LaPalma_2021.dat")
data_EQ = readdlm("EQ_events_all_info5_LaPalma_2021.dat")
# It has the following format:
lon = data_EQ[:,1]
lat = data_EQ[:,2]
depth = -data_EQ[:,3]
Mag = data_EQ[:,10]
# Create a `GeoData` structure from this:
data_all_EQ = GeoData(lon,lat,depth, (Magnitude=Mag,))
# Next, we can write the data to paraview along with the topography. Note that we have to specify that we have `PointData`:
write_paraview(data_all_EQ,"data_all_EQ",PointsData=true)
write_paraview(Topo,"Topo")
# 
# Note that this data is in geographic coordinates, which makes it non-trivial to create slices through the data (see coordinate axis in the plot, where `z` is *not* pointing upwards).
# ## 2. Convert data to cartesian coordinates
# In order to create model setups, it is helpful to first transfer the data to Cartesian.
# This requires us to first determine a *projection point*, that is fixed. Often, it is helpful to use the center of the topography for this. In the present example, we will center the model around La Palma itself:
proj = ProjectionPoint(Lon=-17.84, Lat=28.56)
# Once this is done you can convert the topographic data to the cartesian reference frame
EQ_cart = convert2CartData(data_all_EQ, proj);
Topo_cart = convert2CartData(Topo, proj)
# It is important to realize that the cartesian coordinates of the topographic grid is no longer strictly orthogonal after this conversion. You don't notice that in the current example, as the model domain is rather small.
# In other cases, however, this is quite substantial (e.g., India-Asia collision zone).
# LaMEM needs an orthogonal grid of topography, which we can create with:
Topo_model = CartData(xyz_grid(-35:.1:30,-15:.2:45,0));
# In a next step, the routine `project_CartData` projects a `GeoData` structure to a `CartData` struct
Topo_model = project_CartData(Topo_model, Topo, proj)
# Let's have a look at the data:
write_paraview(EQ_cart,"EQ_cart",PointsData=true)
write_paraview(Topo_model,"Topo_model")
# ## 3. Create a volumetric earthquake plot
# It is useful to plot the earthquake density in 3D, which indicates where most action is happening in the system.
# For this, we first create a 3D grid of the region:
Grid_3D = CartData(xyz_grid(-35:.3:30,-15:.25:45,-50:.5:5))
# Next we check how many earthquakes are around the grid points:
Grid_3D =point_to_nearest_grid(EQ_cart, Grid_3D, radius_factor=3)
# And we can define an array with rock types:
Phases = zeros(Int64,size(Grid_3D.x))
# Points that are below the surface are set to one:
Below = below_surface(Grid_3D, Topo_model);
Phases[Below] .= 1
# Lets assume that the crust is 15 km thick
Phases[NumValue(Grid_3D.z) .< -15] .= 2
# And lets assume that the magma is where we have some earthquake activity:
ind = findall( (Grid_3D.fields.Count .> 75) .&& Phases.>0)
Phases[ind] .= 3 #Magma
# Add rocktypes to the grid:
Grid_3D = addfield(Grid_3D,"Phases",Phases)
# We can save this to paraview format
write_paraview(Grid_3D,"Grid_3D")
# The paraview statefile `/tutorials/LaPalma.pvsm` can be used to reproduce the following plot:
# 
#src Note: The markdown page is generated using:
#src Literate.markdown("tutorials/Tutorial_LaPalma.jl","docs/src/man",keepcomments=true, execute=false, codefence = "```julia" => "```")