Antenna reader example

[1]:
from radiocalibrationtoolkit import *
from healpy.newvisufunc import projview
import plotly.graph_objects as go

# This ensures Plotly output works in multiple places:
# plotly_mimetype: VS Code notebook UI
# notebook: "Jupyter: Export to HTML" command in VS Code
# See https://plotly.com/python/renderers/#multiple-renderers
import plotly.io as pio

pio.renderers.default = "plotly_mimetype+notebook"
[INFO] LFmap: Import successful.
[2]:
layout_settings = dict(
    xaxis=dict(title="azimuth [°]",
               tickprefix="<b>",
               ticksuffix="</b>",
               dtick=30),
    yaxis=dict(
        title="zenith angle [°]",
        tickprefix="<b>",
        ticksuffix="</b>",
        dtick=10,
    ),
    coloraxis=dict(colorbar=dict(
        title=dict(
            text="VEL",
            side="right",
        ),
        tickprefix="<b>",
        ticksuffix="</b>",
    ), ),
    font=dict(
        size=15,
        color="black",
    ),
)

colormap = 'jet'


def get_plotly_cbar_label(quantity):
    if quantity == 'absolute':
        return "<|H|>"
    elif "Phi_amp" in quantity:
        return "|H<sub>Φ</sub>|"
    elif "Theta_amp" in quantity:
        return "|H<sub>ϴ</sub>|"


shift_phi = 0
[3]:
# create antenna instance
antenna_inst = AntennaPattern("./antenna_setup_files/SALLA_EW.xml")

# the NS antenna pattern is shift by 90 degrees clockwise in azimuth, to turn in back, we use shift_phi=-90
# shift_phi = -90
# new_antenna_conventions = {
#     "shift_phi": shift_phi,
# }
# antenna_inst = AntennaPattern("./antenna_setup_files/SALLA_NS.xml",
#                               new_antenna_conventions=new_antenna_conventions)
[INFO] Your keys are: ['EAHTheta_amp', 'EAHTheta_phase', 'EAHPhi_amp', 'EAHPhi_phase', 'absolute'] Use them as the 'quantity'
[4]:
# XML values in dictionary
# antenna_inst.get_raw()
[5]:
# get antenna gain for unpolarized emission
df = antenna_inst.get(frequency=45, quantity="absolute")

# plot
fig = px.imshow(df.T,
                width=600,
                aspect="cube",
                origin='lower',
                color_continuous_scale=colormap)
fig.update_layout(**layout_settings)
fig.update_layout(coloraxis=dict(colorbar=dict(title=dict(
    text=get_plotly_cbar_label("absolute")))))
fig.show()
[6]:
# interpolate
df = antenna_inst.get(frequency=45,
                      quantity='absolute',
                      interp_phi=np.linspace(0, 360, 200),
                      interp_theta=np.linspace(0, 90, 100))

# plot
fig = px.imshow(df.T,
                width=600,
                aspect="cube",
                origin='lower',
                color_continuous_scale=colormap)
fig.update_layout(**layout_settings)
fig.update_layout(coloraxis=dict(colorbar=dict(title=dict(
    text=get_plotly_cbar_label("absolute")))))
fig.show()
[7]:
# antenna gain as a function of frequency, azimuth and zenith in the required polarization as volumetric dataset
volumetric_dataset_df = antenna_inst.get_volumetric_dataset(
    quantity="EAHTheta_amp", frequencies=np.arange(30, 80, 1))
volumetric_dataset_df
[7]:
theta 0.0 1.5 3.0 4.5 6.0 7.5 9.0 10.5 12.0 13.5 ... 76.5 78.0 79.5 81.0 82.5 84.0 85.5 87.0 88.5 90.0
freq phi
30 0.0 1.60419 1.60520 1.60644 1.60790 1.60954 1.61132 1.61319 1.61509 1.61695 1.61869 ... 1.04430 0.99178 0.92724 0.84915 0.75593 0.64597 0.51750 0.36859 0.19698 0.0
3.0 1.60425 1.60522 1.60643 1.60786 1.60948 1.61124 1.61309 1.61497 1.61682 1.61855 ... 1.04482 0.99226 0.92768 0.84954 0.75627 0.64625 0.51773 0.36874 0.19706 0.0
6.0 1.59993 1.60086 1.60205 1.60345 1.60504 1.60678 1.60861 1.61048 1.61231 1.61402 ... 1.04277 0.99030 0.92584 0.84784 0.75475 0.64495 0.51667 0.36799 0.19666 0.0
9.0 1.59122 1.59213 1.59329 1.59467 1.59623 1.59795 1.59977 1.60162 1.60344 1.60514 ... 1.03815 0.98591 0.92171 0.84405 0.75137 0.64205 0.51435 0.36634 0.19578 0.0
12.0 1.57817 1.57905 1.58019 1.58154 1.58309 1.58478 1.58658 1.58842 1.59022 1.59191 ... 1.03097 0.97908 0.91531 0.83818 0.74614 0.63758 0.51076 0.36378 0.19441 0.0
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
79 348.0 1.45920 1.46010 1.45994 1.45871 1.45636 1.45283 1.44810 1.44215 1.43499 1.42669 ... 1.49450 1.53392 1.54726 1.51992 1.43914 1.29492 1.08041 0.79232 0.43093 0.0
351.0 1.47382 1.47481 1.47469 1.47346 1.47107 1.46746 1.46263 1.45653 1.44920 1.44068 ... 1.51165 1.55085 1.56376 1.53564 1.45367 1.30774 1.09096 0.79998 0.43507 0.0
354.0 1.48440 1.48547 1.48540 1.48419 1.48178 1.47814 1.47323 1.46703 1.45956 1.45090 ... 1.52568 1.56459 1.57702 1.54819 1.46519 1.31785 1.09924 0.80598 0.43831 0.0
357.0 1.49092 1.49206 1.49205 1.49087 1.48847 1.48483 1.47988 1.47363 1.46609 1.45734 ... 1.53642 1.57494 1.58688 1.55740 1.47356 1.32515 1.10518 0.81025 0.44062 0.0
360.0 1.49334 1.49456 1.49462 1.49349 1.49113 1.48750 1.48257 1.47631 1.46875 1.45998 ... 1.54367 1.58174 1.59317 1.56314 1.47866 1.32950 1.10867 0.81275 0.44195 0.0

6050 rows × 61 columns

[8]:
# zenith slice
df = volumetric_dataset_df.unstack().stack(level=0).loc[:, 66, :]
# plot
fig = px.imshow(df,
                width=600,
                aspect="cube",
                origin='lower',
                color_continuous_scale=colormap)
fig.update_layout(**layout_settings)
fig.update_layout(
    title="Slice at 66° zenith angle",
    xaxis=dict(title="azimuth [°]"),
    yaxis=dict(title="frequency [MHz]"),
    coloraxis=dict(colorbar=dict(title=dict(
        text=get_plotly_cbar_label("EAHTheta_amp")))),
)

fig.show()

# azimuth slice
df = volumetric_dataset_df.loc[:, 180, :]
# plot
fig = px.imshow(df,
                width=600,
                aspect="cube",
                origin='lower',
                color_continuous_scale=colormap)
fig.update_layout(**layout_settings)
fig.update_layout(
    title="Slice at 180° azimuth",
    xaxis=dict(title="zenith angle [°]", dtick=10),
    yaxis=dict(title="frequency [MHz]"),
    coloraxis=dict(colorbar=dict(title=dict(
        text=get_plotly_cbar_label("EAHTheta_amp")))),
)
fig.show()
[9]:
# some settings
font = "Arial Black"
colorscale = "jet"

scene = dict(
    xaxis=dict(title="<b>Zenith angle [°]</b>", color="black", dtick="30"),
    yaxis=dict(title="<b>Azimuth [°]</b>", color="black", dtick="90"),
    zaxis=dict(
        title="<b>frequency [MHz]</b>",
        tickvals=[40, 50, 60, 70, 80],
    ),
    aspectratio=dict(x=0.75, y=0.75,
                     z=0.75),  # Adjust the aspect ratio as needed
)

# define some functions


def get_volume_data4plots(antenna_inst, quantity):
    """
    Get volumetric data and grid for creating 3D plots.

    Parameters
    ----------
    antenna_inst : AntennaInstance
        An instance of the antenna used for data retrieval.
    quantity : str
        The quantity of interest (e.g., 'absolute', 'EAHTheta_amp').

    Returns
    -------
    tuple
        A tuple containing:
        1. volumetric_dataset_df : pandas.DataFrame
            The volumetric dataset as a pandas DataFrame.
        2. volume_data : numpy.ndarray
            The reshaped volumetric dataset.
        3. PHI : numpy.ndarray
            The meshgrid of PHI values.
        4. FREQ : numpy.ndarray
            The meshgrid of frequency values.
        5. THETA : numpy.ndarray
            The meshgrid of THETA values.
    """

    volumetric_dataset_df = antenna_inst.get_volumetric_dataset(
        quantity=quantity, frequencies=np.arange(30, 80, 3))

    volume_data = volumetric_dataset_df.to_numpy().reshape(
        volumetric_dataset_df.index.get_level_values("freq").nunique(),
        -1,
        volumetric_dataset_df.shape[1],
    )
    PHI, FREQ, THETA = np.meshgrid(
        volumetric_dataset_df.index.get_level_values("phi").unique().values,
        volumetric_dataset_df.index.get_level_values("freq").unique().values,
        volumetric_dataset_df.columns.values,
    )
    return volumetric_dataset_df, volume_data, PHI, FREQ, THETA


def quantity_label2colorbar_dict(quantity_label):
    """
    Create a colorbar dictionary for a given quantity label.

    Parameters
    ----------
    quantity_label : str
        The label for the quantity to be displayed on the colorbar.

    Returns
    -------
    dict
        A dictionary containing colorbar configuration settings.
    """
    return dict(len=0.75,
                title=dict(
                    text="<b>" + quantity_label + " [m]</b>",
                    side="right",
                ),
                tickprefix="<b>",
                ticksuffix="</b>",
                tickfont=dict(size=30))


def create_volume_trace(
        volume_data,
        PHI,
        FREQ,
        THETA,
        quantity,
        opacity=1,
        surface_fill=0.001,
        surface_count=1,
        opacityscale="uniform",
        caps=dict(x_show=False, y_show=False, z_show=False),
):
    """
    Create a go.Volume trace based on the given parameters.

    Parameters
    ----------
    volume_data : numpy.ndarray
        The reshaped volumetric data.
    PHI : numpy.ndarray
        The meshgrid of PHI values.
    FREQ : numpy.ndarray
        The meshgrid of frequency values.
    THETA : numpy.ndarray
        The meshgrid of THETA values.
    quantity : str
        The quantity of interest.
    opacity : float, optional
        The opacity of the volume. Default is 1.
    surface_fill : float, optional
        The surface fill value. Default is 0.001.
    surface_count : int, optional
        The surface count value. Default is 1.
    opacityscale : str, optional
        The opacity scale. Default is "uniform".
    caps : dict, optional
        The caps configuration. Default is dict(x_show=False, y_show=False, z_show=False).

    Returns
    -------
    go.Volume
        A go.Volume trace object.
    """
    return go.Volume(
        y=PHI.flatten(),
        x=THETA.flatten(),
        z=FREQ.flatten(),
        value=np.around(volume_data.flatten(), 3),
        cmin=0,
        cmax=2,
        opacity=opacity,
        # isomax=1.3,
        surface_fill=surface_fill,
        surface_count=surface_count,
        colorscale="jet",
        opacityscale=opacityscale,
        caps=caps,
        colorbar=quantity_label2colorbar_dict(get_plotly_cbar_label(quantity)),
    )


def built_volume_and_slice_plots(
    antenna_inst,
    quantity='absolute',
    slice_types={
        "no_slices": [],
        "slices_z": [45, 60],
        "slices_x": [30, 65],
        "slices_y": [0, 270],
    }):
    """
    Build 3D volume and slice plots for a given antenna instance.

    Parameters
    ----------
    antenna_inst : AntennaInstance
        An instance of the antenna used for data retrieval.
    quantity : str, optional
        The quantity of interest (default is 'absolute').
    slice_types : dict, optional
        Dictionary specifying slice types and their locations (default is None).
    """

    volumetric_dataset_df, volume_data, PHI, FREQ, THETA = get_volume_data4plots(
        antenna_inst, quantity)

    # make plots
    for slice_type, locations in slice_types.items():
        if slice_type == "no_slices":
            fig = go.Figure()
            fig.add_trace(
                create_volume_trace(
                    volume_data,
                    PHI,
                    FREQ,
                    THETA,
                    quantity,
                    opacity=0.8,
                    surface_fill=1,
                    opacityscale="max",
                    surface_count=40,
                    caps=dict(x_show=True, y_show=True, z_show=True),
                ))
        else:
            fig = go.Figure()
            slice_trace = create_volume_trace(
                volume_data,
                PHI,
                FREQ,
                THETA,
                quantity,
            )
            slice_trace.update(
                **{slice_type: dict(show=True, locations=locations)})
            fig.add_trace(slice_trace)

        fig.update_layout(
            scene=scene,
            margin=dict(r=50, b=10, l=10, t=10),
            height=600,
            autosize=False,
            font=dict(family=font, size=18, color="black"),
        )
        fig.show()
[10]:
## show antenna pattern as volumetric data with slices in all axis
built_volume_and_slice_plots(antenna_inst, quantity="absolute")
# built_volume_and_slice_plots(antenna_inst, quantity="EAHTheta_amp")
# built_volume_and_slice_plots(antenna_inst, quantity="EAHPhi_amp")
[11]:
# show antenna pattern as volumetric data with slices in all axis
# version of plots with sliders
# first some functions...


def show_slices(x, y, z, volume_data, quantity, slice_type="x"):
    """
    Show slices of volumetric data using the given parameters.

    Parameters
    ----------
    x : np.ndarray
        The x-coordinate values.
    y : np.ndarray
        The y-coordinate values.
    z : np.ndarray
        The z-coordinate values.
    volume_data : np.ndarray
        The volumetric data.
    slice_type : str, optional
        The type of slice. Valid values are 'x', 'y', or 'z'. Default is 'x'.

    Returns
    -------
    None
        This function does not return anything. It displays the plot.

    """
    quantity_label = get_plotly_cbar_label(quantity)

    x_delta = np.diff(x)[0]
    y_delta = np.diff(y)[0]
    z_delta = np.diff(z)[0]

    cmin = 0
    cmax = 2
    colorscale = "jet"

    if slice_type == "x":
        slider_label = x
        fig = go.Figure(frames=[
            go.Frame(
                data=go.Surface(
                    x=np.ones(x.size) * x[k],
                    y=y,
                    z=(np.ones((z.size, y.size)) * z[:, np.newaxis]).T,
                    surfacecolor=(volume_data[:, :, k]).T,
                    cmin=cmin,
                    cmax=cmax,
                ),
                name=str(k),
            ) for k in range(x.size)
        ])

        # default fig
        fig.add_trace(
            go.Surface(
                x=np.ones(x.size) * x[0],
                y=y,
                z=(np.ones((z.size, y.size)) * z[:, np.newaxis]).T,
                surfacecolor=(volume_data[:, :, 0]).T,
                colorscale=colorscale,
                cmin=cmin,
                cmax=cmax,
                colorbar=quantity_label2colorbar_dict(quantity_label),
            ))

    elif slice_type == "y":
        slider_label = y
        fig = go.Figure(frames=[
            go.Frame(
                data=go.Surface(
                    x=x,
                    y=np.ones(y.size) * y[k],
                    z=np.ones((z.size, x.size)) * z[:, np.newaxis],
                    surfacecolor=(volume_data[:, y.size - 1 - k, :]),
                    cmin=cmin,
                    cmax=cmax,
                ),
                name=str(k),
            ) for k in range(y.size)
        ])

        # default fig
        fig.add_trace(
            go.Surface(
                x=x,
                y=np.ones(y.size) * y[0],
                z=np.ones((z.size, x.size)) * z[:, np.newaxis],
                surfacecolor=(volume_data[:, 0, :]),
                colorscale=colorscale,
                cmin=cmin,
                cmax=cmax,
                colorbar=quantity_label2colorbar_dict(quantity_label),
            ))
    elif slice_type == "z":
        slider_label = z

        fig = go.Figure(frames=[
            go.Frame(
                data=go.Surface(
                    x=x,
                    y=y,
                    z=z[k] * np.ones((y.size, x.size)),
                    surfacecolor=np.flipud(volume_data[z.size - 1 - k]),
                    cmin=cmin,
                    cmax=cmax,
                ),
                name=str(k),
            ) for k in range(z.size)
        ])

        # default fig
        fig.add_trace(
            go.Surface(
                x=x,
                y=y,
                z=z[0] * np.ones((y.size, x.size)),
                surfacecolor=np.flipud(volume_data[z.size - 1]),
                colorscale=colorscale,
                cmin=cmin,
                cmax=cmax,
                colorbar=quantity_label2colorbar_dict(quantity_label),
            ))

    def frame_args(duration):
        return {
            "frame": {
                "duration": duration
            },
            "mode": "immediate",
            "fromcurrent": True,
            "transition": {
                "duration": duration,
                "easing": "linear"
            },
        }

    sliders = [{
        "pad": {
            "b": 10,
            "t": 60
        },
        "len":
        0.9,
        "x":
        0.1,
        "y":
        0,
        "steps": [{
            "args": [[f.name], frame_args(0)],
            "label": str(slider_label[k]),
            "method": "animate",
        } for k, f in enumerate(fig.frames)],
    }]

    # Layout
    fig.update_layout(
        title="<br>Slices in volumetric data",
        width=800,
        height=800,
        scene=dict(
            yaxis=dict(range=[y[0] - y_delta, y[-1] + y_delta],
                       autorange=False),  # Set the y-axis range from -5 to +5
            zaxis=dict(range=[z[0] - z_delta, z[-1] + z_delta],
                       autorange=False),
            xaxis=dict(range=[x[0] - x_delta, x[-1] + x_delta],
                       autorange=False),
            aspectratio=dict(x=1, y=1, z=1),
        ),
        updatemenus=[{
            "buttons": [
                {
                    "args": [None, frame_args(50)],
                    "label": "&#9654;",  # play symbol
                    "method": "animate",
                },
                {
                    "args": [[None], frame_args(0)],
                    "label": "&#9724;",  # pause symbol
                    "method": "animate",
                },
            ],
            "direction":
            "left",
            "pad": {
                "r": 10,
                "t": 70
            },
            "type":
            "buttons",
            "x":
            0.1,
            "y":
            0,
        }],
        sliders=sliders,
    )

    fig.update_layout(
        scene=scene,
        margin=dict(r=50, b=10, l=10, t=10),
        font=dict(family=font, size=18, color="black"),
    )

    fig.show()


def get_adjustable_plot_slices(antenna_inst, quantity="absolute"):
    """
    Generate adjustable 3D slices for a given antenna instance and quantity.

    Parameters
    ----------
    antenna_inst : AntennaInstance
        An instance of the antenna used for data retrieval.
    quantity : str, optional
        The quantity of interest (default is "absolute").

    Returns
    -------
    None
        This function generates adjustable 3D slices and does not return any values.
    """
    volumetric_dataset_df = antenna_inst.get_volumetric_dataset(
        quantity=quantity, frequencies=np.arange(30, 80, 3))

    # create arrays
    volume_data = volumetric_dataset_df.to_numpy().reshape(
        volumetric_dataset_df.index.get_level_values("freq").nunique(),
        -1,
        volumetric_dataset_df.shape[1],
    )
    x = volumetric_dataset_df.columns.values
    y = volumetric_dataset_df.index.get_level_values("phi").unique().values
    z = volumetric_dataset_df.index.get_level_values("freq").unique().values

    # create slices with sliders
    show_slices(x, y, z, volume_data, quantity, slice_type="x")
    show_slices(x, y, z, volume_data, quantity, slice_type="y")
    show_slices(x, y, z, volume_data, quantity, slice_type="z")
[12]:
get_adjustable_plot_slices(antenna_inst, quantity="absolute")
# get_adjustable_plot_slices(antenna_inst, quantity="EAHTheta_amp")
# get_adjustable_plot_slices(antenna_inst, quantity="EAHPhi_amp")
[13]:
# convert to healpy format
# additional shift if NS antenna
update_antenna_conventions = {
    'shift_phi': -90 + shift_phi,
    'flip_theta': True,
    'flip_phi': False,
    'in_degrees': True,
    'add_invisible_sky': True
}

antenna_hpmap_inst = antenna_inst.convert2hp(frequency=45,
                                             **update_antenna_conventions)
[14]:
# hp maps are by default in galactic coordinates, so we need to specify the LST and latitude of the local observer
lst = 18
LATITUDE = -35.206667
rotation_parameters = create_rotation_parameters(lst, LATITUDE)
rotator = create_rotator(lst, LATITUDE, coord=["G", "C"])
[15]:
antenna_hpmap = antenna_hpmap_inst.get_map(rotator=rotator)
[16]:
# galaxy
projview(
    antenna_hpmap,
    cmap='jet',
    return_only_data=False,
    graticule=True,
    graticule_labels=True,
    title='Galactic coordinates',
    xtick_label_color='w',
    # projection_type='cart'
)

# from galaxy to local
projview(
    antenna_hpmap,
    cmap='jet',
    return_only_data=False,
    coord=['G', 'C'],
    rot=rotation_parameters,
    graticule=True,
    graticule_labels=True,
    title='Local coordinates',
    xtick_label_color='w',
    # projection_type='cart'
)
[16]:
<matplotlib.collections.QuadMesh at 0x7f16d39c5570>
../_images/examples_antenna_reader_example_16_1.png
../_images/examples_antenna_reader_example_16_2.png