<jupyter_start><jupyter_text>Simulating Grism Images with ```pyLINEAR```This notebook demonstrates how to use ```pyLINEAR``` to simulate _simple_ grism images, here *simple* refers to the assumption that there is a single SED for a source.  In future notebooks, I will show how this assumption can be addressed on small scales.  The command-line interface (CLI) ```pylinear.py``` can be used to streamline a majority of the operations contained in this notebook.  However, this highlights how the steps can be deconstructed and reapplied to meet your individual needs or "pipeline".  Please check the help menu for the CLI: ```prompt> pylinear -h```; I expect many of the variables will be familiar after you have gone through this notebook. Notebook history:1. initial development: <NAME> Jul. 23, 20202. updates for CLI: <NAME> Aug. 10, 2020<jupyter_code>import pylinear  # import pylinear
from astropy.io import fits   # for reading/writing fits files
import numpy as np            # for random numbers
import glob                   # to make input lists of files
import shutil                 # really just to copy one file
import os                     # operations with file names
segfile = 'seg.fits'           # segmentation we already built
obslst  = 'obs.lst'            # name of the obslst
wcslst  = 'wcs.lst'            # the name of the WCS list we will be creating
beam    = '+1'                 # the beam to simulate
sedlst  = 'sed.lst'            # name of the SED list
seddir  = 'SEDs'               # temporary dir for the SEDs to create
maglim  = 99.0                 # a high magnitude limit (take everything)<jupyter_output><empty_output><jupyter_text>Step 1. Load the sources.The sources are loaded into a ```SourceCollection``` object, which has 2 mandatory inputs:1. ```segfile``` the full path to the segmentation map2. ```obslst```  the full path to a list of direct imagesThe ```obslst``` is an ascii file that takes the form:```    img1.fits     hst_wfc3_f105w    img2.fits     hst_wfc3_f125w    ```The first column is the full path to the direct image and the second indicates the filter associated with this file.  The filters are stored in ```pylinear/sedphot/filters/``` and have an entry in the file ```pylinear/sedphot/filters.py```.  To add more filters to the collection, simply put the ascii file (with two columns of wavelength and transmission) in the appropriate directory and add more entries to the ```filters.py``` file.  **Note, the zeropoint there should be in AB mag units.**Additionally, the ```SourceCollection``` has an optional argument and default: ```detindex=0```, which specifies which row of the ```obs.lst``` represents the image to use as extraction weights.  In principle, this can be any wavelength image you'd like, but in practice this should be as close as possible to the central wavelength of the grism.  So, the default setting of ```detindex=0``` means that the first row of the ```obs.lst``` file will be used as the extraction weight image.<jupyter_code># step 1a. Make the obslst file
with open(obslst,'w') as fp:
    print('img.fits   hst_wfc3_f105w',file=fp)

# here, I am assuming we will work with HST/WFC3 in the F105W filter
# step 1b. Load the sources
sources = pylinear.source.SourceCollection(segfile,obslst,detindex=0,maglim=maglim)<jupyter_output>3: 100%|██████████| 3/3 [00:00<00:00, 131.60it/s]<jupyter_text>Step 2.  Load grism imagesThe grism images will be loaded into a ```GrismCollection``` object, which has one mandatory argument:1. ```img.lst``` This variable changes definition depending on the value of the optional argument ```observed```.  For ```observed=False```, we will be simulating grism images and the file that ```img.lst``` refers to is an ascii file contains the WCS properties of the image(s) to simulate (we will create this file, so the format will be discussed below).  If ```observed=True```, then we are using real observations and so the ```img.lst``` points to a list of files to be loaded as observations (discussed in later notebooks).<jupyter_code># step 2a. create a WCS list

with open(wcslst,'w') as fp:
    # obligatory header info:
    print('# TELESCOPE = HST',file=fp)   # specify the telescope
    print('# INSTRUMENT = WFC3',file=fp) # specify the instrument
    print('# DETECTOR = IR',file=fp)     # specify the detector 
    print('# BLOCKING = ',file=fp)       # specify the blocking filter (only for JWST)
    
    # now put in rows for each image to simulate
    print('simp1   53.000  -27.000   0.  G102',file=fp)
    print('simp2   53.001  -27.001   0.  G102',file=fp)
    print('simp3   53.000  -27.000  45.  G102',file=fp)
    print('simp4   53.001  -27.001  45.  G102',file=fp)
    print('simp5   53.000  -27.000  90.  G102',file=fp)
    print('simp6   53.001  -27.001  90.  G102',file=fp)
    print('simp7   53.000  -27.000 130.  G102',file=fp)
    print('simp8   53.001  -27.001 130.  G102',file=fp)
    print('#file9   53.000  -27.000 175.  G102',file=fp)
    
# the columns are:
# 1. DATASET (recalling DATASET is the ipppssoot id of a given HST image),
# 2. RA (RA of where to center the camera)
# 3. Dec (Dec of where to center the camera)
# 4. ORIENTAT (angle E of N to rotate the camera
# 5. grism (the grism element to simulate)
#
# pyLINEAR will automatically append "_flt.fits" to files as relevant. Rows can be 
# commented out with '#', '%', or ';', such as done in the last entry.  The DATASET
# will be interpreted as a string, so you can make a file root a numeric charater.  
# Finally, dithers are implemented by tweaking the (RA,Dec) fields, as is the 
# difference between rows 1 and 2, for example.
# step 2b. now can load the grism images (in this case will be simulated images)
grisms = pylinear.grism.GrismCollection(wcslst,observed=False)<jupyter_output><empty_output><jupyter_text>Step 3. Make pixel-dispersion tables (PDTs)All table creation is handled by the ```Tabulate``` object, which takes one mandatory argument (the type of table to create).  There are primarily two table types one typically makes: ```pdt``` and ```omt``` for pixel-dispersion and object-masking, respectively.  For simulations, we do not need to mask anything, so we will just be using ```pdt```.  This step is by far the most computationaly intensive step, and so is multithreaded using the ```multiprocessing``` module.  Therefore the optional argument ```ncpu=None``` is used to control how many local computing resources should be used.  If this is set to ```None``` or ```<=0```, then ```pyLINEAR``` will use 1 less than the total number on your system (the thought here is to leave some CPU resources to do minor work).<jupyter_code># step 3a. create a tabulate object:
tabulate = pylinear.modules.Tabulate('pdt',ncpu=0)   # note the ncpu = 0<jupyter_output><empty_output><jupyter_text>Now to run the ```tabulate``` module, one must pass the grism and source collections as well as specify which beam to tabulate.  Here we are only processing a single beam (as indicated above), but one can embed this in a small loop to process multiple beams.  The ```run``` method performs the calculations and returns a list of filenames that it created.<jupyter_code># step 3b. run the tabulate module
tabnames = tabulate.run(grisms,sources,beam)<jupyter_output>Making tables:   0%|          | 0/8 [00:00<?, ?it/s]WARNING: FITSFixedWarning: The WCS transformation has more axes (2) than the image it is associated with (0) [astropy.wcs.wcs]
WARNING: FITSFixedWarning: The WCS transformation has more axes (2) than the image it is associated with (0) [astropy.wcs.wcs]
WARNING: FITSFixedWarning: The WCS transformation has more axes (2) than the image it is associated with (0) [astropy.wcs.wcs]
WARNING: FITSFixedWarning: The WCS transformation has more axes (2) than the image it is associated with (0) [astropy.wcs.wcs]
WARNING: FITSFixedWarning: The WCS transformation has more axes (2) than the image it is associated with (0) [astropy.wcs.wcs]
WARNING: FITSFixedWarning: The WCS transformation has more axes (2) than the image it is associated with (0) [astropy.wcs.wcs]
Making tables: 100%|██████████| 8/8 [00:00<00:00, 12.40it/s]<jupyter_text>Step 4. Simulate the grism exposuresGrism image simulation is handled by the ```Simulate``` module, which has one mandatory argument ```sed.lst```.  This is an ascii file that specifies the spectrum to disperse for each object in the ```segfile```, and has the form COLUMN 1: an integer index for the SEGMENTATION ID and COLUMN 2: the full path to an SED to use.<jupyter_code># 4a. create the Simulate object
simulate = pylinear.modules.Simulate(sedlst,gzip=False)

# pylinear by default wants to gzip the files, but let's disable that for the moment<jupyter_output><empty_output><jupyter_text>Step 4b. Create some SEDs.In principle, this is not really pyLINEAR *per se*, as one may already have the SEDs to use. If you have SEDs, then just make an ```sed.lst``` file and skip this step.  If not, keep reading to cook up some temporary SEDs.<jupyter_code># step 4b. Create some nonsense SEDs

# let's create a directory so the files don't clutter up the workspace
if not os.path.isdir(seddir):
    os.mkdir(seddir)

# create a single wavelength grid (not necessary, just for ease of use later)
wav = np.arange(7000,13001,1)

# create the SED file names
sed1,sed2,sed3 = seddir+'/1.sed',seddir+'/2.sed',seddir+'/3.sed'

# open all the files
with open(sed1,'w') as f1, open(sed2,'w') as f2, open(sed3,'w') as f3,open(sedlst,'w') as fp:
    # do each wavelength 
    for w in wav:
        
        # assume power-law SED for object 1
        print('{} {}'.format(w,(w/10000.)**(-2.)),file=f1)
        
        # assume step-function SED for object 2
        print('{} {}'.format(w,float(w > 10000.)),file=f2)
        
        # assume sine-function SED for object 3 (not meant to be physical)
        print('{} {}'.format(w,1.+np.sin((w-7000.)/6000.*2.*np.pi)),file=f3)
       
    #put the name of each SED in the SEDLIST
    print('1 {}'.format(sed1),file=fp)
    print('2 {}'.format(sed2),file=fp)
    print('3 {}'.format(sed3),file=fp)
# 4c. run the simulation
fltnames = simulate.run(grisms,sources,beam)<jupyter_output>Simulating FLTs:   0%|          | 0/8 [00:00<?, ?it/s]<jupyter_text>step 4d.  Just an FYI: Will be used LaterThe SED for each object is loaded and renormalized such that the bandpass-integrated flux matches that in the direct image.  Therefore, the SED that you specified is ***NOT*** the SED in the grism image (well, the same up to a multiplicative constant).  Therefore ```pyLINEAR``` writes the ***ACTUALLY USED*** SED to a directory: ```simulated_SEDs```.  **At present, this is not configurable.**  Step 5. Inspect the resultsAt this point, ```pyLINEAR``` should have produced a host of noiseless WFC3/IR grism exposures for G102.  I would inspect them and recommend ```ds9```.  ```pyLINEAR``` attempts to emulate the natural ```_flt.fits``` file structure, and as such, these files should be multi-extension fits (MEF) files.  Therefore, I would load these accordingly so you can see the 3 extensions (```SCI```, ```ERR```, and ```DQ```) and ```EXTVER``` accordingly, and verify that they make sense to what you expect.  There are two things to note:1. at this point, you have not specified anything about the noise model, and so the ```ERR``` extension is just set to unity.2. ```pyLINEAR``` does not know anything about the data-quality arrays.  This might be something to implement in the future, but at present the DQAs are initialized to all zeros.  Step 6. (optional) Add photometric noiseWe have taken the philosophical choice that ```pyLINEAR``` does not attempt to understand the noise characteristics of your particular grism.  However, for full understanding of the algorithm, it is important to have some noise in the system.  For the impetus, try removing this step and see how ```pyLINEAR``` fares.  As a robust treatment of the noise would also require implementing the correct sky-background spectrum, this step is also omitted.  Therefore we will emulate a set of observations whose background signal has been removed, but sky noise remains.For the simple tests, let's add a Gaussian random noise field (with zero mean) to our data.  We will start with a very low sky sigma to approximate very bright sources.  If you ratchet this sigma down, then you will need to do one of the following:1. increase the brightness of your sources;2. live with lower S/N of the extractions; or 3. add more simulated grism images.But as our goal is to demonstrate the code, I think the unrealistically low sigma is ok for now.We are making a few simplications here about WFC3/IR, and detectors with multiple chips (eg. WFC3/UVIS, ACS/WFC) will need extra layers of generalization.<jupyter_code>sig = 0.001    # noise RMS in e-/s (note the above commentary)


# lets apply the noise to all the files:
for oldf in glob.glob('simp*_flt.fits'):
    print('Adding noise to {}'.format(oldf))
    
    # let's save the file in case we want to compare
    savefile = oldf+'_noiseless'
    shutil.copyfile(oldf,savefile)

    # open the fits file
    with fits.open(oldf) as hdul:
        sci = hdul[('SCI',1)].data    # the science image
        size = sci.shape              # dimensionality of the image
        
        # update the science extension with random noise
        hdul[('SCI',1)].data = sci + np.random.normal(loc=0.,scale=sig,size=size)
        
        # update the uncertainty extension with the sigma
        hdul[('ERR',1)].data = np.full_like(sci,sig)
        
        # now write to a new file name
        hdul.writeto(oldf,overwrite=True)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09892677228590202, 'fraction_numerical': 0.026064101628093743, 'mean_word_length': 3.872998932764141, 'pattern_counts': {'":': 0, '<': 37, '<?xml version=': 0, '>': 34, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '65837', 'license': 'MIT', 'n_tokens_mistral': 4354, 'n_tokens_neox': 4025, 'n_words': 1830, 'path': 'notebooks/pylinear_simulate_grism.ipynb', 'repo_name': 'Russell-Ryan/pyLINEAR', 'size': '13709'}