diff --git a/.gitmodules b/.gitmodules index 4e8e7b8..e697650 100644 --- a/.gitmodules +++ b/.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 diff --git a/docs/index.rst b/docs/index.rst index 6075f2e..e4455f5 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -75,6 +75,15 @@ Subarrays :maxdepth: 1 subarrays + +ROP Pipeline +============= + +.. toctree:: + :maxdepth: 1 + + rop_pipeline + Reference/API ============= diff --git a/docs/rop_pipeline.rst b/docs/rop_pipeline.rst new file mode 100644 index 0000000..32af422 --- /dev/null +++ b/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 `_. +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 (,). + 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 (,). + 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 diff --git a/iris_pipeline/datamodels/schemas/tmt_ramp.schema.yaml b/iris_pipeline/datamodels/schemas/tmt_ramp.schema.yaml index ee24c50..6240481 100644 --- a/iris_pipeline/datamodels/schemas/tmt_ramp.schema.yaml +++ b/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 diff --git a/iris_pipeline/drsrop_clib/__init__.py b/iris_pipeline/drsrop_clib/__init__.py new file mode 100644 index 0000000..4b300a6 --- /dev/null +++ b/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"] diff --git a/iris_pipeline/drsrop_clib/drsrop_clib.pyx b/iris_pipeline/drsrop_clib/drsrop_clib.pyx new file mode 100644 index 0000000..45b6d7f --- /dev/null +++ b/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 = 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 = 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 = 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) diff --git a/iris_pipeline/iris_readout b/iris_pipeline/iris_readout new file mode 160000 index 0000000..b6d0b70 --- /dev/null +++ b/iris_pipeline/iris_readout @@ -0,0 +1 @@ +Subproject commit b6d0b70a603e5e7cd5546f2a0cadae79fd5031f8 diff --git a/iris_pipeline/pipeline/__init__.py b/iris_pipeline/pipeline/__init__.py index a08bee5..06d30f4 100644 --- a/iris_pipeline/pipeline/__init__.py +++ b/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"] diff --git a/iris_pipeline/pipeline/rop.py b/iris_pipeline/pipeline/rop.py new file mode 100644 index 0000000..8b033ab --- /dev/null +++ b/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' diff --git a/iris_pipeline/readout/generate_ramp.py b/iris_pipeline/readout/generate_ramp.py new file mode 100644 index 0000000..b0943e4 --- /dev/null +++ b/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') diff --git a/iris_pipeline/readout/nonlincorr_step.py b/iris_pipeline/readout/nonlincorr_step.py new file mode 100644 index 0000000..b224e3c --- /dev/null +++ b/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 diff --git a/iris_pipeline/readout/readoutsamp_step.py b/iris_pipeline/readout/readoutsamp_step.py new file mode 100755 index 0000000..333bb4c --- /dev/null +++ b/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 diff --git a/iris_pipeline/readout/tests/run.py b/iris_pipeline/readout/tests/run.py new file mode 100644 index 0000000..8a277dc --- /dev/null +++ b/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") diff --git a/iris_pipeline/readout/tests/sampling.cfg b/iris_pipeline/readout/tests/sampling.cfg new file mode 100644 index 0000000..61f0922 --- /dev/null +++ b/iris_pipeline/readout/tests/sampling.cfg @@ -0,0 +1,9 @@ +name = "rop" +class = "iris_pipeline.pipeline.ROPPipeline" +save_results = True + + [steps] + [[nonlincorr]] + [[readoutsamp]] + mode='mcds' + diff --git a/iris_pipeline/tests/data/drsrop.cfg b/iris_pipeline/tests/data/drsrop.cfg new file mode 100644 index 0000000..61f0922 --- /dev/null +++ b/iris_pipeline/tests/data/drsrop.cfg @@ -0,0 +1,9 @@ +name = "rop" +class = "iris_pipeline.pipeline.ROPPipeline" +save_results = True + + [steps] + [[nonlincorr]] + [[readoutsamp]] + mode='mcds' + diff --git a/iris_pipeline/tests/test_rop.py b/iris_pipeline/tests/test_rop.py new file mode 100644 index 0000000..a2cf743 --- /dev/null +++ b/iris_pipeline/tests/test_rop.py @@ -0,0 +1,9 @@ +import os +import iris_pipeline +iris_pipeline.monkeypatch_jwst_datamodels() + +def test_rop1(): + iris_pipeline.pipeline.ROPPipeline.call('data/sample_ramp_new.fits', config_file="data/drsrop.cfg") + return 1 + +