Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add EngineeringCameraFrame #1085

Merged
merged 11 commits into from
Jun 17, 2019
Merged

Add EngineeringCameraFrame #1085

merged 11 commits into from
Jun 17, 2019

Conversation

maxnoe
Copy link
Member

@maxnoe maxnoe commented Jun 6, 2019

Question: Do we need a case separation for 1 / 2 mirror telescopes?

@watsonjj
Copy link
Contributor

watsonjj commented Jun 6, 2019

This is connected to the proposal in #968, correct?

@maxnoe
Copy link
Member Author

maxnoe commented Jun 6, 2019

Yes, I will make a follow up PR introducing camera frame definition/selection to the fits files and maybe to the plotting.

The idea is to transform into the internal frame on reading and maybe back on display if wished.

@codecov
Copy link

codecov bot commented Jun 6, 2019

Codecov Report

Merging #1085 into master will increase coverage by 0.07%.
The diff coverage is 98.36%.

Impacted file tree graph

@@            Coverage Diff             @@
##           master    #1085      +/-   ##
==========================================
+ Coverage   84.25%   84.33%   +0.07%     
==========================================
  Files         181      182       +1     
  Lines       10979    11041      +62     
==========================================
+ Hits         9250     9311      +61     
- Misses       1729     1730       +1
Impacted Files Coverage Δ
ctapipe/instrument/tests/test_camera.py 100% <ø> (ø) ⬆️
ctapipe/instrument/camera.py 91.62% <100%> (+0.55%) ⬆️
ctapipe/visualization/mpl_camera.py 80.22% <100%> (+0.34%) ⬆️
ctapipe/coordinates/__init__.py 100% <100%> (ø) ⬆️
...tapipe/coordinates/tests/test_engineering_frame.py 100% <100%> (ø)
ctapipe/coordinates/camera_frame.py 85.93% <95%> (+4.11%) ⬆️

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update c154718...0343142. Read the comment docs.

@maxnoe
Copy link
Member Author

maxnoe commented Jun 6, 2019

camera_frames

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


fig, axs = plt.subplots(2, 2, constrained_layout=True, figsize=(6, 3))

axs = axs.ravel()

model = Gaussian(0 * u.m, 0.1 * u.m, 0.1 * u.m, 0.03 * u.m, 0 * u.deg)

cam = CameraGeometry.from_name('FlashCam')
image, *_ = model.generate_image(cam, 1500)

CameraDisplay(cam, ax=axs[0], image=image)
CameraDisplay(
    cam.transform_to(EngineeringCameraFrame()),
    ax=axs[1],
    image=image,
)

cam = CameraGeometry.from_name('LSTCam')
image, *_ = model.generate_image(cam, 1500)

CameraDisplay(cam, ax=axs[2], image=image)
CameraDisplay(
    cam.transform_to(EngineeringCameraFrame()),
    ax=axs[3],
    image=image,
)

axs[0].set_title('CameraFrame')
axs[1].set_title('EngineeringCameraFrame')


plt.show()

@kosack
Copy link
Contributor

kosack commented Jun 7, 2019

This is great! (been needed for a while due to the confusing Monte-Carlo system)

Be sure to document whether the engineering frame is looking at the back of the camera or the front. Also, perhaps we should draw a "North" arrow on the CameraDisplay, and even maybe some text like "Front" or "Back" (or a vector "into the plane" or "out-of the plane" symbol). That would avoid a lot of confusion.

@kosack
Copy link
Contributor

kosack commented Jun 7, 2019

possibly also related to #255

@kosack
Copy link
Contributor

kosack commented Jun 7, 2019

and #875 and #416

perhaps we can close many issues with this

@maxnoe
Copy link
Member Author

maxnoe commented Jun 7, 2019

The question if it's the back or front is important for dual mirror. This is not yet taken into account here!

So we should ask, what we want for this, I guess front?

@watsonjj
Copy link
Contributor

watsonjj commented Jun 7, 2019

In CHEC, our (engineering) mapping is defined as if we are looking at the camera from the front

if self.geom.pix_type.startswith("hex"):
rr = sqrt(aa * 2 / 3 / sqrt(3)) + 2 * PIXEL_EPSILON
r = sqrt(area * 2 / 3 / sqrt(3)) + 2 * PIXEL_EPSILON
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

these double names were just because single-character variable names were forbidden in the style guide since they are difficult to search for in a text editor. But I guess if we use modern IDEs with true syntax support (rather than search/replace), it's ok.

@kosack
Copy link
Contributor

kosack commented Jun 7, 2019

draw a "North" arrow on the CameraDisplay

Actually the more I think of it this is also dependent on the pointing direction of the telescope, so it's not so easy... I guess in the EngineeringFrame we assume the telescope is parked, so "Up" might be more useful?

@maxnoe
Copy link
Member Author

maxnoe commented Jun 7, 2019

I verified this for single mirror telescopes using FACT and pedestals of stars we took some time ago in 2015, the positions are queried from the yale bright star catalogue (V/50 in Vizier):

fact_star_pedestal

from ctapipe.instrument import CameraGeometry
from ctapipe.visualization import CameraDisplay
import matplotlib.pyplot as plt
from ctapipe.coordinates import EngineeringCameraFrame, CameraFrame
import astropy.units as u
import numpy as np
from astroquery.vizier import Vizier
from astropy.coordinates import SkyCoord, AltAz
from fact.instrument.constants import LOCATION, FOCAL_LENGTH_MM
from astropy.time import Time

print('done')

vizier = Vizier(
    catalog='V/50',
    columns=['_RAJ2000', '_DEJ2000', 'Vmag', 'Name']
)
t = vizier.query_region('Aldebaran', radius=2.5 * u.deg, catalog='V/50')[0]
stars = SkyCoord(ra=t['_RAJ2000'], dec=t['_DEJ2000'], frame='icrs')

obstime = Time('2015-11-16T00:35:11')
altaz = AltAz(obstime=obstime, location=LOCATION)
pointing = SkyCoord.from_name('Aldebaran').transform_to(altaz)

camera_frame = CameraFrame(
    location=LOCATION,
    obstime=obstime,
    focal_length=FOCAL_LENGTH_MM * u.mm,
    telescope_pointing=pointing,
)
engineering_frame = EngineeringCameraFrame(
    location=LOCATION,
    obstime=obstime,
    focal_length=FOCAL_LENGTH_MM * u.mm,
    telescope_pointing=pointing,
)

stars_cam = stars.transform_to(camera_frame)
stars_eng = stars.transform_to(engineering_frame)

print(stars_cam.x)
print(stars_eng.x)


fig, axs = plt.subplots(1, 2, constrained_layout=True, figsize=(6, 3))

cam = CameraGeometry.from_name('FACT')
image = np.genfromtxt('./aldebaran_pedvar.txt')

d1 = CameraDisplay(cam, ax=axs[0], image=image, cmap='inferno')
d2 = CameraDisplay(
    cam.transform_to(engineering_frame),
    ax=axs[1],
    image=image,
    cmap='inferno',
)

axs[0].plot(stars_cam.x.value, stars_cam.y.value, 'wo', mfc='none', ms=25, mew=2)
axs[1].plot(stars_eng.x.value, stars_eng.y.value, 'wo', mfc='none', ms=25, mew=2)

d1.set_limits_minmax(0, 0.3 * image.max())
d2.set_limits_minmax(0, 0.3 * image.max())

axs[0].set_title('CameraFrame')
axs[1].set_title('EngineeringCameraFrame')

plt.show()

@watsonjj
Copy link
Contributor

watsonjj commented Jun 7, 2019

I have had a play around with this PR, along with star field data (Draco constellation) that was recently taken with CHEC. This was done to figure-out/confirm the transformation needed for 2M telescopes. I achieved the following:

image

The notebook for this can be found at:
https://github.com/watsonjj/CHECLabPySB/blob/master/CHECLabPySB/d190607_pointing/engineering.ipynb

CHEC's engineering camera frame is defined by looking at the camera from the secondary mirror. The top of the camera is pointed to the top of the telescope in park position. To convert from this frame to the ctapipe CameraFrame, one must rotate 90 degrees anti-clockwise.

I achieved the frame transformations with the following code:

from ctapipe.coordinates.camera_frame import CameraFrame, CartesianRepresentation, AffineTransform, frame_transform_graph

CAMERA_TO_ENGINEERING2M_MATRIX = np.array([
    [0, 1, 0],
    [-1, 0, 0],
    [0, 0, 1]
])
ENGINEERING2M_FRAME_TRANSFORM_OFFSET = CartesianRepresentation(0, 0, 0, unit=u.m)

ENGINEERING2M_TO_CAMERA_MATRIX = np.array([
    [0, -1, 0],
    [1, 0, 0],
    [0, 0, 1]
])


class Engineering2MCameraFrame(CameraFrame):
    pass

@frame_transform_graph.transform(AffineTransform, CameraFrame, Engineering2MCameraFrame)
def camera_to_engineering2M(from_coord, to_frame):
    return CAMERA_TO_ENGINEERING2M_MATRIX, ENGINEERING2M_FRAME_TRANSFORM_OFFSET

@frame_transform_graph.transform(AffineTransform, Engineering2MCameraFrame, CameraFrame)
def engineering2M_to_camera(from_coord, to_frame):
    return ENGINEERING2M_TO_CAMERA_MATRIX, ENGINEERING2M_FRAME_TRANSFORM_OFFSET

As you can see, two different matrices are required for the transformation back and forth in this case.

My solution was to define a new camera frame Engineering2MCameraFrame...however it would be nicer to avoid having two different EngineeringCameraFrames for the two types of telescopes... I couldn't immediately see how to achieve that. Do you have any ideas?

@maxnoe
Copy link
Member Author

maxnoe commented Jun 8, 2019

That should be easy. Add a frame attribute n_mirrors and return the one or the other matrix dependent on it.

@maxnoe
Copy link
Member Author

maxnoe commented Jun 10, 2019

I added the dual mirror case, @watsonjj can you check if it now works with this PR out of the box using EngineeringFrame(n_mirrors=2)?

@maxnoe maxnoe force-pushed the add_engineering_frame branch from c38dca0 to e98b2ef Compare June 10, 2019 15:57
@watsonjj
Copy link
Contributor

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants