mod_util_functions/plot_model_mod.py at main · joaopaulomo/mod_util_functions · GitHub

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
from typing import List
from examples.seismic import Model

import matplotlib
from matplotlib import pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable

import numpy as np


def plot_density(model: Model, show: bool = True):
    """
    Plot a two-dimensional density field from a seismic Model
    object.

    Parameters
    ----------
    model : Model
        Object that holds the density model.
    """
    rho = 1/model.b.data
    rho = rho[
        model.nbl:-model.nbl,
        model.nbl:-model.nbl
    ]
    domain_size = 1.e-3 * np.array(model.domain_size)
    extent = [model.origin[0], model.origin[0] + domain_size[0],
              model.origin[1] + domain_size[1], model.origin[1]]

    plt.figure()
    plot = plt.imshow(np.transpose(rho), cmap='cividis', vmin=np.min(rho), vmax=np.max(rho),
                      extent=extent)
    plt.xlabel('X position (km)')
    plt.ylabel('Depth (km)')

    ax = plt.gca()
    div = make_axes_locatable(ax)
    cax = div.append_axes('right', '5%', 0.05)
    cbar = plt.colorbar(plot, cax=cax)
    cbar.set_label('Density ($g/$cm$^3$)')

    if show:
        plt.show()


import numpy as np
import matplotlib.pyplot as plt
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable

def plot_image2(data, model, vmin=None, vmax=None, colorbar=True, cmap="gray",
                show: bool = True, clip_percent=None, clip_low=None):
    """
    Works like plot_image, but adjusts the image aspect ratio and clips high values
    based on a percentile, while keeping the range dynamic for better contrast.

    Parameters
    ----------
    data : ndarray
        Image data to plot.
    model : SeismicModel
        Reference model to get the aspect ratio.
    cmap : str
        Choice of colormap. Defaults to gray scale for images as a seismic convention.
    clip_percent : float, optional
        Percentage threshold for clipping high values.
    clip_low : float, optional
        Percentage threshold for clipping low values.
    """
    plt.figure()

    # Calculate domain size and extent for the plot
    domain_size = 1.e-3 * np.array(model.domain_size)
    extent = [
        model.origin[0], model.origin[0] + domain_size[0],
        model.origin[1] + domain_size[1], model.origin[1]
    ]

    print("Original data shape:", data.shape)
    data = np.asarray(data)


    if clip_percent is not None or clip_low is not None:
        low_value = np.percentile(data, clip_low) if clip_low is not None else np.min(data)
        high_value = np.percentile(data, clip_percent) if clip_percent is not None else np.max(data)
        print(f"Clipping low values below {low_value:.5f} and high values above {high_value:.5f}.")
        data = np.clip(data, low_value, high_value)

        if vmin is None:
            vmin = low_value
        if vmax is None:
            vmax = high_value
    else:
        print("Clipping not applied; plotting data as is.")

    plot = plt.imshow(
        np.transpose(data),
        vmin=vmin,
        vmax=vmax,
        cmap=cmap,
        extent=extent
    )

    # Create aligned colorbar on the right
    if colorbar:
        ax = plt.gca()
        divider = make_axes_locatable(ax)
        cax = divider.append_axes("right", size="5%", pad=0.05)
        plt.colorbar(plot, cax=cax)

    if show:
        plt.show()


def plot_graph(xdata: np.ndarray, ydata: np.ndarray, xlabel: str, ylabel: str,
               title: str = '', show: bool = True):
    plt.figure()
    plt.plot(xdata, ydata)
    plt.xlabel(xlabel)
    plt.ylabel(ylabel)
    plt.title(title)

    if show:
        plt.show()


def plot_array(arr: np.ndarray, xlabel: str, ylabel: str,
               title: str = '', show: bool = True):
    x = np.arange(1, 1 + len(arr))
    plot_graph(x, arr, xlabel, ylabel, title, show)
    plt.xticks(x[::int(np.ceil(len(arr)/10))])


def plot_trace(traces: List[np.ndarray], ymin, ymax, styles: List[str], legends: List[str],
                xlabel: str, ylabel: str, show: bool = True):
    plt.figure(figsize = (8, 9))

    for i in range(len(traces)):
        y = np.linspace(ymin, ymax, traces[i].size)
        plt.plot(traces[i], y, styles[i], linewidth=2)

    plt.legend(legends, fontsize=15)
    plt.ylabel(xlabel)
    plt.xlabel(ylabel)
    plt.gca().invert_yaxis()

    if show:
        plt.show()