"Snowball" masking with snowblind

grizli/aws/field_tiles.py

@@ -494,7 +494,10 @@ def make_all_tile_images(root, force=False, ref_tile=(8,8), cleanup=True, zoom_l
             if f in all_filters:
                 filters.append(f)
                 break
-
+
+        if ('cos' in root) & ('f150w2-clear' in all_filters) & (len(filters) < 3):
+            filters = ['f814w','f150w2-clear','f444w-clear']
+
         split_tiles(root, ref_tile=ref_tile, 
                     filters=filters,
                     zoom_levels=zoom_levels,

grizli/aws/visit_processor.py

@@ -273,7 +273,9 @@ def s3_put_exposure(flt_file, product, assoc, remove_old=True, verbose=True):
         print(f'Add {file}_{extension} ({len(rows)}) to exposure_files table')
 
 
-def make_visit_mosaic(assoc, base_path=ROOT_PATH, version='v7.0', pixscale=0.08, vmax=0.5, skip_existing=True, sync=True, clean=False, verbose=True, **kwargs):
+snowblind_kwargs = dict(require_prefix='jw', max_fraction=0.3, new_jump_flag=1024, min_radius=4, growth_factor=1.5, unset_first=True, verbose=True)
+
+def make_visit_mosaic(assoc, base_path=ROOT_PATH, version='v7.0', pixscale=0.08, vmax=0.5, skip_existing=True, sync=True, clean=False, verbose=True, snowblind_kwargs=snowblind_kwargs, **kwargs):
     """
     Make a mosaic of the exposures from a visit with a tangent point selected
     from the sky tile grid
@@ -421,6 +423,7 @@ def make_visit_mosaic(assoc, base_path=ROOT_PATH, version='v7.0', pixscale=0.08,
                  pixfrac=0.8,
                  res=res,
                  weight_type=weight_type,
+                 snowblind_kwargs=snowblind_kwargs,
     )
 
     files = glob.glob(f'{assoc}*_sci.fits*')
@@ -1676,7 +1679,7 @@ def process_visit(assoc, clean=True, sync=True, max_dt=4, combine_same_pa=False,
         add_shifts_log(assoc=assoc, remove_old=True, verbose=True)
         add_wcs_log(assoc=assoc, remove_old=True, verbose=True)
 
-    if do_make_visit_mosaic:
+    if (do_make_visit_mosaic) & (with_db):
         make_visit_mosaic(assoc, sync=sync,
                           base_path=ROOT_PATH,
                           clean=False,
@@ -2022,7 +2025,7 @@ def query_exposures(ra=53.16, dec=-27.79, size=1., pixel_scale=0.1, theta=0,  fi
     return header, wcs, SQL, res
 
 
-def cutout_mosaic(rootname='gds', product='{rootname}-{f}', ra=53.1615666, dec=-27.7910651, size=5*60, theta=0., filters=['F160W'], ir_scale=0.1, ir_wcs=None, res=None, half_optical=True, kernel='point', pixfrac=0.33, make_figure=True, skip_existing=True, clean_flt=True, gzip_output=True, s3output='s3://grizli-v2/HST/Pipeline/Mosaic/', split_uvis=True, extra_query='', extra_wfc3ir_badpix=True, fix_niriss=True, scale_nanojy=10, verbose=True, weight_type='err', rnoise_percentile=99, get_dbmask=True, niriss_ghost_kwargs={}, scale_photom=True, context='jwst_0989.pmap', calc_wcsmap=False, make_exptime_map=False, expmap_sample_factor=4, keep_expmap_small=True, **kwargs):
+def cutout_mosaic(rootname='gds', product='{rootname}-{f}', ra=53.1615666, dec=-27.7910651, size=5*60, theta=0., filters=['F160W'], ir_scale=0.1, ir_wcs=None, res=None, half_optical=True, kernel='point', pixfrac=0.33, make_figure=True, skip_existing=True, clean_flt=True, gzip_output=True, s3output='s3://grizli-v2/HST/Pipeline/Mosaic/', split_uvis=True, extra_query='', extra_wfc3ir_badpix=True, fix_niriss=True, scale_nanojy=10, verbose=True, weight_type='err', rnoise_percentile=99, get_dbmask=True, niriss_ghost_kwargs={}, snowblind_kwargs=None, scale_photom=True, context='jwst_0989.pmap', calc_wcsmap=False, make_exptime_map=False, expmap_sample_factor=4, keep_expmap_small=True, **kwargs):
     """
     Make mosaic from exposures defined in the exposure database
 
@@ -2296,6 +2299,7 @@ def cutout_mosaic(rootname='gds', product='{rootname}-{f}', ra=53.1615666, dec=-
                                      context=context,
                                      get_dbmask=get_dbmask,
                                      niriss_ghost_kwargs=niriss_ghost_kwargs,
+                                     snowblind_kwargs=snowblind_kwargs,
                                      weight_type=weight_type,
                                      calc_wcsmap=calc_wcsmap)
 
@@ -2939,7 +2943,138 @@ def run_one(clean=2, sync=True):
             fp.write(f'{time.ctime()} {assoc}\n')
 
         process_visit(assoc, clean=clean, sync=sync)
+
+
+def snowblind_exposure_mask(assoc, file, verbose=True, clean=True, new_jump_flag=1024, min_radius=4, growth_factor=1.5, unset_first=True, run_again=False, force=False):
+    """
+    Run snowblind on a previously processed exposures
+    """
+    import time
+    import boto3
+    import astropy.io.fits as pyfits
+
+    if 0:
+        db.execute("drop table exposure_snowblind")
+        db.execute("""create table exposure_snowblind as 
+        select assoc, dataset || '_' || extension || '.fits' as file, instrume 
+        from exposure_files 
+        where instrume in ('NIRCAM','NIRISS')
+        group by assoc, dataset, extension, instrume
+        """)
+        db.execute("alter table exposure_snowblind add column status int default 70")
+        db.execute("alter table exposure_snowblind add column maskfrac real default 0.")
+        db.execute("alter table exposure_snowblind add column npix int default 0")
+        db.execute("alter table exposure_snowblind add column utime double precision default 0.")
+
+    status = db.SQL(f"""select * from exposure_snowblind
+    where assoc = '{assoc}' and file = '{file}'
+    """)
+    if len(status) == 0:
+        msg = f"File {assoc} {file} not found"
+        utils.log_comment(utils.LOGFILE, msg, verbose=verbose)
+        return False
+
+    if (status['status'].max() != 0) & (~force):
+        msg = "{assoc}  {file}  status={status} locked".format(**status[0])
+        utils.log_comment(utils.LOGFILE, msg, verbose=verbose)
+        return False
+
+    # set lock status
+    db.execute(f"""update exposure_snowblind set status = 1 
+    where assoc = '{assoc}' and file = '{file}'
+    """)
+
+    # Run the snowblind mask    
+    s3 = boto3.resource('s3')
+    s3_client = boto3.client('s3')
+    bkt = s3.Bucket('grizli-v2')
+
+    local_file = os.path.basename(file)
+    s3_prefix = f'HST/Pipeline/{assoc}/Prep/{file}'
 
+    if os.path.exists(local_file):
+        clean = False
+    else:
+        msg = f'Fetch {s3_prefix}'
+        utils.log_comment(utils.LOGFILE, msg, verbose=verbose)
+        bkt.download_file(s3_prefix, local_file,
+                          ExtraArgs={"RequestPayer": "requester"})
+
+    do_run = True
+    with pyfits.open(local_file) as im:
+        dq = im['DQ'].data & new_jump_flag
+        maskfrac = (dq > 0).sum() / dq.size
+        if 'SNOWBLND' in im['SCI'].header:
+            do_run = run_again
+
+    if do_run:        
+        dq, maskfrac = utils.jwst_snowblind_mask(local_file,
+                                       new_jump_flag=new_jump_flag,
+                                       min_radius=min_radius,
+                                       growth_factor=growth_factor,
+                                       unset_first=unset_first)
+
+    with pyfits.open(local_file, mode='update') as im:
+        if unset_first:
+            im['DQ'].data -= im['DQ'].data & new_jump_flag
+
+        im['DQ'].data |= dq.astype(im['DQ'].data.dtype)
+        im['SCI'].header['SNOWMASK'] = (True, 'Snowball mask applied')
+        im['SCI'].header['SNOWBLND'] = (True, 'Mask with snowblind')
+        im['SCI'].header['SNOWBALF'] = (maskfrac,
+                             'Fraction of masked pixels in snowball mask')
+
+        im.flush()
+        npix = ((im['DQ'].data & new_jump_flag) > 0).sum()
+
+    # Upload back to s3
+    msg = f'Send {file} > s3://grizli-v2/{s3_prefix}'
+    utils.log_comment(utils.LOGFILE, msg, verbose=verbose)
+    bkt.upload_file(local_file, s3_prefix, ExtraArgs={'ACL': 'public-read'},
+                    Callback=None, Config=None)
 
+    # Update status
+    db.execute(f"""update exposure_snowblind
+    set status = 2, maskfrac = {maskfrac}, npix = {npix}, utime = {time.time()}
+    where assoc = '{assoc}' and file = '{file}'
+    """)
+
+    if clean:
+        os.remove(local_file)
+
+    return True
+
+
+def snowblind_batch():
+    import os
+    os.chdir('/GrizliImaging')
+
+    db.execute("""update exposure_snowblind set status = 70 
+    where (file not like 'jw01895%%' and assoc not like 'j03%%m27%%') and status = 0
+    """)
+
+    db.execute("""update exposure_snowblind set status = 0 
+    where (file like 'jw01895%%' and (assoc like 'j03%%m27%%' or assoc like 'j12%%p62%%')) and status = 70
+    """)
+
+    files = db.SQL("""select * from exposure_snowblind
+    where (file like '%%_rate.fits') and status = 0
+    order by random()
+    """)
+
+    files = db.SQL("""select * from exposure_snowblind
+    where (file like '%%_rate.fits') and status = 0 and assoc like '%%f444w%%'
+    order by random()
+    """)
+
+    print(f'Run {len(files)} files ....\n')
+    for row in files:
+        snowblind_exposure_mask(row['assoc'], row['file'], verbose=True, clean=True, 
+                                new_jump_flag=1024,
+                                min_radius=4,
+                                growth_factor=1.5, unset_first=True, run_again=False)
+        # break
+
+
 if __name__ == '__main__':
     run_all()

grizli/data/README.hotpix.txt

@@ -0,0 +1,16 @@
+# SW additional hot pixels
+
+jwst_nircam_newhot* from Chris Willott
+
+Oringal filenames:
+
+    jwst_nircam_newhot_a1_0062_extra20231129.fits     jwst_nircam_newhot_NRCA1_extra20231129.fits     
+    jwst_nircam_newhot_a2_0068_extra20231129.fits     jwst_nircam_newhot_NRCA2_extra20231129.fits     
+    jwst_nircam_newhot_a3_0067_extra20231129.fits     jwst_nircam_newhot_NRCA3_extra20231129.fits     
+    jwst_nircam_newhot_a4_0066_extra20231129.fits     jwst_nircam_newhot_NRCA4_extra20231129.fits     
+    jwst_nircam_newhot_b1_0060_extra20231129.fits     jwst_nircam_newhot_NRCB1_extra20231129.fits     
+    jwst_nircam_newhot_b2_0064_extra20231129.fits     jwst_nircam_newhot_NRCB2_extra20231129.fits     
+    jwst_nircam_newhot_b3_0069_extra20231129.fits     jwst_nircam_newhot_NRCB3_extra20231129.fits     
+    jwst_nircam_newhot_b4_0061_extra20231129.fits     jwst_nircam_newhot_NRCB4_extra20231129.fits     
+    jwst_nircam_newhot_blong_0065_extra20231129.fits  jwst_nircam_newhot_NRCBLONG_extra20231129.fits 
+    jwst_nircam_newhot_along_0063_extra20231129.fits  jwst_nircam_newhot_NRCALONG_extra20231129.fits  

grizli/data/auto_script_defaults.yml

@@ -141,6 +141,11 @@ visit_prep_args: &prep_args
         mask_bit: 1024
         instruments: [NIRCAM, NIRISS]
         max_fraction: 0.3
+        snowblind_kwargs:
+            new_jump_flag: 1024
+            min_radius: 4
+            growth_factor: 1.5
+            unset_first: True
     angle_background_kwargs:
         threshold: 1.8
         detection_background: True

grizli/data/nircam_edge_bpix.py

@@ -0,0 +1,26 @@
+"""
+Mask pixels at the edge of NIRCam LW
+"""
+import astropy.io.fits as pyfits
+import glob
+import numpy as np
+import os
+
+os.system('gunzip nrc_badpix_231206_*fits.gz')
+
+nircam_edge = 8
+
+files = glob.glob('nrc_badpix_231206_*fits')
+files.sort()
+for file in files:
+    with pyfits.open(file, mode='update') as im:
+        dq = np.zeros_like(im[0].data)
+        dq[:nircam_edge,:] |= 1024
+        dq[-nircam_edge:,:] |= 1024
+        dq[:,:nircam_edge] |= 1024
+        dq[:,-nircam_edge:] |= 1024
+        im[0].data |= (dq > 0)
+
+        im.flush()
+
+os.system('gzip nrc_badpix_231206_*fits')

grizli/jwst_utils.py

@@ -1276,7 +1276,11 @@ def initialize_jwst_image(filename, verbose=True, max_dq_bit=14, orig_keys=ORIG_
     elif img[0].header['OINSTRUM'] == 'NIRCAM':
         _det = img[0].header['DETECTOR']
 
-        bpfiles = [os.path.join(os.path.dirname(__file__), 
+        bpfiles = [os.path.join(os.path.dirname(__file__),
+                   f'data/nrc_badpix_231206_{_det}.fits.gz')]
+        bpfiles += [os.path.join(os.path.dirname(__file__),
+                   f'data/nrc_badpix_20230710_{_det}.fits.gz')]
+        bpfiles += [os.path.join(os.path.dirname(__file__), 
                    f'data/nrc_badpix_230120_{_det}.fits.gz')]
         bpfiles += [os.path.join(os.path.dirname(__file__), 
                    f'data/nrc_lowpix_0916_{_det}.fits.gz')]

grizli/prep.py

@@ -3377,8 +3377,10 @@ def oneoverf_column_correction(visit, thresholds=[10,1.5], dilate_iter=[10,2], i
 
             _im.flush()
 
+SNOWBLIND_KWARGS = dict(new_jump_flag=1024, min_radius=4,
+                        growth_factor=1.5, unset_first=True)
 
-def mask_snowballs(visit, snowball_erode=3, snowball_dilate=18, mask_bit=1024, instruments=['NIRCAM','NIRISS'], max_fraction=0.3, unset4=False, **kwargs):
+def mask_snowballs(visit, snowball_erode=3, snowball_dilate=18, mask_bit=1024, instruments=['NIRCAM','NIRISS'], max_fraction=0.3, unset4=False, snowblind_kwargs=SNOWBLIND_KWARGS, **kwargs):
     """
     Mask JWST IR snowballs
 
@@ -3405,7 +3407,10 @@ def mask_snowballs(visit, snowball_erode=3, snowball_dilate=18, mask_bit=1024, i
 
     unset4 : bool
         Unset DQ=4 bit for flagged CRs
-
+
+    snowblind_kwargs : dict
+        First try to use `grizli.utils.jwst_snowblind_mask` to flag snowballs
+
     Returns
     -------
     Updates `DQ` extension of files in ``visit['files']`` and sets some
@@ -3427,7 +3432,26 @@ def mask_snowballs(visit, snowball_erode=3, snowball_dilate=18, mask_bit=1024, i
 
             if _instrume not in instruments:
                 continue
-
+
+            if snowblind_kwargs is not None:
+                sdq, sfrac = utils.jwst_snowblind_mask(_file, **snowblind_kwargs)
+                if sdq is not None:
+                    # Close and reopen
+                    _im.close()
+                    with pyfits.open(_file, mode='update') as _xim:
+                        _xim['DQ'].data |= sdq.astype(_xim['DQ'].data.dtype)
+                        _xim['SCI'].header['SNOWMASK'] = (True, 'Snowball mask applied')
+                        _xim['SCI'].header['SNOWBLND'] = (True, 'Mask with snowblind')
+                        _xim['SCI'].header['SNOWBALF'] = (sfrac,
+                                          'Fraction of masked pixels in snowball mask')
+
+                        if unset4:
+                            _xim['DQ'].data -= (_xim['DQ'].data & 4)
+
+                        _xim.flush()
+
+                    continue
+
             crs = (_im['DQ'].data & 4)
 
             _erode = nd.binary_erosion(crs > 0, iterations=snowball_erode)