finalizing slit plot and coordinate retrieval

gris/GrisFitsFile.py

@@ -120,7 +120,7 @@ class GrisFitsFile(FitsFile):
         Y = y_slitcenter - np.sin(angle) * delta_pix * arcsec_per_pixel
 
         UserWarning(
-            "Using FITS file coordinates, they are probably wrong by at least 100 arcsec!"
+            "Using FITS file coordinates, they are probably wrong by at least 100 arcsec!, Also the angle is off!"
         )
 
         return X, Y
@@ -374,41 +374,55 @@ class GrisFitsFile(FitsFile):
     def spatial_coords(self):
         wcs = WCS(self.header)
         pixel_shape = wcs.pixel_shape
-        n_x = pixel_shape[0]
         n_y = pixel_shape[1]
-        n_dim = len(pixel_shape)
-
-        slit_list = np.zeros((n_x, 4))
-        slit_list[:, 0] = np.arange(n_x)
-        slit_list[:, 1:] = 0
-
+        slit_list = np.zeros((n_y, 4))
+        slit_list[:, 1] = np.arange(n_y)
+        slit_list[:, 0] = 0
         converted = wcs.pixel_to_world_values(slit_list)
         degrees = u.degree * converted[:, :2]
         arcsec = degrees.to(u.arcsec)
+
         return arcsec
 
     def plot_slit(self):
         m = self.full_disk_map
 
-        fig = plt.figure()  # noqa F841
-        ax = plt.subplot(projection=m)
-        im = m.plot()  # noqa:F841
-
+        # retrieve different data sources
         X, Y = self.spatial_coords[:, 0], self.spatial_coords[:, 1]
         x_slitpos, y_slitpos = self.slit_coords * u.arcsec
-        slit = SkyCoord(X, Y, frame=m.coordinate_frame)
-        ax.plot_coord(slit, 'g-')
-        header = self.header
-        slitpos = SkyCoord(x_slitpos, y_slitpos, frame=m.coordinate_frame)
-        ax.plot_coord(slitpos, 'ro-')
-        x_ref, y_ref = header["CRVAL1"] * u.arcsec, header["CRVAL2"] * u.arcsec
-        ref_point = SkyCoord(x_ref, y_ref, frame=m.coordinate_frame)
-        ax.plot_coord(ref_point, 'bo-')
-
-        for k, v in header.items():
-            if k.startswith("C"): print(f"{k}:{v}")
-
-        print(self.obs_time)
+        x_ref, y_ref = self.header["CRVAL1"] * u.arcsec, self.header["CRVAL2"] * u.arcsec
+
+        # determine fov limits
+        x_min = min([min(X), min(x_slitpos), x_ref]) - 20 * u.arcsec
+        x_max = max([max(X), max(x_slitpos), x_ref]) + 20 * u.arcsec
+        y_min = min([min(Y), min(y_slitpos), y_ref]) - 20 * u.arcsec
+        y_max = max([max(Y), max(y_slitpos), y_ref]) + 20 * u.arcsec
+
+        # Create submap of FOV and plot
+        submap = m.submap(SkyCoord(x_min, y_min, frame=m.coordinate_frame),
+                          SkyCoord(x_max, y_max, frame=m.coordinate_frame))
+        fig = plt.figure()  # noqa F841
+        ax = plt.subplot(projection=submap)
+        im = submap.plot()  # noqa:F841
+
+        # create SkyCoord objects and overplot
+        slitpos = SkyCoord(x_slitpos, y_slitpos, frame=submap.coordinate_frame)
+        slit_wcs = SkyCoord(X, Y, frame=submap.coordinate_frame)
+        ref_point = SkyCoord(x_ref, y_ref, frame=submap.coordinate_frame)
+        ax.plot_coord(ref_point, 'bo', label="(CRVAL1,CRVAL2)")
+        ax.plot_coord(slitpos, 'r--', label="Slitpos")
+        ax.plot_coord(slit_wcs, 'g-', label="WCS Coordinates")
+
+        # Create Legend and annotations
+        title = "Gris Slit {0}_{1:03d}_{2:03d}_{3:03d}"
+        title = title.format(self.obs_time.strftime('%Y%m%d'), self.runnumber, self.mapnumber, self.step_index)
+        plt.title(title)
+
+        annotation = f"Spectrum taken at \n{self.obs_time.strftime('%Y-%m-%d %H:%M:%S')}"
+        annotation += f"\nHMI Reference taken at \n{submap.date.strftime('%Y-%m-%d %H:%M:%S')}"
+        plt.text(0.05, 0.8, annotation, horizontalalignment='left', transform=ax.transAxes)
+        plt.tight_layout(pad=3.5)
+        _ = plt.legend()
 
         return fig, ax