Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
def cam_to_tel():
# Coordinates in any fram can be given as a numpy array of the xyz positions
# e.g. in this case the position on pixels in the camera
pix_x = np.ones(2048) * u.m
pix_y = np.ones(2048) * u.m
# first define the camera frame
camera_frame = CameraFrame(focal_length=15 * u.m)
# create a coordinate in that frame
camera_coord = SkyCoord(x=pix_x, y=pix_y, frame=camera_frame)
# then use transform to function to convert to a new system making sure
# to give the required values for the conversion (these are not checked
# yet)
telescope_coord = camera_coord.transform_to(TelescopeFrame())
# Print coordinates in the new frame
print("Telescope Coordinate", telescope_coord)
# Transforming back is then easy
camera_coord2 = telescope_coord.transform_to(camera_frame)
# We can easily check the distance between 2 coordinates in the same frame
# In this case they should be the same
def get_pixel_coords(self, tel_id):
"""Get pixel coords in telescope frame for telescope with id `tel_id`"""
# memoize transformation
if tel_id not in self.pixels_in_tel_frame:
telescope = self.source.subarray.tel[tel_id]
cam = telescope.camera.geometry
camera_frame = CameraFrame(
focal_length=telescope.optics.equivalent_focal_length,
rotation=cam.cam_rotation,
)
cam_coords = SkyCoord(x=cam.pix_x, y=cam.pix_y, frame=camera_frame)
tel_coord = cam_coords.transform_to(TelescopeFrame())
self.pixels_in_tel_frame[tel_id] = tel_coord
coords = self.pixels_in_tel_frame[tel_id]
return coords.fov_lon, coords.fov_lat
)
camera_frame = CameraFrame(
focal_length=focal_length, telescope_pointing=pointing
)
cog_coord = SkyCoord(x=moments.x, y=moments.y, frame=camera_frame,)
cog_coord = cog_coord.transform_to(horizon_frame)
p2_coord = SkyCoord(x=p2_x, y=p2_y, frame=camera_frame)
p2_coord = p2_coord.transform_to(horizon_frame)
# re-project from sky to a "fake"-parallel-pointing telescope
# then recalculate the psi angle
if self.divergent_mode:
camera_frame_parallel = CameraFrame(
focal_length=focal_length, telescope_pointing=array_pointing
)
cog_sky_to_parallel = cog_coord.transform_to(camera_frame_parallel)
p2_sky_to_parallel = p2_coord.transform_to(camera_frame_parallel)
angle_psi_corr = np.arctan2(
cog_sky_to_parallel.y - p2_sky_to_parallel.y,
cog_sky_to_parallel.x - p2_sky_to_parallel.x,
)
self.corrected_angle_dict[tel_id] = angle_psi_corr
circle = HillasPlane(
p1=cog_coord,
p2=p2_coord,
telescope_position=subarray.positions[tel_id],
weight=moments.intensity * (moments.length / moments.width),
)
def cam_to_nom():
pix_x = np.ones(2048) * u.m
pix_y = np.ones(2048) * u.m
pointing_direction = SkyCoord(alt=70 * u.deg, az=180 * u.deg, frame=AltAz())
camera_frame = CameraFrame(
focal_length=15 * u.m, telescope_pointing=pointing_direction
)
camera_coord = SkyCoord(pix_x, pix_y, frame=camera_frame)
# In this case we bypass the telescope system
nominal_frame = NominalFrame(origin=AltAz(alt=75 * u.deg, az=180 * u.deg))
nom_coord = camera_coord.transform_to(nominal_frame)
horizon = camera_coord.transform_to(AltAz())
print("Nominal Coordinate", nom_coord)
print("Horizon coordinate", horizon)
for tel_id, moments in hillas_dict.items():
focal_length = subarray.tel[tel_id].optics.equivalent_focal_length
pointing = SkyCoord(
alt=pointing_alt[tel_id],
az=pointing_az[tel_id],
frame='altaz'
)
hf = HorizonFrame(
array_direction=pointing,
pointing_direction=pointing
)
cf = CameraFrame(
focal_length=focal_length,
array_direction=pointing,
pointing_direction=pointing
)
cog_coord = SkyCoord(x=moments.x, y=moments.y, frame=cf)
cog_coord = cog_coord.transform_to(hf)
cog_direction = spherical_to_cartesian(1, cog_coord.alt, cog_coord.az)
cog_direction = np.array(cog_direction).ravel()
weights.append(self.hillas_planes[tel_id].weight)
tels.append(self.hillas_planes[tel_id].pos)
dirs.append(cog_direction)
# minimising the test function
self.hillas_planes = {}
k = next(iter(telescopes_pointings))
horizon_frame = telescopes_pointings[k].frame
for tel_id, moments in hillas_dict.items():
# we just need any point on the main shower axis a bit away from the cog
p2_x = moments.x + 0.1 * u.m * np.cos(moments.psi)
p2_y = moments.y + 0.1 * u.m * np.sin(moments.psi)
focal_length = subarray.tel[tel_id].optics.equivalent_focal_length
pointing = SkyCoord(
alt=telescopes_pointings[tel_id].alt,
az=telescopes_pointings[tel_id].az,
frame=horizon_frame,
)
camera_frame = CameraFrame(
focal_length=focal_length, telescope_pointing=pointing
)
cog_coord = SkyCoord(x=moments.x, y=moments.y, frame=camera_frame,)
cog_coord = cog_coord.transform_to(horizon_frame)
p2_coord = SkyCoord(x=p2_x, y=p2_y, frame=camera_frame)
p2_coord = p2_coord.transform_to(horizon_frame)
# re-project from sky to a "fake"-parallel-pointing telescope
# then recalculate the psi angle
if self.divergent_mode:
camera_frame_parallel = CameraFrame(
focal_length=focal_length, telescope_pointing=array_pointing
)
cog_sky_to_parallel = cog_coord.transform_to(camera_frame_parallel)
disp.set_limits_percent(70)
plt.suptitle("Sample {:03d}".format(ii))
plt.pause(0.01)
if args.write:
plt.savefig('CT{:03d}_EV{:010d}_S{:02d}.png'
.format(args.tel, event.r0.event_id, ii))
else:
# display integrated event:
im = event.r0.tel[args.tel].adc_sums[args.channel]
peds, gains = get_mc_calibration_coeffs(event, args.tel)
im = apply_mc_calibration(im, peds, gains, args.tel)
disp.image = im
if args.hillas:
clean_mask = ctapipe.image.cleaning.tailcuts_clean(geom, im, 1, picture_thresh=10, boundary_thresh=5)
camera_coord = CameraFrame(x=x,y=y,z=np.zeros(x.shape)*u.m)
nom_coord = camera_coord.transform_to(NominalFrame(array_direction=[70*u.deg,0*u.deg],
pointing_direction=[70*u.deg,0*u.deg],
focal_length=tel['TelescopeTable_VersionFeb2016'][tel['TelescopeTable_VersionFeb2016']['TelID']==args.tel]['FL'][0]*u.m))
image = np.asanyarray(im * clean_mask, dtype=np.float64)
nom_x = nom_coord.x
nom_y = nom_coord.y
hillas = reco.hillas_parameters(x,y,im * clean_mask)
hillas_nom = reco.hillas_parameters(nom_x,nom_y,im * clean_mask)
print (hillas)
print (hillas_nom)
hole_radius=0 * u.m,
):
"""Create an efficient negative log_likelihood function that does
not rely on astropy units internally by defining needed values as closures
in this function
"""
# get all the neeed values and transform them into appropriate units
optics = telescope_description.optics
mirror_area = optics.mirror_area.to_value(u.m ** 2)
mirror_radius = np.sqrt(mirror_area / np.pi)
focal_length = optics.equivalent_focal_length
cam = telescope_description.camera.geometry
camera_frame = CameraFrame(focal_length=focal_length, rotation=cam.cam_rotation)
cam_coords = SkyCoord(x=cam.pix_x, y=cam.pix_y, frame=camera_frame)
tel_coords = cam_coords.transform_to(TelescopeFrame())
# Use only a subset of pixels, indicated by mask:
pixel_x = tel_coords.fov_lon.to_value(u.rad)
pixel_y = tel_coords.fov_lat.to_value(u.rad)
if mask is not None:
pixel_x = pixel_x[mask]
pixel_y = pixel_y[mask]
image = image[mask]
pedestal = pedestal[mask]
pixel_diameter = 2 * (
np.sqrt(cam.pix_area[0] / np.pi) / focal_length * u.rad
).to_value(u.rad)
geom, image, picture_thresh=tailcuts[0], boundary_thresh=tailcuts[1]
)
# TODO: correct this hack for values over 90
altval = event.mcheader.run_array_direction[1]
if altval > Angle(90, unit=u.deg):
warnings.warn("Altitude over 90 degrees")
altval = Angle(90, unit=u.deg)
telescope_pointing = SkyCoord(
alt=altval, az=event.mcheader.run_array_direction[0], frame=AltAz()
)
camera_coord = SkyCoord(
x=x,
y=y,
frame=CameraFrame(
focal_length=teldes.optics.equivalent_focal_length,
rotation=geom.pix_rotation,
telescope_pointing=telescope_pointing,
),
)
nom_coord = camera_coord.transform_to(NominalFrame(origin=telescope_pointing))
x = nom_coord.delta_az.to(u.deg)
y = nom_coord.delta_alt.to(u.deg)
if cleaning:
img = image * clean_mask
else:
img = image
muon_ring_fit = MuonRingFitter(fit_method="kundu_chaudhuri")