CHEC camera geometry converter classes broken for Prod4

[STSpencer]

The chec_to_2d_array() camera geometry class no longer works with CHEC Prod4 simulations. This is what you get for a muon ring if you try to use it:

The cameradisplay class works fine though using event.inst.subarray.tel[tel_id].camera

Is there a workaround to extract a 2d array as an image from the cameradisplay class? I've tried using the new geometry with the old 2d array extractor to no avail.

[kosack]

Perhaps we should consider replacing all of that functionality with the ones implemented by the ML group? Isn't there now a general library for doing it using multiple methods?

[kosack]

I guess the problem is just that the pixel ordering has changed, so the method used now is far too simplistic, and probably needs to be replaced with a version that doesn't assume the pixels are in some specific order (similar to how the hex1d_to_rect2d method).

Not sure the best way to proceed - we don't want to break reading older prods for now, and it's not easy to detect which prod is being used.

I guess the problem is just that the pixel ordering has changed, so the method used now is far too simplistic, and probably needs to be replaced with a version that doesn't assume the pixels are in some specific order (similar to how the hex1d_to_rect2d method).

[kosack]

currently the pixel transform is hardcoded (but can be replaced with another):

def chec_transformation_map():
    # By default, pixels map to the last element of input_img_ext (i.e. NaN)
    img_map = np.full([8 * 6, 8 * 6], -1, dtype=int)

    # Map values
    img_map[:8, 8:-8] = np.arange(8 * 8 * 4).reshape([8, 8 * 4])
    img_map[8:40, :] = np.arange(32 * 48).reshape([32, 48]) + 256
    img_map[-8:, 8:-8] = np.arange(8 * 8 * 4).reshape([8, 8 * 4]) + 1792

    return img_map

You specify the transformation when calling it - so you could define a new one for Prod4 (again not sure if there is a good way to do that automatically)

chec_to_2d_array(input_img, img_map=chec_transformation_map())

[maxnoe]

IMHO we should start thinking about requiring pixel positions in additional coordinate systems. E.g. rows and cols for square pixels and axial coordinates (see https://www.redblobgames.com/grids/hexagons/#coordinates) for the hexagonal ones.

Or a way to correctly calculate those in ctapipe. (I have python/numpy code to go from cartesian to axial coordinates here https://github.com/maxnoe/pyhexgrid)

[maxnoe]

What do you think about this:

from ctapipe.instrument import CameraGeometry
from ctapipe.visualization import CameraDisplay
from ctapipe.image.toymodel import Gaussian
import matplotlib.pyplot as plt
import astropy.units as u
import numpy as np

cam = CameraGeometry.from_name('CHEC')
model = Gaussian(0.05 * u.m, 0.03 * u.m, 0.05 * u.m, 0.02*u.m, '30d')
img, _, _ = model.generate_image(cam, 2000, 3)

size = np.sqrt(cam.pix_area)

def pos_to_index(pos, size):
    rnd = np.round((pos / size).to_value(u.dimensionless_unscaled), 1)
    unique = np.sort(np.unique(rnd))
    mask = np.append(np.diff(unique) > 0.5, True)
    bins = np.append(unique[mask] - 0.5, unique[-1] + 0.5)
    return np.digitize(rnd, bins) - 1

col = pos_to_index(cam.pix_x, size)
row = pos_to_index(cam.pix_y, size)

img_sq = np.full((row.max() + 1, col.max() + 1), np.nan)
img_sq[row, col] = img

fig, (ax1, ax2) = plt.subplots(1, 2)

disp = CameraDisplay(cam, image=img, ax=ax1)
disp.add_colorbar(ax=ax1)

img = ax2.imshow(img_sq[::-1])

fig.colorbar(img, ax=ax2)

plt.show()

Basically binning the pixels on a grid where the bin width is at least half the pixel size.

[watsonjj]

I have a similar approach implemented here: https://github.com/watsonjj/CHECLabPy/blob/2eb4a1845df68c9f35ea378ad3fb8fab1294c290/CHECLabPy/utils/mapping.py#L176

[STSpencer]

Thanks @MaxNoe and @watsonjj ; that's done the trick for now.

[kosack]

that's also not far from how the hex binning is done, I think, so looks good. Perhaps open a PR?

[maxnoe]

This was fixed in #1728