quantity conversion. profile shape and bad hdu bug

nso-level2.cabal

@@ -218,6 +218,7 @@ test-suite spec
       Test.FrameSpec
       Test.ParseSpec
       Test.QualifySpec
+      Test.QuantitySpec
       Paths_nso_level2
   autogen-modules:
       Paths_nso_level2

src/NSO/Image/Asdf.hs

@@ -94,7 +94,7 @@ data InversionTree = InversionTree
   { fileuris :: Fileuris
   , meta :: InversionTreeMeta
   , quantities :: HDUSection QuantityGWCS (Quantities (DataTree QuantityMeta))
-  , profiles :: HDUSection ProfileGWCS (Profiles ProfileTree)
+  , profiles :: ProfileSection
   }
   deriving (Generic, ToAsdf)
 
@@ -146,6 +146,25 @@ instance (ToAsdf hdus, ToAsdf gwcs) => ToAsdf (HDUSection gwcs hdus) where
       ]
 
 
+data ProfileSection = ProfileSection
+  { axes :: [AxisLabel]
+  , wcs :: ProfileGWCS
+  , hdus :: Profiles ProfileTree
+  }
+
+
+instance ToAsdf ProfileSection where
+  toValue section =
+    mconcat
+      -- merge the fields from both
+      [ Object
+          [ ("axes", toNode section.axes)
+          , ("wcs", toNode section.wcs)
+          ]
+      , toValue section.hdus
+      ]
+
+
 -- Quantities ------------------------------------------------
 
 quantitiesSection :: NonEmpty FrameQuantitiesMeta -> HDUSection QuantityGWCS (Quantities (DataTree QuantityMeta))
@@ -246,26 +265,25 @@ instance (KnownText ref) => KnownText (Ref ref) where
 
 -- Profiles ------------------------------------------------
 
-profilesSection :: NonEmpty FrameProfilesMeta -> HDUSection ProfileGWCS (Profiles ProfileTree)
+profilesSection :: NonEmpty FrameProfilesMeta -> ProfileSection
 profilesSection frames =
-  let shape = (head frames).shape
-      wcs = (head frames).profiles.orig630.wcs :: WCSHeader ProfileAxes
-   in HDUSection
+  let wcs = (head frames).profiles.orig630.wcs :: WCSHeader ProfileAxes
+   in ProfileSection
         { wcs = profileGWCS wcs
         , axes = ["frameY", "slitX", "wavelength", "stokes"]
-        , shape = shape.axes
-        , hdus = profilesTree shape frames
+        , hdus = profilesTree frames
         }
 
 
-profilesTree :: Shape Profile -> NonEmpty FrameProfilesMeta -> Profiles ProfileTree
-profilesTree shape frames =
-  let ps = fmap (.profiles) frames
+profilesTree :: NonEmpty FrameProfilesMeta -> Profiles ProfileTree
+profilesTree frames =
+  let frame = head frames
+      ps = fmap (.profiles) frames
    in Profiles
-        { orig630 = profileTree shape $ fmap (.orig630) ps
-        , orig854 = profileTree shape $ fmap (.orig854) ps
-        , fit630 = profileTree shape $ fmap (.fit630) ps
-        , fit854 = profileTree shape $ fmap (.fit854) ps
+        { orig630 = profileTree frame.shape630 $ fmap (.orig630) ps
+        , orig854 = profileTree frame.shape854 $ fmap (.orig854) ps
+        , fit630 = profileTree frame.shape630 $ fmap (.fit630) ps
+        , fit854 = profileTree frame.shape854 $ fmap (.fit854) ps
         }
 
 

src/NSO/Image/Frame.hs

@@ -63,7 +63,8 @@ data FrameQuantitiesMeta = FrameQuantitiesMeta
 
 
 data FrameProfilesMeta = FrameProfilesMeta
-  { shape :: Shape Profile
+  { shape630 :: Shape Profile
+  , shape854 :: Shape Profile
   , profiles :: Profiles ProfileHeader
   }
   deriving (Generic)
@@ -145,21 +146,26 @@ frameMetaFromL2Fits path slice wpo wpf l1 fits = runParser $ do
   qshape <- parseHeader @(Shape Quantity) qh
   quants <- parseQuantities
 
-  ph <- headerFor @Orig630
-  pshape <- parseHeader @(Shape Profile) ph
+  p630 <- headerFor @Orig630
+  shape630 <- parseHeader @(Shape Profile) p630
+
+  p854 <- headerFor @Orig854
+  shape854 <- parseHeader @(Shape Profile) p854
+
   profs <- parseProfiles
   pure $
     L2FrameMeta
       { path
       , primary
       , quantities = FrameQuantitiesMeta{quantities = quants, shape = qshape}
-      , profiles = FrameProfilesMeta{profiles = profs, shape = pshape}
+      , profiles = FrameProfilesMeta{profiles = profs, shape630, shape854}
       }
  where
   headerFor :: forall a es. (HDUOrder a, Parser :> es) => Eff es Header
   headerFor =
     let HDUIndex index = hduIndex @a
-     in case fits.extensions !? index of
+        extIndex = index - 1
+     in case fits.extensions !? extIndex of
           Nothing -> parseFail $ "Missing HDU at " ++ show index
           Just (Image img) -> pure img.header
           Just _ -> parseFail $ "Expected ImageHDU at " ++ show index
@@ -224,7 +230,8 @@ frameMeta frame path =
   profilesMeta ps =
     FrameProfilesMeta
       { profiles = profileHeaders ps
-      , shape = Shape $ addDummyY $ dataCubeAxes ps.orig630.image
+      , shape630 = Shape $ addDummyY $ dataCubeAxes ps.orig630.image
+      , shape854 = Shape $ addDummyY $ dataCubeAxes ps.orig854.image
       }
 
   addDummyY :: Axes Row -> Axes Row
@@ -283,10 +290,10 @@ instance HDUOrder GasPressure where
 instance HDUOrder Density where
   hduIndex = 11
 instance HDUOrder Orig630 where
-  hduIndex = 11
-instance HDUOrder Orig854 where
   hduIndex = 12
-instance HDUOrder Fit630 where
+instance HDUOrder Orig854 where
   hduIndex = 13
-instance HDUOrder Fit854 where
+instance HDUOrder Fit630 where
   hduIndex = 14
+instance HDUOrder Fit854 where
+  hduIndex = 15

src/NSO/Image/Quantity.hs

@@ -6,6 +6,7 @@ module NSO.Image.Quantity where
 
 import Control.Exception (Exception)
 import Data.ByteString (ByteString)
+import Data.Fixed (mod')
 import Data.Massiv.Array as M (Index, Ix2 (..), IxN (..), Sz (..), map)
 import Data.Text (pack)
 import Data.Time.Format.ISO8601 (iso8601Show)
@@ -234,9 +235,9 @@ quantityHDUs qs = runPureEff $ do
   temperature <- dataHDU @Temperature qs.temperature
   electronPressure <- dataHDU @ElectronPressure $ convertData dyneCmToNm qs.electronPressure
   microTurbulence <- dataHDU @Microturbulence $ convertData cmsToKms qs.microTurbulence
-  magStrength <- dataHDU @MagStrength qs.magStrength
+  magStrength <- dataHDU @MagStrength $ convertData gaussToTesla qs.magStrength
   magInclination <- dataHDU @MagInclination qs.magInclination
-  magAzimuth <- dataHDU @MagAzimuth qs.magAzimuth
+  magAzimuth <- dataHDU @MagAzimuth $ convertData forcePositive360 qs.magAzimuth
   geoHeight <- dataHDU @GeoHeight qs.geoHeight
   gasPressure <- dataHDU @GasPressure $ convertData dyneCmToNm qs.gasPressure
   density <- dataHDU @Density $ convertData gcmToKgm qs.density
@@ -250,12 +251,6 @@ quantityHDUs qs = runPureEff $ do
       , image = DataCube $ M.map f da.array
       }
 
-  cmsToKms = (/ 100000)
-
-  gcmToKgm = (* 1000)
-
-  dyneCmToNm = (/ 10)
-
   dataHDU
     :: forall info es
      . (ToHeader info)
@@ -266,6 +261,26 @@ quantityHDUs qs = runPureEff $ do
     pure $ QuantityHDU $ ImageHDU{header = toHeader q.header, dataArray = addDummyAxis darr}
 
 
+cmsToKms :: Float -> Float
+cmsToKms = (/ 100000)
+
+
+gcmToKgm :: Float -> Float
+gcmToKgm = (* 1000)
+
+
+dyneCmToNm :: Float -> Float
+dyneCmToNm = (/ 10)
+
+
+gaussToTesla :: Float -> Float
+gaussToTesla = (/ 10000)
+
+
+forcePositive360 :: Float -> Float
+forcePositive360 deg = deg `mod'` 360
+
+
 data QuantityAxes alt = QuantityAxes
   { depth :: QuantityAxis alt Depth
   , slitX :: QuantityAxis alt X
@@ -335,7 +350,7 @@ wcsDepth :: (Monad m) => m (QuantityAxis alt n)
 wcsDepth = do
   let crpix = Key 12
       crval = Key 0
-      cdelt = Key 0.1
+      cdelt = Key (-0.1)
       cunit = Key ""
       ctype = Key "LOGTAU"
   let keys = WCSAxisKeywords{..}

test/Test/FrameSpec.hs

@@ -23,6 +23,9 @@ spec = do
         M.toLists (computeAs P f0.array) `shouldBe` ([[[0, 1, 2, 3], [1, 2, 3, 4]]] :: [[[Float]]])
 
 
+
+
+
 sample4D :: DataCube [Q, Depth, FrameY, SlitX]
 sample4D = DataCube $ M.makeArray @D @Ix4 Seq (Sz (1 :> 2 :> 3 :. 4)) sumIndex
  where

test/Test/QuantitySpec.hs

@@ -0,0 +1,16 @@
+module Test.QuantitySpec where
+
+import Data.Massiv.Array as M
+import NSO.Image.Quantity (forcePositive360)
+import NSO.Prelude
+import Skeletest
+
+
+spec :: Spec
+spec =
+  describe "Quantities" $ do
+    describe "convertData" $ do
+      it "should move negative angles to positive" $ do
+        let arr = M.fromLists' @P @Ix1 @Float Seq [-100, -200, 0, 30, 300, 400]
+        let arr2 = M.map forcePositive360 arr
+        M.toLists (computeAs P arr2) `shouldBe` [260, 160, 0, 30, 300, 40]