-
Notifications
You must be signed in to change notification settings - Fork 3
/
Copy pathplot_example.jl
179 lines (147 loc) · 5.52 KB
/
plot_example.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
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
#' # Plotting data and models with SlimPlotting
#' ---
#' title: Overview of SlimPlotting utilities
#' author: Mathias Louboutin
#' date: April 2023
#' ---
#' This example script is written using [Weave.jl](https://github.com/JunoLab/Weave.jl) and can be converted to different format for documentation and usage
#' This example is converted to a markdown file for the documentation.
#' # Import SlimPlotting, SegyIO to read seismic data, JLD2 for hdf5-like files
using PythonPlot, SlimPlotting, SegyIO, JLD2
#' # Initialize all needed data
#' Close all figures if any existing
plotclose("all")
#' Path to the files and data used for these examples
data_path = dirname(pathof(SlimPlotting)) * "/../data/";
#' Read the data
# Pure array
vp = Float32.(segy_read("$(data_path)2dVP.sgy").data);
dm = diff(vp, dims = 1);
shot = Float32.(segy_read("$(data_path)2dshot.segy").data);
xloc = get_header(segy_read("$(data_path)2dshot.segy"), "GroupX")
fslice = JLD2.load("$(data_path)2dfslice.jld");
#' # Create structures to mimic JUDI-like inputs
#' In the future this should be instead converted into an extension rather than implicit knowledge of the structure
# Dummy structures to check plot with metadata
struct geometry
xloc::Any
end
struct shotrec
data::Any
dt::Any
geometry::Any
end
struct Phys
data::Any
d::Any
end
## Make physical objects
dmp = Phys(dm, (10, 20))
vpp = Phys(vp, (10, 20))
fslicep = Phys(fslice["Freq"][1, :, :], (12.5, 12.5))
shotp = shotrec([shot], 0.008, geometry([xloc]));
#' # Model perturbation
#' We plot here a model perturbation (i.e a Reverse-time Migrated image) and compare a few colormaps:
#' - The perceptually accurate `Greys` colormap from colorcet
#' - The standard matplotlib `Greys` colormap
#' - Another perceptually accurate `Greys` colormap from colorcet
figure(figsize = (10, 10))
subplot(311)
plot_simage(dmp; new_fig = false, name = "Colorcet Greys")
subplot(312)
plot_simage(dm, (10, 20); cmap = "Greys", new_fig = false, name = "Greys")
subplot(313)
plot_simage(dm, (10, 20); cmap = :cet_CET_L2, new_fig = false, name = "Colorcet Greys 2")
tight_layout();
display(gcf());
#' # Velocity
#' We plot here a velocity model and compare a few colormaps:
#' - The perceptually accurate `jet` colormap from colorcet named `cet_rainbow4`
#' - The ColorSchemes `vik` colormap
#' - The matplotlib `jet`
figure(figsize = (10, 10))
subplot(311)
plot_velocity(vpp; new_fig = false, name = "colorcet jet", cmap = "cet_rainbow4")
subplot(312)
plot_velocity(vp, (10, 20); cmap = :vik, new_fig = false, name = "ColorSchemes's vik")
subplot(313)
plot_velocity(vp, (10, 20); cmap = "jet", new_fig = false, name = "matplotlib jet")
tight_layout();
display(gcf());
#' # Frequency slice
#' We plot here a frequency slice for a seismic dataset and compare a few colormaps:
#' - The perceptually accurate `bwr` colormap from colorcet named `cet_CET_D1A`
#' - The standard matplotlib `bwr` colormap
#' - The matplotlib "bwr
# Frequency slice
figure(figsize = (10, 5))
subplot(131)
plot_fslice(fslice["Freq"][1, :, :], (12.5, 12.5); new_fig = false, name = "colorcet bwr")
subplot(132)
plot_fslice(fslicep; cmap = :bwr, new_fig = false, name = "bwr")
subplot(133)
plot_fslice(fslicep; cmap = "bwr", new_fig = false, name = "Matplotlib bwr")
tight_layout();
display(gcf());
#' # Shot record
#' ## Seismic blue-white-red
#'
#' We plot here a frequency slice for a seismic dataset and compare a few colormaps for the `bwr` colormap:
#' - The standard matplotlib `bwr` colormap
#' - The perceptually accurate `bwr` colormap from colorcet named `cet_CET_D1A` and `cet_CET_D1`
# Shot record
figure(figsize = (10, 5))
subplot(131)
plot_sdata(shotp; new_fig = false, name = "matplotlib seismic", cmap = "bwr")
subplot(132)
plot_sdata(shot, (12.5, 0.008); cmap = :cet_CET_D1A, new_fig = false, name = "Colorcet bwr")
subplot(133)
plot_sdata(shot, (12.5, 0.008); cmap = :cet_CET_D1, new_fig = false, name = "Colorcet bwr 2")
tight_layout();
display(gcf());
#' ## Seismic greys
#' We plot here a frequency slice for a seismic dataset and compare a few colormaps for the `greys` colormap:
#' - The standard matplotlib `gray` colormap
#' - The perceptually accurate `greys` colormap from colorcet named `cet_CET_L1`
# Shot record
figure(figsize = (10, 5))
subplot(121)
plot_sdata(shotp; new_fig = false, name = "colorcet gray", cmap = "cet_CET_L1")
subplot(122)
plot_sdata(shot, (12.5, 0.008); cmap = "gray", new_fig = false, name = "Greys")
tight_layout();
display(gcf());
#' # Compare shot records
#' One of the main visual representation of FWI inversion is to compare the true shot record with the synthetic data from
#' the current velocity model. A good way to visualize this difference is to overlay the two shot records alternating the traces
#' between each shots with a different colormap to check the alignment of the events. We show below how to do this with the
#' `compare_shots` function
figure(figsize = (10, 5))
subplot(131)
compare_shots(shotp, shotp; new_fig = false, name = "Overlap compare")
subplot(132)
compare_shots(
shotp,
shotp;
new_fig = false,
cmap = ("bwr", "RdBu"),
name = "Overlap compare custom cmap",
)
subplot(133)
compare_shots(
shotp,
shotp;
side_by_side = true,
new_fig = false,
name = "Side by side compare",
)
tight_layout();
display(gcf());
#' # Wiggle traces
#' We finally show the traditional wiggle plot for a shot record used in seismic.
# Wiggle plot
figure(figsize=(5, 5))
subplot(111)
wiggle_plot(shot[1:5:end, 1:10:end], xloc[1:10:end], 0:0.02:4.6; new_fig=false)
tight_layout();
display(gcf())