PySDM_examples.Shima_et_al_2009.spectrum_plotter

  1import matplotlib
  2import numpy as np
  3from matplotlib import pyplot
  4from open_atmos_jupyter_utils import show_plot
  5from packaging import version
  6from PySDM_examples.Shima_et_al_2009.error_measure import error_measure
  7
  8from PySDM.physics.constants import si
  9
 10_matplotlib_version_3_3_3 = version.parse("3.3.0")
 11_matplotlib_version_actual = version.parse(matplotlib.__version__)
 12
 13
 14class SpectrumColors:
 15    def __init__(self, begining="#2cbdfe", end="#b317b1"):
 16        self.b = begining
 17        self.e = end
 18
 19    def __call__(self, value: float):
 20        bR, bG, bB = int(self.b[1:3], 16), int(self.b[3:5], 16), int(self.b[5:7], 16)
 21        eR, eG, eB = int(self.e[1:3], 16), int(self.e[3:5], 16), int(self.e[5:7], 16)
 22        R = bR + int((eR - bR) * value)
 23        G = bG + int((eG - bG) * value)
 24        B = bB + int((eB - bB) * value)
 25        result = f"#{hex(R)[2:4]}{hex(G)[2:4]}{hex(B)[2:4]}"
 26        return result
 27
 28
 29class SpectrumPlotter:
 30    def __init__(self, settings, title=None, grid=True, legend=True, log_base=10):
 31        self.settings = settings
 32        self.format = "pdf"
 33        self.colors = SpectrumColors()
 34        self.smooth = False
 35        self.smooth_scope = 2
 36        self.legend = legend
 37        self.grid = grid
 38        self.title = title
 39        self.xlabel = "particle radius [µm]"
 40        self.ylabel = "dm/dlnr [g/m^3/(unit dr/r)]"
 41        self.log_base = log_base
 42        self.ax = pyplot
 43        self.fig = pyplot
 44        self.finished = False
 45
 46    def finish(self):
 47        if self.finished:
 48            return
 49        self.finished = True
 50        if self.grid:
 51            self.ax.grid()
 52
 53        base_arg = {
 54            "base"
 55            + (
 56                "x" if _matplotlib_version_actual < _matplotlib_version_3_3_3 else ""
 57            ): self.log_base
 58        }
 59        if self.title is not None:
 60            try:
 61                self.ax.title(self.title)
 62            except TypeError:
 63                self.ax.set_title(self.title)
 64        try:
 65            self.ax.xscale("log", **base_arg)
 66            self.ax.xlabel(self.xlabel)
 67            self.ax.ylabel(self.ylabel)
 68        except AttributeError:
 69            self.ax.set_xscale("log", **base_arg)
 70            self.ax.set_xlabel(self.xlabel)
 71            self.ax.set_ylabel(self.ylabel)
 72        if self.legend:
 73            self.ax.legend()
 74
 75    def show(self):
 76        self.finish()
 77        pyplot.tight_layout()
 78        show_plot()
 79
 80    def save(self, file):
 81        self.finish()
 82        pyplot.savefig(file, format=self.format)
 83
 84    def plot(
 85        self, spectrum, t, label=None, color=None, title=None, add_error_to_label=False
 86    ):
 87        error = self.plot_analytic_solution(self.settings, t, spectrum, title)
 88        if label is not None and add_error_to_label:
 89            label += f" error={error:.4g}"
 90        self.plot_data(self.settings, t, spectrum, label, color)
 91        return error
 92
 93    def plot_analytic_solution(self, settings, t, spectrum, title):
 94        if t == 0:
 95            analytic_solution = settings.spectrum.size_distribution
 96        else:
 97
 98            def analytic_solution(x):
 99                return settings.norm_factor * settings.kernel.analytic_solution(
100                    x=x, t=t, x_0=settings.X0, N_0=settings.n_part
101                )
102
103        volume_bins_edges = self.settings.formulae.trivia.volume(
104            settings.radius_bins_edges
105        )
106        dm = np.diff(volume_bins_edges)
107        dr = np.diff(settings.radius_bins_edges)
108
109        pdf_m_x = volume_bins_edges[:-1] + dm / 2
110        pdf_m_y = analytic_solution(pdf_m_x)
111
112        pdf_r_x = settings.radius_bins_edges[:-1] + dr / 2
113        pdf_r_y = pdf_m_y * dm / dr * pdf_r_x
114
115        x = pdf_r_x * si.metres / si.micrometres
116        y_true = (
117            pdf_r_y
118            * self.settings.formulae.trivia.volume(radius=pdf_r_x)
119            * settings.rho
120            / settings.dv
121            * si.kilograms
122            / si.grams
123        )
124
125        self.ax.plot(x, y_true, color="black")
126
127        if spectrum is not None:
128            y = spectrum * si.kilograms / si.grams
129            error = error_measure(y, y_true, x)
130            self.title = (
131                title or f"error measure: {error:.2f}"
132            )  # TODO #327 relative error
133            return error
134        return None
135
136    def plot_data(self, settings, t, spectrum, label, color):
137        if self.smooth:
138            scope = self.smooth_scope
139            if t != 0:
140                new = np.copy(spectrum)
141                for _ in range(2):
142                    for i in range(scope, len(spectrum) - scope):
143                        new[i] = np.mean(spectrum[i - scope : i + scope + 1])
144                    scope = 1
145                    for i in range(scope, len(spectrum) - scope):
146                        spectrum[i] = np.mean(new[i - scope : i + scope + 1])
147
148            x = settings.radius_bins_edges[:-scope]
149            dx = np.diff(x)
150            self.ax.plot(
151                (x[:-1] + dx / 2) * si.metres / si.micrometres,
152                spectrum[:-scope] * si.kilograms / si.grams,
153                label=label or f"t = {t}s",
154                color=color
155                or self.colors(t / (self.settings.output_steps[-1] * self.settings.dt)),
156            )
157        else:
158            self.ax.step(
159                settings.radius_bins_edges[:-1] * si.metres / si.micrometres,
160                spectrum * si.kilograms / si.grams,
161                where="post",
162                label=label or f"t = {t}s",
163                color=color
164                or self.colors(t / (self.settings.output_steps[-1] * self.settings.dt)),
165            )
class SpectrumColors:
15class SpectrumColors:
16    def __init__(self, begining="#2cbdfe", end="#b317b1"):
17        self.b = begining
18        self.e = end
19
20    def __call__(self, value: float):
21        bR, bG, bB = int(self.b[1:3], 16), int(self.b[3:5], 16), int(self.b[5:7], 16)
22        eR, eG, eB = int(self.e[1:3], 16), int(self.e[3:5], 16), int(self.e[5:7], 16)
23        R = bR + int((eR - bR) * value)
24        G = bG + int((eG - bG) * value)
25        B = bB + int((eB - bB) * value)
26        result = f"#{hex(R)[2:4]}{hex(G)[2:4]}{hex(B)[2:4]}"
27        return result
SpectrumColors(begining='#2cbdfe', end='#b317b1')
16    def __init__(self, begining="#2cbdfe", end="#b317b1"):
17        self.b = begining
18        self.e = end
b
e
class SpectrumPlotter:
 30class SpectrumPlotter:
 31    def __init__(self, settings, title=None, grid=True, legend=True, log_base=10):
 32        self.settings = settings
 33        self.format = "pdf"
 34        self.colors = SpectrumColors()
 35        self.smooth = False
 36        self.smooth_scope = 2
 37        self.legend = legend
 38        self.grid = grid
 39        self.title = title
 40        self.xlabel = "particle radius [µm]"
 41        self.ylabel = "dm/dlnr [g/m^3/(unit dr/r)]"
 42        self.log_base = log_base
 43        self.ax = pyplot
 44        self.fig = pyplot
 45        self.finished = False
 46
 47    def finish(self):
 48        if self.finished:
 49            return
 50        self.finished = True
 51        if self.grid:
 52            self.ax.grid()
 53
 54        base_arg = {
 55            "base"
 56            + (
 57                "x" if _matplotlib_version_actual < _matplotlib_version_3_3_3 else ""
 58            ): self.log_base
 59        }
 60        if self.title is not None:
 61            try:
 62                self.ax.title(self.title)
 63            except TypeError:
 64                self.ax.set_title(self.title)
 65        try:
 66            self.ax.xscale("log", **base_arg)
 67            self.ax.xlabel(self.xlabel)
 68            self.ax.ylabel(self.ylabel)
 69        except AttributeError:
 70            self.ax.set_xscale("log", **base_arg)
 71            self.ax.set_xlabel(self.xlabel)
 72            self.ax.set_ylabel(self.ylabel)
 73        if self.legend:
 74            self.ax.legend()
 75
 76    def show(self):
 77        self.finish()
 78        pyplot.tight_layout()
 79        show_plot()
 80
 81    def save(self, file):
 82        self.finish()
 83        pyplot.savefig(file, format=self.format)
 84
 85    def plot(
 86        self, spectrum, t, label=None, color=None, title=None, add_error_to_label=False
 87    ):
 88        error = self.plot_analytic_solution(self.settings, t, spectrum, title)
 89        if label is not None and add_error_to_label:
 90            label += f" error={error:.4g}"
 91        self.plot_data(self.settings, t, spectrum, label, color)
 92        return error
 93
 94    def plot_analytic_solution(self, settings, t, spectrum, title):
 95        if t == 0:
 96            analytic_solution = settings.spectrum.size_distribution
 97        else:
 98
 99            def analytic_solution(x):
100                return settings.norm_factor * settings.kernel.analytic_solution(
101                    x=x, t=t, x_0=settings.X0, N_0=settings.n_part
102                )
103
104        volume_bins_edges = self.settings.formulae.trivia.volume(
105            settings.radius_bins_edges
106        )
107        dm = np.diff(volume_bins_edges)
108        dr = np.diff(settings.radius_bins_edges)
109
110        pdf_m_x = volume_bins_edges[:-1] + dm / 2
111        pdf_m_y = analytic_solution(pdf_m_x)
112
113        pdf_r_x = settings.radius_bins_edges[:-1] + dr / 2
114        pdf_r_y = pdf_m_y * dm / dr * pdf_r_x
115
116        x = pdf_r_x * si.metres / si.micrometres
117        y_true = (
118            pdf_r_y
119            * self.settings.formulae.trivia.volume(radius=pdf_r_x)
120            * settings.rho
121            / settings.dv
122            * si.kilograms
123            / si.grams
124        )
125
126        self.ax.plot(x, y_true, color="black")
127
128        if spectrum is not None:
129            y = spectrum * si.kilograms / si.grams
130            error = error_measure(y, y_true, x)
131            self.title = (
132                title or f"error measure: {error:.2f}"
133            )  # TODO #327 relative error
134            return error
135        return None
136
137    def plot_data(self, settings, t, spectrum, label, color):
138        if self.smooth:
139            scope = self.smooth_scope
140            if t != 0:
141                new = np.copy(spectrum)
142                for _ in range(2):
143                    for i in range(scope, len(spectrum) - scope):
144                        new[i] = np.mean(spectrum[i - scope : i + scope + 1])
145                    scope = 1
146                    for i in range(scope, len(spectrum) - scope):
147                        spectrum[i] = np.mean(new[i - scope : i + scope + 1])
148
149            x = settings.radius_bins_edges[:-scope]
150            dx = np.diff(x)
151            self.ax.plot(
152                (x[:-1] + dx / 2) * si.metres / si.micrometres,
153                spectrum[:-scope] * si.kilograms / si.grams,
154                label=label or f"t = {t}s",
155                color=color
156                or self.colors(t / (self.settings.output_steps[-1] * self.settings.dt)),
157            )
158        else:
159            self.ax.step(
160                settings.radius_bins_edges[:-1] * si.metres / si.micrometres,
161                spectrum * si.kilograms / si.grams,
162                where="post",
163                label=label or f"t = {t}s",
164                color=color
165                or self.colors(t / (self.settings.output_steps[-1] * self.settings.dt)),
166            )
SpectrumPlotter(settings, title=None, grid=True, legend=True, log_base=10)
31    def __init__(self, settings, title=None, grid=True, legend=True, log_base=10):
32        self.settings = settings
33        self.format = "pdf"
34        self.colors = SpectrumColors()
35        self.smooth = False
36        self.smooth_scope = 2
37        self.legend = legend
38        self.grid = grid
39        self.title = title
40        self.xlabel = "particle radius [µm]"
41        self.ylabel = "dm/dlnr [g/m^3/(unit dr/r)]"
42        self.log_base = log_base
43        self.ax = pyplot
44        self.fig = pyplot
45        self.finished = False
settings
format
colors
smooth
smooth_scope
legend
grid
title
xlabel
ylabel
log_base
ax
fig
finished
def finish(self):
47    def finish(self):
48        if self.finished:
49            return
50        self.finished = True
51        if self.grid:
52            self.ax.grid()
53
54        base_arg = {
55            "base"
56            + (
57                "x" if _matplotlib_version_actual < _matplotlib_version_3_3_3 else ""
58            ): self.log_base
59        }
60        if self.title is not None:
61            try:
62                self.ax.title(self.title)
63            except TypeError:
64                self.ax.set_title(self.title)
65        try:
66            self.ax.xscale("log", **base_arg)
67            self.ax.xlabel(self.xlabel)
68            self.ax.ylabel(self.ylabel)
69        except AttributeError:
70            self.ax.set_xscale("log", **base_arg)
71            self.ax.set_xlabel(self.xlabel)
72            self.ax.set_ylabel(self.ylabel)
73        if self.legend:
74            self.ax.legend()
def show(self):
76    def show(self):
77        self.finish()
78        pyplot.tight_layout()
79        show_plot()
def save(self, file):
81    def save(self, file):
82        self.finish()
83        pyplot.savefig(file, format=self.format)
def plot( self, spectrum, t, label=None, color=None, title=None, add_error_to_label=False):
85    def plot(
86        self, spectrum, t, label=None, color=None, title=None, add_error_to_label=False
87    ):
88        error = self.plot_analytic_solution(self.settings, t, spectrum, title)
89        if label is not None and add_error_to_label:
90            label += f" error={error:.4g}"
91        self.plot_data(self.settings, t, spectrum, label, color)
92        return error
def plot_analytic_solution(self, settings, t, spectrum, title):
 94    def plot_analytic_solution(self, settings, t, spectrum, title):
 95        if t == 0:
 96            analytic_solution = settings.spectrum.size_distribution
 97        else:
 98
 99            def analytic_solution(x):
100                return settings.norm_factor * settings.kernel.analytic_solution(
101                    x=x, t=t, x_0=settings.X0, N_0=settings.n_part
102                )
103
104        volume_bins_edges = self.settings.formulae.trivia.volume(
105            settings.radius_bins_edges
106        )
107        dm = np.diff(volume_bins_edges)
108        dr = np.diff(settings.radius_bins_edges)
109
110        pdf_m_x = volume_bins_edges[:-1] + dm / 2
111        pdf_m_y = analytic_solution(pdf_m_x)
112
113        pdf_r_x = settings.radius_bins_edges[:-1] + dr / 2
114        pdf_r_y = pdf_m_y * dm / dr * pdf_r_x
115
116        x = pdf_r_x * si.metres / si.micrometres
117        y_true = (
118            pdf_r_y
119            * self.settings.formulae.trivia.volume(radius=pdf_r_x)
120            * settings.rho
121            / settings.dv
122            * si.kilograms
123            / si.grams
124        )
125
126        self.ax.plot(x, y_true, color="black")
127
128        if spectrum is not None:
129            y = spectrum * si.kilograms / si.grams
130            error = error_measure(y, y_true, x)
131            self.title = (
132                title or f"error measure: {error:.2f}"
133            )  # TODO #327 relative error
134            return error
135        return None
def plot_data(self, settings, t, spectrum, label, color):
137    def plot_data(self, settings, t, spectrum, label, color):
138        if self.smooth:
139            scope = self.smooth_scope
140            if t != 0:
141                new = np.copy(spectrum)
142                for _ in range(2):
143                    for i in range(scope, len(spectrum) - scope):
144                        new[i] = np.mean(spectrum[i - scope : i + scope + 1])
145                    scope = 1
146                    for i in range(scope, len(spectrum) - scope):
147                        spectrum[i] = np.mean(new[i - scope : i + scope + 1])
148
149            x = settings.radius_bins_edges[:-scope]
150            dx = np.diff(x)
151            self.ax.plot(
152                (x[:-1] + dx / 2) * si.metres / si.micrometres,
153                spectrum[:-scope] * si.kilograms / si.grams,
154                label=label or f"t = {t}s",
155                color=color
156                or self.colors(t / (self.settings.output_steps[-1] * self.settings.dt)),
157            )
158        else:
159            self.ax.step(
160                settings.radius_bins_edges[:-1] * si.metres / si.micrometres,
161                spectrum * si.kilograms / si.grams,
162                where="post",
163                label=label or f"t = {t}s",
164                color=color
165                or self.colors(t / (self.settings.output_steps[-1] * self.settings.dt)),
166            )