wcs graphs are working!!!

Author: seanhess

src/NSO/Image/Asdf.hs

@@ -47,8 +47,6 @@ import Text.Casing (quietSnake)
 -- TODO: fix gwcs
 -- TODO: labeled meta.axes for profiles
 
-
-
 asdfDocument :: Id Inversion -> Dataset -> [Dataset] -> PixelsPerBin -> UTCTime -> L1Asdf -> Frames L2FitsMeta -> Document
 asdfDocument inversionId dscanon dsets bin now l1asdf metas =
   let frames = Frames $ NE.sort metas.frames
@@ -155,7 +153,7 @@ data QuantitiesSection = QuantitiesSection
   , gwcs :: QuantityGWCS
   }
 instance ToAsdf QuantitiesSection where
-  schema _ = "tag:sunpy.org:ndcube/ndcube/ndcollection-1.0.0"
+  schema _ = "tag:sunpy.org:ndcube/ndcollection-1.0.0"
   toValue section =
     Object
       [ ("meta", toNode meta)
@@ -285,7 +283,7 @@ newtype ArmsProfileAxes f = ArmsProfileAxes (Arms (ArmProfile AlignedAxes))
 
 
 instance ToAsdf ProfilesSection where
-  schema _ = "tag:sunpy.org:ndcube/ndcube/ndcollection-1.0.0"
+  schema _ = "tag:sunpy.org:ndcube/ndcollection-1.0.0"
   toValue section =
     Object
       [ ("meta", toNode meta)

src/NSO/Image/GWCS.hs

@@ -8,7 +8,8 @@ import Data.Massiv.Array qualified as M
 import Debug.Trace
 import NSO.Image.Fits.Profile
 import NSO.Image.Fits.Quantity
-import NSO.Image.GWCS.L1GWCS (HPLat, HPLon, L1GWCS (..), L1WCSTransform (..), Time, Zero, l1FixInputs, l1ScaleAxes, l1WCSTransform)
+import NSO.Image.GWCS.L1GWCS (HPLat, HPLon, L1GWCS (..), L1WCSTransform (..), Time, Zero)
+import NSO.Image.GWCS.L1GWCS qualified as L1
 import NSO.Image.Headers (Observation (..), Obsgeo (..))
 import NSO.Image.Headers.Types (Degrees (..), Key (..), Meters (..), PixelsPerBin (..))
 import NSO.Image.Headers.WCS (PC (..), PCXY (..), WCSAxisKeywords (..), WCSCommon (..), WCSHeader (..), Wav, X, Y, toWCSAxis)
@@ -56,16 +57,17 @@ transformQuantity
 transformQuantity bin l1trans axes =
   fullTransform
  where
-  -- dropUnusedZeros fullTransform
-
   fullTransform :: Transform (Pix Depth, Pix X, Pix Y) (Linear Depth, HPLon, HPLat, Time, Zero Wav, Zero Stokes)
   fullTransform =
     let mid = middleFrame axes
      in transformOpticalDepth (toWCSAxis mid.depth.keys) <&> spaceTimeTransform
 
   spaceTimeTransform :: Transform (Pix X, Pix Y) (HPLon, HPLat, Time, Zero Wav, Zero Stokes)
-  spaceTimeTransform = l1FixInputs |> l1ScaleAxes bin |> l1WCSTransform l1trans
+  spaceTimeTransform = L1.fixInputs |> L1.scaleZeroAxes bin |> L1.varyingTransform l1trans
+
 
+-- inverseTransform :: Transform (Linear Depth, HPLon, HPLat, Time, Zero Wav, Zero Stokes) (Pix Depth, Pix X, Pix Y)
+-- inverseTransform = _
 
 -- dropUnusedZeros :: Transform inp (a, b, c, d, Zero z1, Zero z2) -> Transform inp (w, x, y, z)
 -- dropUnusedZeros (Transform t) = Transform t
@@ -209,7 +211,7 @@ quantityGWCS bin l1gwcs primaries quants =
         { name = "optical_depth"
         , axes =
             NE.fromList
-              [ FrameAxis 0 "optical_depth" "optical_depth" Pixel
+              [ FrameAxis 0 "optical depth" "phys.absorption.opticalDepth" Pixel
               ]
         }
 
@@ -270,8 +272,8 @@ celestialFrame n helioFrame =
     , referenceFrame = helioFrame
     , axes =
         NE.fromList
-          [ FrameAxis n "helioprojective longitude" (AxisType "pos.helioprojective.lon") Unit.Arcseconds
-          , FrameAxis (n + 1) "helioprojective latitude" (AxisType "pos.helioprojective.lat") Unit.Arcseconds
+          [ FrameAxis n "pos.helioprojective.lon" (AxisType "pos.helioprojective.lon") Unit.Arcseconds
+          , FrameAxis (n + 1) "pos.helioprojective.lat" (AxisType "pos.helioprojective.lat") Unit.Arcseconds
           ]
     }
 

src/NSO/Image/GWCS/L1GWCS.hs

@@ -5,6 +5,7 @@ import NSO.Image.Headers.WCS (Wav, X, Y)
 import NSO.Prelude hiding (identity)
 import NSO.Types.Common (Stokes)
 import Telescope.Asdf
+import Telescope.Asdf.Core (Quantity (..), Unit (..))
 import Telescope.Asdf.GWCS as GWCS
 import Telescope.Data.Parser
 
@@ -64,8 +65,8 @@ instance ToAsdf L1WCSTransform where
   toValue (L1WCSTransform n) = n.value
 
 
-l1WCSTransform :: L1WCSTransform -> Transform (Scale X, Zero Wav, Pix Y, Zero Stokes) (HPLon, HPLat, Time, Zero Wav, Zero Stokes)
-l1WCSTransform t = originalL1Transform |> reorderOutput
+varyingTransform :: L1WCSTransform -> Transform (Scale X, Zero Wav, Pix Y, Zero Stokes) (HPLon, HPLat, Time, Zero Wav, Zero Stokes)
+varyingTransform t = originalL1Transform |> reorderOutput
  where
   originalL1Transform :: Transform (Scale X, Zero Wav, Pix Y, Zero Stokes) (HPLon, Wav, HPLat, Time, Stokes)
   originalL1Transform = transform t
@@ -74,16 +75,47 @@ l1WCSTransform t = originalL1Transform |> reorderOutput
   reorderOutput = transform $ Mapping [0, 2, 3, 1, 4]
 
 
--- "fix" inputs to unused input axes in the l1 transform: wavlength and stokes
-l1FixInputs :: Transform (Pix X, Pix Y) (Pix X, Pix Wav, Pix Y, Pix Stokes)
-l1FixInputs = transform $ Mapping [0, 0, 1, 0]
-
-
-l1ScaleAxes :: PixelsPerBin -> Transform (Pix X, Pix Wav, Pix Y, Pix Stokes) (Scale X, Zero Wav, Pix Y, Zero Stokes)
-l1ScaleAxes (PixelsPerBin p) = binX <&> zero <&> identity @(Pix Y) <&> zero
+scaleZeroAxes :: PixelsPerBin -> Transform (Pix X, Pix Wav, Pix Y, Pix Stokes) (Scale X, Zero Wav, Pix Y, Zero Stokes)
+scaleZeroAxes (PixelsPerBin p) = binX <&> zero <&> identity @(Pix Y) <&> zero
  where
   binX :: Transform (Pix X) (Scale X)
   binX = scale (fromIntegral p)
 
   zero :: (ToAxes a) => Transform (Pix a) (Zero a)
-  zero = transform $ Scale 0
+  zero = transform $ Const1D $ Quantity Pixel (Integer 0)
+
+
+fixInputs :: Transform (Pix X, Pix Y) (Pix X, Pix Wav, Pix Y, Pix Stokes)
+fixInputs = transform $ Mapping [0, 0, 1, 0]
+
+
+data Const1D = Const1D Quantity
+instance ToAsdf Const1D where
+  schema _ = "!transform/constant-1.2.0"
+  toValue (Const1D q) =
+    Object
+      [ ("dimensions", toNode $ Integer 1)
+      , ("value", toNode q)
+      ]
+instance FromAsdf Const1D where
+  parseValue val = do
+    o <- parseValue @Object val
+    Const1D <$> o .: "value"
+
+-- fixWavStokes0 :: (ToAxes out) => Transform (Pix X, Pix Wav, Pix Y, Pix Stokes) out -> Transform (Pix X, Pix Y) out
+-- fixWavStokes0 (Transform t) = transform $ FixInputs t [(1, 0), (3, 0)]
+--
+--
+-- -- This ONLY works with integers for some reason (doesn't like quantities)
+-- -- combine with Const1D or Scale or something to convert it to a quantity?
+-- data FixInputs = FixInputs Transformation [(Int, Int)]
+-- instance ToAsdf FixInputs where
+--   schema _ = "!transform/fix_inputs-1.2.0"
+--   toValue (FixInputs t fs) =
+--     Object [("forward", toNode [toNode t, toNode forward])]
+--    where
+--     forward =
+--       Object
+--         [ ("keys", toNode $ fmap fst fs)
+--         , ("values", toNode $ fmap snd fs)
+--         ]