diff --git a/package/lib/fits.py b/package/lib/fits.py
index b75a266..549b0b5 100755
--- a/package/lib/fits.py
+++ b/package/lib/fits.py
@@ -144,8 +144,13 @@ def save_Stokes(
         Only returned if return_hdul is True.
     """
     # Create new WCS object given the modified images
+    wcs = WCS(header_stokes).deepcopy()
     new_wcs = WCS(header_stokes).celestial.deepcopy()
-    header = remove_stokes_axis_from_header(header_stokes).copy() if header_stokes["NAXIS"] > 2 else header_stokes.copy()
+    header = wcs.to_header().copy()
+    header["NAXIS"] = header_stokes["NAXIS"]
+    for i in range(wcs.wcs.naxis):
+        header["NAXIS%d" % (i + 1)] = header_stokes["NAXIS%d" % (i + 1)]
+    header = remove_stokes_axis_from_header(header).copy() if header_stokes["NAXIS"] > 2 else header
 
     if data_mask.shape != (1, 1):
         vertex = clean_ROI(data_mask)
@@ -155,8 +160,8 @@ def save_Stokes(
     for key, val in list(new_wcs.to_header().items()) + [("NAXIS", 2), ("NAXIS1", new_wcs.array_shape[1]), ("NAXIS2", new_wcs.array_shape[0])]:
         header[key] = val
 
-    header["TELESCOP"] = (header_stokes["TELESCOP"] if "TELESCOP" in list(header_stokes.keys()) else "HST", "telescope used to acquire data")
-    header["INSTRUME"] = (header_stokes["INSTRUME"] if "INSTRUME" in list(header_stokes.keys()) else "FOC", "identifier for instrument used to acuire data")
+    header["TELESCOP"] = (header_stokes["TELESCOP"] if "TELESCOP" in list(header_stokes.keys()) else "HST", "Telescope used to acquire data")
+    header["INSTRUME"] = (header_stokes["INSTRUME"] if "INSTRUME" in list(header_stokes.keys()) else "FOC", "identifier for instrument used to acquire data")
     header["PHOTPLAM"] = (header_stokes["PHOTPLAM"], "Pivot Wavelength")
     header["PHOTBW"] = (header_stokes["PHOTBW"], "RMS Bandwidth of the Filter and Detector")
     header["PHOTFLAM"] = (header_stokes["PHOTFLAM"], "Inverse Sensitivity in DN/sec/cm**2/Angst")
diff --git a/package/lib/utils.py b/package/lib/utils.py
index 452a66c..efba45d 100755
--- a/package/lib/utils.py
+++ b/package/lib/utils.py
@@ -303,7 +303,8 @@ def remove_stokes_axis_from_header(header):
     wcs = WCS(header).deepcopy()
     new_wcs = WCS(header).celestial.deepcopy()
     new_head = header.copy()
-    del new_head["NAXIS%d" % (new_wcs.wcs.naxis + 1)]
+    if "NAXIS%d" % (new_wcs.wcs.naxis + 1) in new_head.keys():
+        del new_head["NAXIS%d" % (new_wcs.wcs.naxis + 1)]
     new_head["NAXIS"] = new_wcs.wcs.naxis
     for i, k in enumerate(new_wcs.array_shape[::-1]):
         new_head["NAXIS%d" % (i + 1)] = k