How to use the ctapipe.coordinates.CameraFrame function in ctapipe

To help you get started, we’ve selected a few ctapipe examples, based on popular ways it is used in public projects.

Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.

github cta-observatory / ctapipe / examples / coordinate_transformations.py View on Github external
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
github cta-observatory / ctapipe / ctapipe / tools / muon_reconstruction.py View on Github external
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
github cta-observatory / ctapipe / ctapipe / reco / HillasReconstructor.py View on Github external
)

            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),
            )
github cta-observatory / ctapipe / examples / coordinate_transformations.py View on Github external
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)
github cta-observatory / ctapipe / ctapipe / reco / HillasReconstructor.py View on Github external
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
github cta-observatory / ctapipe / ctapipe / reco / HillasReconstructor.py View on Github external
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)
github cta-observatory / ctapipe / examples / read_hessio_single_tel.py View on Github external
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)
github cta-observatory / ctapipe / ctapipe / image / muon / intensity_fitter.py View on Github external
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)
github cta-observatory / ctapipe / ctapipe / image / muon / muon_reco_functions.py View on Github external
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")