8053365388

Committed 26 Feb 2024 05:51PM UTC coverage: 62.974% (-0.3%) from 63.244%

Build # 8053365388

Build Type

push

github

Committed by

ghiggi

Commit Message

Drop support for python 3.8

Run Details

3723 of 5912 relevant lines covered (62.97%)

0.63 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

0.0

/gpm_api/visualization/plot_3d.py

# -----------------------------------------------------------------------------.
# MIT License

# Copyright (c) 2024 GPM-API developers
#
# This file is part of GPM-API.

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

# -----------------------------------------------------------------------------.
"""This module contains functions for 3D visualization of GPM-API RADAR data."""
import numpy as np

# TODO:
# - Isosurface contour buggy at low reflectivity
#   --> Should I replace values as 0-1 at each round?
# - 3D terrain
# - surface on bin surface height


def create_pyvista_2d_surface(data_array, spacing=(1, 1, 1), origin=(0, 0, 0)):
    """Create pyvista ImageData object from 2D xr.DataArray."""
    import pyvista as pv

    dimensions = (data_array.shape[0], data_array.shape[1], 1)
    surf = pv.ImageData(
        dimensions=dimensions,
        spacing=spacing,
        origin=origin,
    )
    data = data_array.values
    data = data[:, ::-1]
    surf.point_data.set_array(data.flatten(order="F"), name=data_array.name)
    surf.set_active_scalars(data_array.name)
    return surf


def create_pyvista_3d_volume(data_array, spacing=(1, 1, 0.25), origin=(0, 0, 0)):
    """Create pyvista ImageData object from 3D xr.DataArray."""
    import pyvista as pv

    # Remove vertical areas without values
    data_array = data_array.gpm_api.slice_range_with_valid_data()

    # Create ImageData object
    # - TODO: scale (factor)
    dimensions = data_array.shape
    vol = pv.ImageData(
        dimensions=dimensions,
        spacing=spacing,
        origin=origin,
    )
    data = data_array.values
    data = data[:, ::-1, ::-1]
    vol.point_data.set_array(data.flatten(order="F"), name=data_array.name)
    vol.set_active_scalars(data_array.name)
    return vol


def get_slider_button_positions(i, num_buttons, spacing_factor=0.04):
    """Return the pointa and pointb parameters for pl.add_slider_widget."""
    if num_buttons < 1:
        raise ValueError("Number of buttons must be at least 1")

    # Define margin
    min_pointa = 0.025
    max_pointb = 0.98

    # Define allowable buttons width
    total_width = max_pointb - min_pointa - spacing_factor * (num_buttons - 1)

    # Define button width
    button_width = total_width / num_buttons

    # Define pointa and pointb
    start_x = min_pointa + i * (button_width)
    if i > 0:
        start_x = start_x + spacing_factor
    end_x = start_x + button_width
    pointa = (start_x, 0.1)
    pointb = (end_x, 0.1)
    return pointa, pointb


class OpacitySlider:
    """Opacity Slider for pyvista pl.add_slider_widget."""

    def __init__(self, pl_actor):
        self.pl_actor = pl_actor

    def __call__(self, value):
        self.pl_actor.GetProperty().SetOpacity(value)


def add_3d_isosurfaces(
    vol,
    pl,
    isovalues=[30, 40, 50],
    opacities=[0.3, 0.5, 1],
    method="contour",
    style="surface",
    add_sliders=False,
    **mesh_kwargs,
):
    """Add 3D isosurface to a pyvista plotter object.

    If add_sliders=True, isosurface opacity can be adjusted interactively.

    """
    # Checks
    if len(isovalues) != len(opacities):
        raise ValueError("Expected same number of isovalues and opacities values.")
    # TODO: check there are values larger than max isovalues

    # Define opacity dictionary
    dict_opacity = dict(zip(isovalues, opacities))

    # Precompute isosurface
    dict_isosurface = {isovalue: vol.contour([isovalue], method=method) for isovalue in isovalues}
    n_isosurfaces = len(dict_isosurface)
    # Add isosurfaces
    for i, (isovalue, isosurface) in enumerate(dict_isosurface.items()):
        pl_actor = pl.add_mesh(
            isosurface, opacity=dict_opacity[isovalue], style=style, **mesh_kwargs
        )
        if add_sliders:
            # Define opacity slider
            opacity_slider = OpacitySlider(pl_actor)
            # Define slicer button position
            pointa, pointb = get_slider_button_positions(i=i, num_buttons=n_isosurfaces)
            # Add slider widget
            pl.add_slider_widget(
                callback=opacity_slider,
                rng=[0, 1],
                value=dict_opacity[isovalue],
                title=f"Isovalue={isovalue}",
                pointa=pointa,
                pointb=pointb,
                style="modern",
            )


class IsosurfaceSlider:
    def __init__(self, vol, method="contour", isovalue=None):
        """Define pyvista slider for 3D isosurfaces."""

        self.vol = vol
        vmin, vmax = vol.get_data_range()
        self.vmin = vmin
        self.vmax = vmax
        self.method = method
        if isovalue is None:
            isovalue = vmin + (vmax - vmin) / 2
        self.isovalue = isovalue
        self.isosurface = vol.contour([isovalue], method=method)

    def __call__(self, value):
        self.isovalue = value
        self.update()

    def update(self):
        result = self.vol.contour([self.isovalue], method=self.method)
        self.isosurface.copy_from(result)


def add_3d_isosurface_slider(vol, pl, method="contour", isovalue=None, **mesh_kwargs):
    """Add a 3D isosurface slider enabling to slide through the 3D volume."""
    isosurface_slider = IsosurfaceSlider(vol, method=method, isovalue=isovalue)
    isosurface = isosurface_slider.isosurface
    pl.add_mesh(isosurface, **mesh_kwargs)
    pl.add_slider_widget(
        callback=isosurface_slider,
        rng=vol.get_data_range(),
        value=isosurface_slider.isovalue,
        title="Isosurface",
        pointa=(0.4, 0.9),
        pointb=(0.9, 0.9),
        style="modern",
    )


class OrthogonalSlicesSlider:
    def __init__(self, vol, x=1, y=1, z=1):
        """Define pyvista sliders for 3D orthogonal slices."""

        self.vol = vol
        self.slices = vol.slice_orthogonal(x=x, y=y, z=z)
        # Set default parameters
        self.kwargs = {
            "x": x,
            "y": y,
            "z": z,
        }

    def __call__(self, param, value):
        self.kwargs[param] = value
        self.update()

    def update(self):
        # This is where you call your simulation
        result = self.vol.slice_orthogonal(**self.kwargs)
        self.slices[0].copy_from(result[0])
        self.slices[1].copy_from(result[1])
        self.slices[2].copy_from(result[2])
        return


def add_3d_orthogonal_slices(vol, pl, x=None, y=None, z=None, add_sliders=False, **mesh_kwargs):
    """Add 3D orthogonal slices with interactive sliders."""
    # Define bounds
    x_rng = vol.bounds[0:2]
    y_rng = vol.bounds[2:4]
    z_rng = vol.bounds[4:6]

    # Define default values if not provided
    # - If value is 0, means no slice plotted !
    if x is None:
        x = int(np.diff(x_rng) / 2)
    if y is None:
        y = int(np.diff(y_rng) / 2)
    if z is None:
        z = int(z_rng[0] + 0.01)

    # Define orthogonal slices (and sliders)
    if add_sliders:
        orthogonal_slices_slider = OrthogonalSlicesSlider(vol)
        orthogonal_slices = orthogonal_slices_slider.slices
    else:
        orthogonal_slices = vol.slice_orthogonal(x=x, y=y, z=z, progress_bar=False)

    # Display orthogonal slices
    pl.add_mesh(orthogonal_slices, **mesh_kwargs)

    # Add slider widgets
    if add_sliders:
        pl.add_slider_widget(
            callback=lambda value: orthogonal_slices_slider("x", int(value)),
            rng=x_rng,
            value=x,
            title="Along-Track",
            pointa=(0.025, 0.1),
            pointb=(0.31, 0.1),
            style="modern",
        )
        pl.add_slider_widget(
            callback=lambda value: orthogonal_slices_slider("y", int(value)),
            rng=y_rng,
            value=y,
            title="Cross-Track",
            pointa=(0.35, 0.1),
            pointb=(0.64, 0.1),
            style="modern",
        )
        pl.add_slider_widget(
            callback=lambda value: orthogonal_slices_slider("z", value),
            rng=z_rng,
            value=z,
            title="Elevation",
            pointa=(0.67, 0.1),
            pointb=(0.98, 0.1),
            style="modern",
        )

1	# -----------------------------------------------------------------------------.
2	# MIT License
3
4	# Copyright (c) 2024 GPM-API developers
5	#
6	# This file is part of GPM-API.
7
8	# Permission is hereby granted, free of charge, to any person obtaining a copy
9	# of this software and associated documentation files (the "Software"), to deal
10	# in the Software without restriction, including without limitation the rights
11	# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
12	# copies of the Software, and to permit persons to whom the Software is
13	# furnished to do so, subject to the following conditions:
14	#
15	# The above copyright notice and this permission notice shall be included in all
16	# copies or substantial portions of the Software.
17	#
18	# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
19	# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
20	# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
21	# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
22	# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
23	# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
24	# SOFTWARE.
25
26	# -----------------------------------------------------------------------------.
27	"""This module contains functions for 3D visualization of GPM-API RADAR data."""	×
28	import numpy as np	×
29
30	# TODO:
31	# - Isosurface contour buggy at low reflectivity
32	# --> Should I replace values as 0-1 at each round?
33	# - 3D terrain
34	# - surface on bin surface height
35
36
37	def create_pyvista_2d_surface(data_array, spacing=(1, 1, 1), origin=(0, 0, 0)):	×
38	"""Create pyvista ImageData object from 2D xr.DataArray."""
39	import pyvista as pv	×
40
41	dimensions = (data_array.shape[0], data_array.shape[1], 1)	×
42	surf = pv.ImageData(	×
43	dimensions=dimensions,
44	spacing=spacing,
45	origin=origin,
46	)
47	data = data_array.values	×
48	data = data[:, ::-1]	×
49	surf.point_data.set_array(data.flatten(order="F"), name=data_array.name)	×
50	surf.set_active_scalars(data_array.name)	×
51	return surf	×
52
53
54	def create_pyvista_3d_volume(data_array, spacing=(1, 1, 0.25), origin=(0, 0, 0)):	×
55	"""Create pyvista ImageData object from 3D xr.DataArray."""
56	import pyvista as pv	×
57
58	# Remove vertical areas without values
59	data_array = data_array.gpm_api.slice_range_with_valid_data()	×
60
61	# Create ImageData object
62	# - TODO: scale (factor)
63	dimensions = data_array.shape	×
64	vol = pv.ImageData(	×
65	dimensions=dimensions,
66	spacing=spacing,
67	origin=origin,
68	)
69	data = data_array.values	×
70	data = data[:, ::-1, ::-1]	×
71	vol.point_data.set_array(data.flatten(order="F"), name=data_array.name)	×
72	vol.set_active_scalars(data_array.name)	×
73	return vol	×
74
75
76	def get_slider_button_positions(i, num_buttons, spacing_factor=0.04):	×
77	"""Return the pointa and pointb parameters for pl.add_slider_widget."""
78	if num_buttons < 1:	×
79	raise ValueError("Number of buttons must be at least 1")	×
80
81	# Define margin
82	min_pointa = 0.025	×
83	max_pointb = 0.98	×
84
85	# Define allowable buttons width
86	total_width = max_pointb - min_pointa - spacing_factor * (num_buttons - 1)	×
87
88	# Define button width
89	button_width = total_width / num_buttons	×
90
91	# Define pointa and pointb
92	start_x = min_pointa + i * (button_width)	×
93	if i > 0:	×
94	start_x = start_x + spacing_factor	×
95	end_x = start_x + button_width	×
96	pointa = (start_x, 0.1)	×
97	pointb = (end_x, 0.1)	×
98	return pointa, pointb	×
99
100
101	class OpacitySlider:	×
102	"""Opacity Slider for pyvista pl.add_slider_widget."""
103
104	def __init__(self, pl_actor):	×
105	self.pl_actor = pl_actor	×
106
107	def __call__(self, value):	×
108	self.pl_actor.GetProperty().SetOpacity(value)	×
109
110
111	def add_3d_isosurfaces(	×
112	vol,
113	pl,
114	isovalues=[30, 40, 50],
115	opacities=[0.3, 0.5, 1],
116	method="contour",
117	style="surface",
118	add_sliders=False,
119	**mesh_kwargs,
120	):
121	"""Add 3D isosurface to a pyvista plotter object.
122
123	If add_sliders=True, isosurface opacity can be adjusted interactively.
124
125	"""
126	# Checks
127	if len(isovalues) != len(opacities):	×
128	raise ValueError("Expected same number of isovalues and opacities values.")	×
129	# TODO: check there are values larger than max isovalues
130
131	# Define opacity dictionary
132	dict_opacity = dict(zip(isovalues, opacities))	×
133
134	# Precompute isosurface
135	dict_isosurface = {isovalue: vol.contour([isovalue], method=method) for isovalue in isovalues}	×
136	n_isosurfaces = len(dict_isosurface)	×
137	# Add isosurfaces
138	for i, (isovalue, isosurface) in enumerate(dict_isosurface.items()):	×
139	pl_actor = pl.add_mesh(	×
140	isosurface, opacity=dict_opacity[isovalue], style=style, **mesh_kwargs
141	)
142	if add_sliders:	×
143	# Define opacity slider
144	opacity_slider = OpacitySlider(pl_actor)	×
145	# Define slicer button position
146	pointa, pointb = get_slider_button_positions(i=i, num_buttons=n_isosurfaces)	×
147	# Add slider widget
148	pl.add_slider_widget(	×
149	callback=opacity_slider,
150	rng=[0, 1],
151	value=dict_opacity[isovalue],
152	title=f"Isovalue={isovalue}",
153	pointa=pointa,
154	pointb=pointb,
155	style="modern",
156	)
157
158
159	class IsosurfaceSlider:	×
160	def __init__(self, vol, method="contour", isovalue=None):	×
161	"""Define pyvista slider for 3D isosurfaces."""
162
163	self.vol = vol	×
164	vmin, vmax = vol.get_data_range()	×
165	self.vmin = vmin	×
166	self.vmax = vmax	×
167	self.method = method	×
168	if isovalue is None:	×
169	isovalue = vmin + (vmax - vmin) / 2	×
170	self.isovalue = isovalue	×
171	self.isosurface = vol.contour([isovalue], method=method)	×
172
173	def __call__(self, value):	×
174	self.isovalue = value	×
175	self.update()	×
176
177	def update(self):	×
178	result = self.vol.contour([self.isovalue], method=self.method)	×
179	self.isosurface.copy_from(result)	×
180
181
182	def add_3d_isosurface_slider(vol, pl, method="contour", isovalue=None, **mesh_kwargs):	×
183	"""Add a 3D isosurface slider enabling to slide through the 3D volume."""
184	isosurface_slider = IsosurfaceSlider(vol, method=method, isovalue=isovalue)	×
185	isosurface = isosurface_slider.isosurface	×
186	pl.add_mesh(isosurface, **mesh_kwargs)	×
187	pl.add_slider_widget(	×
188	callback=isosurface_slider,
189	rng=vol.get_data_range(),
190	value=isosurface_slider.isovalue,
191	title="Isosurface",
192	pointa=(0.4, 0.9),
193	pointb=(0.9, 0.9),
194	style="modern",
195	)
196
197
198	class OrthogonalSlicesSlider:	×
199	def __init__(self, vol, x=1, y=1, z=1):	×
200	"""Define pyvista sliders for 3D orthogonal slices."""
201
202	self.vol = vol	×
203	self.slices = vol.slice_orthogonal(x=x, y=y, z=z)	×
204	# Set default parameters
205	self.kwargs = {	×
206	"x": x,
207	"y": y,
208	"z": z,
209	}
210
211	def __call__(self, param, value):	×
212	self.kwargs[param] = value	×
213	self.update()	×
214
215	def update(self):	×
216	# This is where you call your simulation
217	result = self.vol.slice_orthogonal(**self.kwargs)	×
218	self.slices[0].copy_from(result[0])	×
219	self.slices[1].copy_from(result[1])	×
220	self.slices[2].copy_from(result[2])	×
221	return	×
222
223
224	def add_3d_orthogonal_slices(vol, pl, x=None, y=None, z=None, add_sliders=False, **mesh_kwargs):	×
225	"""Add 3D orthogonal slices with interactive sliders."""
226	# Define bounds
227	x_rng = vol.bounds[0:2]	×
228	y_rng = vol.bounds[2:4]	×
229	z_rng = vol.bounds[4:6]	×
230
231	# Define default values if not provided
232	# - If value is 0, means no slice plotted !
233	if x is None:	×
234	x = int(np.diff(x_rng) / 2)	×
235	if y is None:	×
236	y = int(np.diff(y_rng) / 2)	×
237	if z is None:	×
238	z = int(z_rng[0] + 0.01)	×
239
240	# Define orthogonal slices (and sliders)
241	if add_sliders:	×
242	orthogonal_slices_slider = OrthogonalSlicesSlider(vol)	×
243	orthogonal_slices = orthogonal_slices_slider.slices	×
244	else:
245	orthogonal_slices = vol.slice_orthogonal(x=x, y=y, z=z, progress_bar=False)	×
246
247	# Display orthogonal slices
248	pl.add_mesh(orthogonal_slices, **mesh_kwargs)	×
249
250	# Add slider widgets
251	if add_sliders:	×
252	pl.add_slider_widget(	×
253	callback=lambda value: orthogonal_slices_slider("x", int(value)),
254	rng=x_rng,
255	value=x,
256	title="Along-Track",
257	pointa=(0.025, 0.1),
258	pointb=(0.31, 0.1),
259	style="modern",
260	)
261	pl.add_slider_widget(	×
262	callback=lambda value: orthogonal_slices_slider("y", int(value)),
263	rng=y_rng,
264	value=y,
265	title="Cross-Track",
266	pointa=(0.35, 0.1),
267	pointb=(0.64, 0.1),
268	style="modern",
269	)
270	pl.add_slider_widget(	×
271	callback=lambda value: orthogonal_slices_slider("z", value),
272	rng=z_rng,
273	value=z,
274	title="Elevation",
275	pointa=(0.67, 0.1),
276	pointb=(0.98, 0.1),
277	style="modern",
278	)

ghiggi / gpm_api / 8053365388

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous