ROP pipeline prototype

.gitmodules

@@ -2,3 +2,6 @@
 	url = https://github.com/astropy/astropy-helpers.git
 	path = astropy_helpers
 	branch = refs/heads/v3.1
+[submodule "iris_pipeline/iris_readout"]
+	path = iris_pipeline/iris_readout
+	url = https://github.com/arunsurya77/iris_readout

docs/index.rst

@@ -75,6 +75,15 @@ Subarrays
    :maxdepth: 1
 
    subarrays
+
+ROP Pipeline
+=============
+
+.. toctree::
+   :maxdepth: 1
+
+   rop_pipeline
+
 
 Reference/API
 =============

docs/rop_pipeline.rst

@@ -0,0 +1,97 @@
+DRS-ROP Pipeline
+==========================
+
+The DRS_ROP Pipeline works on detector level readouts. The current steps implemented in the pipeline are
+- Non-linearity Correction
+- Detector Readout Sampling
+
+
+Non-linearity Correction
+------------------------
+Non-linearity correction step corrects for the non-linear response of the detector to incoming flux. This step is executed before sampling algorithms.
+
+
+Detector Readout Sampling 
+------------------------
+The H4RG detecors are readout in non-destructive reads and sampling algorithms are used to estimate the accumulated electrons in the detector for an integration time. The sampling algorithms currently implemented in the pipeline are
+- Correlated Double Sampling
+- Multi Correlated Double Sampling
+- Up-the-Ramp Sampling
+
+Requirements
+------------
+You need following packages to run ROP pipeline
+
+libcfitsio
+cython
+
+Install the iris_pipeline using 
+
+.. code-block:: ini
+
+	git clone --recurse-submodules https://github.com/arunsurya77/iris_pipeline
+
+
+
+Running the Examples
+---------------------
+There is a example run in the iris_pipeline/readout/tests directory. The sample ramp is given in the sample_ramp_new.fits. This file can be extracted from `Figshare <https://figshare.com/articles/sample_ramp_new_fits/12462491>`_.
+There is a minimal pytest script that is in 'iris_pipeline/tests'. This also requires the 'sample_ramp_new.fits' from the link above to be downloaded into iris_pipeline/tests/data.
+
+The sampling.cfg gives the configurations for the pipeline
+
+``sampling.cfg``:
+
+.. code-block:: ini
+
+name = "rop"
+class = "iris_pipeline.pipeline.ROPPipeline"
+save_results = True
+
+    [steps]
+      [[nonlincorr]]
+      [[readoutsamp]]
+       mode='mcds'
+
+
+The sampling mode is set by the ``mode`` keyword which can be ``mcds`` or ``utr``. MCDS algorithm also requires the group number, the number of reads to be co-added. This is currently hardcoded in this version.
+
+
+Execute the pipeline from the command line
+------------------------------------------
+
+We can use ``tmtrun`` from a terminal to execute the pipeline:
+
+::
+
+   tmtrun sampling.cfg sample_ramp.fits
+
+or run using python script
+::
+
+   python run.py
+
+
+here is the output log:
+
+.. code:: bash
+
+		2020-06-09 19:46:04,763 - stpipe.rop - INFO - ROPPipeline instance created.
+		2020-06-09 19:46:04,764 - stpipe.rop.nonlincorr - INFO - NonlincorrStep instance created.
+		2020-06-09 19:46:04,765 - stpipe.rop.readoutsamp - INFO - ReadoutsampStep instance created.
+		2020-06-09 19:46:04,801 - stpipe.rop - INFO - Step rop running with args ('sample_ramp_new.fits',).
+		2020-06-09 19:46:05,405 - stpipe.rop - INFO - Prefetching reference files for dataset: 'sample_ramp_new.fits' reftypes = ['nonlincoeff']
+		2020-06-09 19:46:07,018 - stpipe.rop - INFO - Prefetch for NONLINCOEFF reference file is '/home/arun/crds_cache/references/tmt/iris/tmt_iris_nonlin_coeff.fits'.
+		2020-06-09 19:46:07,019 - stpipe.rop - INFO - Starting ROP Pipeline ...
+		2020-06-09 19:46:07,144 - stpipe.rop.nonlincorr - INFO - Step nonlincorr running with args (<TMTRampModel(1, 4, 10, 10) from sample_ramp_new.fits>,).
+		2020-06-09 19:46:07,283 - stpipe.rop.nonlincorr - INFO - Step nonlincorr done
+		2020-06-09 19:46:07,321 - stpipe.rop.readoutsamp - INFO - Step readoutsamp running with args (<TMTRampModel(1, 4, 10, 10) from sample_ramp_new.fits>,).
+		2020-06-09 19:46:07,343 - stpipe.rop.readoutsamp - INFO - MCDS Sampling Selected
+		(10, 10)
+		2020-06-09 19:46:07,410 - stpipe.rop.readoutsamp - INFO - Step readoutsamp done
+		2020-06-09 19:46:07,528 - stpipe.rop - INFO - Saved model in sample_ramp_new_rop.fits
+		2020-06-09 19:46:07,529 - stpipe.rop - INFO - Step rop done
+
+
+
+This creates a sample_ramp_new_rop.fits file in the working directory that is the processed 

iris_pipeline/datamodels/schemas/tmt_ramp.schema.yaml

@@ -25,14 +25,4 @@ allOf:
       fits_hdu: ERR
       default: 0.0
       datatype: float32
-    zeroframe:
-      title: Zeroframe array
-      fits_hdu: ZEROFRAME
-      default: 0.0
-      ndim: 3
-      datatype: float32
-    group:
-      $ref: http://jwst.stsci.edu/schemas/group.schema.yaml
-    int_times:
-      $ref: http://jwst.stsci.edu/schemas/int_times.schema.yaml
 $schema: http://stsci.edu/schemas/fits-schema/fits-schema

iris_pipeline/drsrop_clib/__init__.py

@@ -0,0 +1,2 @@
+from .drsrop_clib import uptheramp_c,mcds_c,nonlin_c
+__all__ = ["uptheramp_c,mcds_c,nonlin_c"]

iris_pipeline/drsrop_clib/drsrop_clib.pyx

@@ -0,0 +1,18 @@
+cdef extern from "../iris_readout/src/drsrop.c":
+    float* uptheramp (int* arr, int* time, int a, int b, int c)
+    float* mcds (int* arr, int* time, int a, int b, int c, int num_coadd)
+    float* nonlin_corr (int* arr, int* time, int a, int b, int c, int x0, int y0, float* c0, float* c1, float* c2, float* c3, float* c4)
+import numpy as np
+cimport numpy as np
+
+def uptheramp_c(int[:,:,:] arr_f, int[:] arr_time):
+    cdef float[:,:] results = <float[:arr_f.shape[1],:arr_f.shape[2]]> uptheramp(&arr_f[0,0,0], &arr_time[0], arr_f.shape[0], arr_f.shape[1], arr_f.shape[2])
+    return np.asarray(results)
+
+def mcds_c(int[:,:,:] arr_f, int[:] arr_time, int num_coadd):
+    cdef float[:,:] results = <float[:arr_f.shape[1],:arr_f.shape[2]]> mcds(&arr_f[0,0,0], &arr_time[0], arr_f.shape[0], arr_f.shape[1], arr_f.shape[2], num_coadd)
+    return np.asarray(results)
+
+def nonlin_c(int[:,:,:] arr_f, int[:] arr_time, float[:,:] arr_c0, float[:,:] arr_c1, float[:,:] arr_c2, float[:,:] arr_c3, float[:,:] arr_c4):
+    cdef float[:,:,:] results = <float[:arr_f.shape[0],:arr_f.shape[1],:arr_f.shape[2]]> nonlin_corr(&arr_f[0,0,0], &arr_time[0], arr_f.shape[0], arr_f.shape[1], arr_f.shape[2],0,0, &arr_c0[0,0],&arr_c1[0,0],&arr_c2[0,0],&arr_c3[0,0],&arr_c4[0,0] )
+    return np.asarray(results)

iris_pipeline/pipeline/__init__.py

@@ -1,4 +1,4 @@
 from .image2 import Image2Pipeline
 from .preprocess_flatfield import PreprocessFlatfield
-
-__all__ = ["Image2Pipeline", "PreprocessFlatfield"]
+from .rop import ROPPipeline
+__all__ = ["Image2Pipeline", "PreprocessFlatfield","ROPPipeline"]

iris_pipeline/pipeline/rop.py

@@ -0,0 +1,51 @@
+#!/usr/bin/env python
+import logging
+from collections import defaultdict
+import os.path as op
+
+from jwst import datamodels
+from jwst.associations.load_as_asn import LoadAsLevel2Asn
+from jwst.stpipe import Pipeline
+from iris_pipeline import datamodels
+
+# step imports
+from ..readout import readoutsamp_step
+from ..readout import nonlincorr_step
+__all__ = ['rop']
+
+# Define logging
+log = logging.getLogger()
+log.setLevel(logging.DEBUG)
+
+
+class ROPPipeline(Pipeline):
+    """
+    Detector1Pipeline for IRIS
+    """
+
+    spec = """
+        save_calibrated_ramp = boolean(default=False)
+    """
+
+    # Define aliases to steps
+    step_defs = {
+
+         "nonlincorr": nonlincorr_step.NonlincorrStep,
+                 "readoutsamp": readoutsamp_step.ReadoutsampStep,
+                 }
+
+    # start the actual processing
+    def process(self, input):
+        log.info('Starting ROP Pipeline ...')
+        # open the input as a RampModel
+        input = datamodels.TMTRampModel(input)
+        input = self.nonlincorr(input)
+        input = self.readoutsamp(input)
+        return input
+
+    def setup_output(self, input):
+        # Determine the proper file name suffix to use later
+        if input.meta.cal_step.ramp_fit == 'COMPLETE':
+            self.suffix = 'rate'
+        else:
+            self.suffix = 'ramp'

iris_pipeline/readout/generate_ramp.py

@@ -0,0 +1,58 @@
+import pytest
+import numpy as np
+
+from jwst.ramp_fitting.ramp_fit import ramp_fit
+from jwst.datamodels import dqflags
+from jwst.datamodels import MIRIRampModel
+from jwst.datamodels import RampModel
+from jwst.datamodels import GainModel, ReadnoiseModel
+
+
+def setup_inputs(ngroups=4, readnoise=10, nints=1,
+                 nrows=4096, ncols=4096, nframes=1, grouptime=1.0,gain=1, deltatime=1):
+        arr=np.array([np.zeros([4096,4096]),np.ones([4096,4096])*4000,np.ones([4096,4096])*8000,np.ones([4096,4096])*12000],dtype=np.float64)
+        arr=np.array([arr])        
+
+
+        times = np.array(list(range(ngroups)),dtype=np.float64) * deltatime
+        gain = np.ones(shape=(nrows, ncols), dtype=np.float64) * gain
+        err = np.ones(shape=(nints, ngroups, nrows, ncols), dtype=np.float64)
+        data = np.zeros(shape=(nints, ngroups, nrows, ncols), dtype=np.float64)
+        data=arr
+        pixdq = np.zeros(shape=(nrows, ncols), dtype= np.float64)
+        read_noise = np.full((nrows, ncols), readnoise, dtype=np.float64)
+        gdq = np.zeros(shape=(nints, ngroups, nrows, ncols), dtype=np.int32)
+        model1 = RampModel(data=data, err=err, pixeldq=pixdq, groupdq=gdq, times=times)
+        model1.meta.instrument.name='MIRI'
+        model1.meta.instrument.detector='MIRIMAGE'
+        model1.meta.instrument.filter='F480M'
+        model1.meta.observation.date='2015-10-13'
+        model1.meta.exposure.type='MIR_IMAGE'
+        model1.meta.exposure.group_time = deltatime
+        model1.meta.subarray.name='FULL'
+        model1.meta.subarray.xstart=1
+        model1.meta.subarray.ystart = 1
+        model1.meta.subarray.xsize = 10
+        model1.meta.subarray.ysize = 10
+        model1.meta.exposure.frame_time =deltatime
+        model1.meta.exposure.ngroups = ngroups
+        model1.meta.exposure.group_time = deltatime
+        model1.meta.exposure.nframes = 1
+        model1.meta.exposure.groupgap = 0
+        model1.times=times
+        gain = GainModel(data=gain)
+        gain.meta.instrument.name='MIRI'
+        gain.meta.subarray.xstart = 1
+        gain.meta.subarray.ystart = 1
+        gain.meta.subarray.xsize = 10
+        gain.meta.subarray.ysize = 10
+        rnModel = ReadnoiseModel(data=read_noise)
+        rnModel.meta.instrument.name='MIRI'
+        rnModel.meta.subarray.xstart = 1
+        rnModel.meta.subarray.ystart = 1
+        rnModel.meta.subarray.xsize = 10
+        rnModel.meta.subarray.ysize = 10
+        return model1, gdq, rnModel, pixdq, err, gain
+
+model1, gdq, rnModel, pixdq, err, gain = setup_inputs()
+model1.to_fits('test_ramp.fits')

iris_pipeline/readout/nonlincorr_step.py

@@ -0,0 +1,54 @@
+#! /usr/bin/env python
+
+from jwst.stpipe import Step
+from jwst import datamodels
+from ..drsrop_clib import uptheramp_c,mcds_c,nonlin_c
+import numpy as np
+import iris_pipeline
+from astropy.io import fits
+
+__all__ = ["NonlincorrStep"]
+
+
+class NonlincorrStep(Step):
+    """
+    ReadoutsampStep:  Sampling
+    """
+
+    spec = """
+        sigma = float(default=3.0)  # Clipping threshold
+        maxiters = integer(default=None)  # Number of clipping iterations
+    """
+    reference_file_types = ["nonlincoeff"]
+    def process(self, input):
+
+        """
+        Step for Nonlinearity Correction
+        """
+
+        # Load the input data model
+        with datamodels.open(input) as input_model:
+            nonlin_coeff_file=self.get_reference_file(input_model, "nonlincoeff")
+            nonlin_coeff_data=fits.getdata(nonlin_coeff_file).astype(np.float32)                      
+            c0= nonlin_coeff_data[0]
+            c1= nonlin_coeff_data[1]
+            c2= nonlin_coeff_data[2]
+            c3= nonlin_coeff_data[3]
+            c4= nonlin_coeff_data[4]
+            # Get the reference file names
+            input_data = input_model.data[:,:,:,:].astype(np.int32)
+            ramp_list=[]
+            for it in range(len(input_data)):
+                ramp_data=input_data[it]
+                num_reads=len(ramp_data)
+                time_arr=np.arange(0,num_reads,1).astype(np.int32)
+                result=nonlin_c(ramp_data,time_arr,c0,c1,c2,c3,c4)
+                result=result.astype(np.int32)
+                ramp_list.append(result)            
+            #result = uptheramp_c(result,time_arr)
+            output_data=np.array(ramp_list).astype(np.int32)
+            result = iris_pipeline.datamodels.TMTRampModel(data=output_data)
+            result.update(input_model)
+
+
+        return result

iris_pipeline/readout/readoutsamp_step.py

@@ -0,0 +1,56 @@
+#! /usr/bin/env python
+
+from jwst.stpipe import Step
+from jwst import datamodels
+import iris_pipeline
+from ..drsrop_clib import uptheramp_c,mcds_c,nonlin_c
+import numpy as np
+
+__all__ = ["ReadoutsampStep"]
+
+
+class ReadoutsampStep(Step):
+    """
+    ReadoutsampStep:  Sampling
+    """
+
+    spec = """
+        sigma = float(default=3.0)  # Clipping threshold
+        maxiters = integer(default=None)  # Number of clipping iterations
+        mode= string(default='utr')
+    """
+
+    def process(self, input):
+
+        """
+        Sampling step.
+        Currently 3 Sampling modes
+         -- CDS
+         -- MCDS
+         -- UTR
+        """
+
+        # Load the input data model
+        with datamodels.open(input) as input_model:
+
+            # Get the reference file names
+            input_data = input_model.data[:,:,:,:].astype(np.int32)
+            im_list=[]
+            for it in range(len(input_data)):
+                ramp_data=input_data[it]
+                num_reads=len(ramp_data)
+                time_arr=np.arange(0,num_reads,1).astype(np.int32)                              
+                if(self.mode=='mcds'):
+                        self.log.info("MCDS Sampling Selected")
+                        result = mcds_c(ramp_data,time_arr,2)
+                else:
+                        self.log.info("UTR Sampling Selected")
+                        result = utr(ramp_data,time_arr,2)
+                im_list.append(result)
+            result=np.sum(im_list,axis=0) 
+            print(result.shape)
+            result = iris_pipeline.datamodels.IRISImageModel(data=result)
+            result.update(input_model)
+
+
+        return result

iris_pipeline/readout/tests/run.py

@@ -0,0 +1,7 @@
+import os
+import iris_pipeline
+import jwst.datamodels
+print(iris_pipeline.__file__)
+iris_pipeline.monkeypatch_jwst_datamodels()
+#from iris_pipeline.pipeline import ROPPipeline
+iris_pipeline.pipeline.ROPPipeline.call('sample_ramp_new.fits', config_file="sampling.cfg")