l1 may 28 dload + imgdirx check

download_from_db.py

@@ -19,7 +19,7 @@ text = file.read()
 
 text_split = get_file_list(text)
 
-target_directory = '/scratch/slam/sinjan/solo_attic_fits/fits_files/'
+target_directory = '/data/slam/home/sinjan/fits_files/'
 
 os.chdir(target_directory)
 
@@ -33,5 +33,6 @@ for file in text_split:
   else:
     subprocess.call(["wget", "--user", f"{username}","--password", f"{password}", f"https://www2.mps.mpg.de/services/proton/phi/imgdb/{file}"])
     subprocess.call(["gunzip", f"{file}"])
+    subprocess.call(["chmod", "a+r", f"{file}"])
 
 print("Download and unpacking complete")

file_names.txt

-solo_L1_phi-hrt-ilam_20200529T205144_V202106111602C_0065260200.fits.gz
+l1/2020-04-17/solo_L1_phi-hrt-ilam_20200417T112613_V202106111548C_0024150001.fits.gz
+l1/2020-04-17/solo_L1_phi-hrt-ilam_20200417T114945_V202106111548C_0024150005.fits.gz
+l1/2020-04-17/solo_L1_phi-hrt-ilam_20200417T125349_V202106111549C_0024150006.fits.gz
+l1/2020-04-17/solo_L1_phi-hrt-ilam_20200417T125349_V202106111549C_0064151006.fits.gz
+l1/2020-04-17/solo_L1_phi-hrt-ilam_20200417T132752_V202106111549C_0024150010.fits.gz
+l1/2020-04-17/solo_L1_phi-hrt-ilam_20200417T151154_V202106111549C_0024150002.fits.gz
+l1/2020-04-17/solo_L1_phi-hrt-ilam_20200417T151455_V202106111549C_0024150003.fits.gz
+l1/2020-04-17/solo_L1_phi-hrt-ilam_20200417T163923_V202106111549C_0024150004.fits.gz
+l1/2020-04-17/solo_L1_phi-hrt-ilam_20200417T165722_V202106111549C_0024150012.fits.gz
+l1/2020-04-17/solo_L1_phi-hrt-ilam_20200417T174538_V202106111549C_0024150020.fits.gz
+l1/2020-04-17/solo_L1_phi-hrt-ilam_20200417T182538_V202106111550C_0024150021.fits.gz
+l1/2020-04-17/solo_L1_phi-hrt-ilam_20200417T190538_V202106111550C_0024150022.fits.gz
+l1/2020-04-17/solo_L1_phi-hrt-ilam_20200417T194538_V202106111550C_0024150023.fits.gz
+l1/2020-04-17/solo_L1_phi-hrt-ilam_20200417T202538_V202106111551C_0024150024.fits.gz
+l1/2020-04-17/solo_L1_phi-hrt-ilam_20200417T210538_V202106111551C_0024150025.fits.gz
+l1/2020-04-17/solo_L1_phi-hrt-ilam_20200417T214539_V202106111551C_0024150026.fits.gz
+l1/2020-04-17/solo_L1_phi-hrt-ilam_20200417T222538_V202106111552C_0024150027.fits.gz
+l1/2020-04-17/solo_L1_phi-hrt-ilam_20200417T230538_V202106111552C_0024150028.fits.gz
+l1/2020-04-17/solo_L1_phi-hrt-ilam_20200417T234538_V202106111552C_0024150029.fits.gz

hrt_pipe.py

@@ -9,15 +9,16 @@ from operator import itemgetter
 from utils import *
 
 
-def get_data(path, scaling = True):
+def get_data(path, scaling = True, bit_convert_scale = True):
     """load science data from path"""
     try:
         data, header = load_fits(path)
 
-        if scaling:
-      
-            data /=  256. #conversion from 24.8bit to 32bit
+        if bit_convert_scale: #conversion from 24.8bit to 32bit
+                data /=  256.
 
+        if scaling:
+            
             accu = header['ACCACCUM']*header['ACCROWIT']*header['ACCCOLIT'] #getting the number of accu from header
 
             data /= accu
@@ -72,7 +73,7 @@ def demod_hrt(data,pmp_temp):
     
 
 
-def phihrt_pipe(data_f, dark_f = '', flat_f = '', scale_data = True, norm_f = True, clean_f = False, 
+def phihrt_pipe(data_f, dark_f = '', flat_f = '', scale_data = True, bit_flat = True, norm_f = True, clean_f = False, 
                 sigma = 59, flat_states = 24, prefilter_f = None,flat_c = True, dark_c = True, fs_c = True,
                 demod = True, norm_stokes = True, out_dir = './',  out_demod_file = False,  
                 correct_ghost = False, ItoQUV = False, ctalk_params = None, rte = False, out_rte_filename = ''):
@@ -193,11 +194,11 @@ def phihrt_pipe(data_f, dark_f = '', flat_f = '', scale_data = True, norm_f = Tr
 
             if scale_data:
 
-                if hdr_arr[scan]['BITPIX'] == 16:
+                #if hdr_arr[scan]['BITPIX'] == 16:
 
-                    print(f"This scan: {data_f[scan]} has a bits per pixel of: 16 \nPerforming the extra scaling")
+                #    print(f"This scan: {data_f[scan]} has a bits per pixel of: 16 \nPerforming the extra scaling")
 
-                    data_arr[scan] *= 81920/127 #conversion factor if 16 bits
+                data_arr[scan] *= 81920/127 #conversion factor if 16 bits
 
             if 'IMGDIRX' in hdr_arr[scan] and hdr_arr[scan]['IMGDIRX'] == 'YES':
                 print("This scan has been flipped in the Y axis to conform to orientation standards.")
@@ -251,6 +252,21 @@ def phihrt_pipe(data_f, dark_f = '', flat_f = '', scale_data = True, norm_f = Tr
 
         print(f"Data shape is {data.shape}")
 
+        #--------
+        # test if the scans have same IMGDIRX keyword
+        #--------
+        if all('IMGDIRX' in hdr for hdr in hdr_arr):
+            first_imgdirx = hdr_arr[0]['IMGDIRX']
+            result = all(hdr['IMGDIRX'] == first_imgdirx for hdr in hdr_arr)
+            if (result):
+                print(f"All the scans have the same IMGDIRX keyword: {first_imgdirx}")
+                imgdirx_flipped = str(first_imgdirx)
+            else:
+                print("The scans have different IMGDIRX keywords! Please fix \n Ending Process")
+                exit()
+        else:
+            print("Not all the scans contain the 'IMGDIRX' keyword! Assuming not flipped")
+
     elif isinstance(data_f, str):
         #case when data f is just one file
         data, header = get_data(data_f, scaling = scale_data)
@@ -284,6 +300,13 @@ def phihrt_pipe(data_f, dark_f = '', flat_f = '', scale_data = True, norm_f = Tr
         pmp_temp = str(header['HPMPTSP1'])
         print(f"Data PMP Set Point Temperature is {pmp_temp}")
 
+        if 'IMGDIRX' in header:
+            imgdirx_flipped = str(header['IMGDIRX'])
+
+        else:
+            print("Scan header does not contain the 'IMGDIRX' keyword! Assuming not flipped")
+
+
     else:
         printc("ERROR, data_f argument is neither a string nor list containing strings: {} \n Ending Process",data_f,color=bcolors.FAIL)
         exit()
@@ -324,51 +347,51 @@ def phihrt_pipe(data_f, dark_f = '', flat_f = '', scale_data = True, norm_f = Tr
 
         start_time = time.time()
     
-        try:
-            flat, header_flat = get_data(flat_f)
+        #try:
+        flat, header_flat = get_data(flat_f, bit_convert_scale=bit_flat)
 
-            print(f"Flat field shape is {flat.shape}")
-            
-            if header_flat['BITPIX'] == 16:
-    
-                print("Number of bits per pixel is: 16")
+        print(f"Flat field shape is {flat.shape}")
+        
+        if header_flat['BITPIX'] == 16:
 
-                flat *= 614400/128
+            print("Number of bits per pixel is: 16")
 
-            flat = np.moveaxis(flat, 0,-1) #so that it is [y,x,24]
-            flat = flat.reshape(2048,2048,6,4) #separate 24 images, into 6 wavelengths, with each 4 pol states
-            flat = np.moveaxis(flat, 2,-1)
-            
-            print(flat.shape)
+            flat *= 614400/128
 
-            _, _, _, cpos_f = fits_get_sampling(flat_f,verbose = True)
+        flat = np.moveaxis(flat, 0,-1) #so that it is [y,x,24]
+        flat = flat.reshape(2048,2048,6,4) #separate 24 images, into 6 wavelengths, with each 4 pol states
+        flat = np.moveaxis(flat, 2,-1)
+        
+        print(flat.shape)
+
+        _, _, _, cpos_f = fits_get_sampling(flat_f,verbose = True)
 
-            print(f"The continuum position of the flat field is at {cpos_f} index position")
+        print(f"The continuum position of the flat field is at {cpos_f} index position")
+        
+        if flat_f[-62:] == 'solo_L0_phi-hrt-flat_0667134081_V202103221851C_0162201100.fits':
             
-            if flat_f[-62:] == 'solo_L0_phi-hrt-flat_0667134081_V202103221851C_0162201100.fits':
-                
-                print("This flat needs to be rolled")
-                
-                flat = np.roll(flat, 1, axis = -1)
+            print("This flat needs to be rolled")
+            
+            flat = np.roll(flat, 1, axis = -1)
 
-                cpos_f = cpos_arr[0]
+            cpos_f = cpos_arr[0]
 
-            #if continuum wavelength of the flat is not the same as the data, attempt to roll
-            if cpos_f != cpos_arr[0]:
-                print("The flat field continuum position is not the same as the data, please check your input data. \n Ending Process")
+        #if continuum wavelength of the flat is not the same as the data, attempt to roll
+        if cpos_f != cpos_arr[0]:
+            print("The flat field continuum position is not the same as the data, please check your input data. \n Ending Process")
 
-                exit()
+            exit()
 
-            flat_pmp_temp = str(header_flat['HPMPTSP1'])
+        flat_pmp_temp = str(header_flat['HPMPTSP1'])
 
-            print(f"Flat PMP Temperature Set Point: {flat_pmp_temp}")
-                
-            printc('--------------------------------------------------------------',bcolors.OKGREEN)
-            printc(f"------------ Load flats time: {np.round(time.time() - start_time,3)} seconds",bcolors.OKGREEN)
-            printc('--------------------------------------------------------------',bcolors.OKGREEN)
+        print(f"Flat PMP Temperature Set Point: {flat_pmp_temp}")
+            
+        printc('--------------------------------------------------------------',bcolors.OKGREEN)
+        printc(f"------------ Load flats time: {np.round(time.time() - start_time,3)} seconds",bcolors.OKGREEN)
+        printc('--------------------------------------------------------------',bcolors.OKGREEN)
 
-        except Exception:
-            printc("ERROR, Unable to open flats file: {}",flat_f,color=bcolors.FAIL)
+        # except Exception:
+        #     printc("ERROR, Unable to open/process flats file: {}",flat_f,color=bcolors.FAIL)
 
 
     else:
@@ -451,8 +474,6 @@ def phihrt_pipe(data_f, dark_f = '', flat_f = '', scale_data = True, norm_f = Tr
 
         #demod the flats
 
-        
-
         flat_demod, demodM = demod_hrt(flat, flat_pmp_temp)
 
         norm_factor = norm_factor = np.mean(flat_demod[512:1536,512:1536,0,0])
@@ -640,6 +661,9 @@ def phihrt_pipe(data_f, dark_f = '', flat_f = '', scale_data = True, norm_f = Tr
 
         field_stop = np.where(field_stop > 0,1,0)
 
+        if imgdirx_flipped == 'NO':
+            field_stop = field_stop[:,::-1]
+
         data *= field_stop[start_row:start_row + data_size[0],start_col:start_col + data_size[1],np.newaxis, np.newaxis, np.newaxis]
 
         printc('--------------------------------------------------------------',bcolors.OKGREEN)
@@ -795,7 +819,7 @@ def phihrt_pipe(data_f, dark_f = '', flat_f = '', scale_data = True, norm_f = Tr
 
         #check if the output directory exists, if not, create it
         if not os.path.exists(out_dir): 
-            print(f"{out_dir} does not exist, creating it now")
+            print(f"{out_dir} does not exist -->>>>>>> Creating it")
             os.makedirs(out_dir)
         
         if isinstance(data_f, list):
@@ -811,7 +835,7 @@ def phihrt_pipe(data_f, dark_f = '', flat_f = '', scale_data = True, norm_f = Tr
             #print(uniq_scan_DIDs)
             #print(scan_name_list)
             if uniq_scan_DIDs == []:
-                print("The scan's DIDs are all unique")
+                print("The scans' DIDs are all unique")
 
             else:
 
@@ -891,8 +915,11 @@ def phihrt_pipe(data_f, dark_f = '', flat_f = '', scale_data = True, norm_f = Tr
             #must invert each scan independently, as cmilos only takes in one dataset at a time
 
              #get wave_axis from the header information of the science scans
+            if cpos_arr[0] == 0:
+                shift_w =  wave_axis[3] - wavelength
+            elif cpos_arr[0] == 5:
+                shift_w =  wave_axis[2] - wavelength
 
-            shift_w =  wave_axis[3] - wavelength
             wave_axis = wave_axis - shift_w
 
             print('It is assumed the wavelength array is given by the header')

run.py

 from hrt_pipe import phihrt_pipe
 import numpy as np
 
-sciencedata_fits_filenames = ['solo_L0_phi-hrt-ilam_20210421T120003_V202106080929C_0144210101.fits', 'solo_L0_phi-hrt-ilam_20210424T120003_V202106141014C_0144240101.fits',
-  'solo_L0_phi-hrt-ilam_20210425T120002_V202106141020C_0144250101.fits', 'solo_L0_phi-hrt-ilam_20210426T120002_V202106162118C_0144260101.fits',
-  'solo_L0_phi-hrt-ilam_20210427T120002_V202106162052C_0144270101.fits', 'solo_L0_phi-hrt-ilam_20210427T120002_V202106171444C_0144270101.fits', 
-  'solo_L0_phi-hrt-ilam_20210427T120002_V202106171517C_0144270101.fits'] ##['solo_L1_phi-hrt-ilam_20210223T170002_V202106111612C_0142230201.fits']
+sciencedata_fits_filenames = ['solo_L1_phi-hrt-ilam_20200528T171109_V202106111600C_0045140102.fits']
+
+# ['solo_L0_phi-hrt-ilam_20210421T120003_V202106080929C_0144210101.fits', 'solo_L0_phi-hrt-ilam_20210424T120003_V202106141014C_0144240101.fits',
+#   'solo_L0_phi-hrt-ilam_20210425T120002_V202106141020C_0144250101.fits', 'solo_L0_phi-hrt-ilam_20210426T120002_V202106162118C_0144260101.fits',
+#   'solo_L0_phi-hrt-ilam_20210427T120002_V202106162052C_0144270101.fits', 'solo_L0_phi-hrt-ilam_20210427T120002_V202106171444C_0144270101.fits', 
+#   'solo_L0_phi-hrt-ilam_20210427T120002_V202106171517C_0144270101.fits'] ##['solo_L1_phi-hrt-ilam_20210223T170002_V202106111612C_0142230201.fits']
   
 #['solo_L0_phi-hrt-ilam_0667414748_V202103221851C_0142230201.fits']#['solo_L0_phi-hrt-ilam_0667414905_V202103221851C_0142230602.fits', 'solo_L0_phi-hrt-ilam_0667415054_V202103221851C_0142230603.fits', 'solo_L0_phi-hrt-ilam_0667415205_V202103221851C_0142230604.fits', 'solo_L0_phi-hrt-ilam_0667415354_V202103221851C_0142230605.fits', 'solo_L0_phi-hrt-ilam_0667415505_V202103221851C_0142230606.fits', 'solo_L0_phi-hrt-ilam_0667415654_V202103221851C_0142230607.fits', 'solo_L0_phi-hrt-ilam_0667415805_V202103221851C_0142230608.fits']#['../fits_files/solo_L0_phi-hrt-ilam_0667414748_V202103221851C_0142230201.fits']
 #sciencedata_fits_filenames = ['solo_L0_phi-hrt-ilam_0667414748_V202103221851C_0142230201.fits']
 #sciencedata_fits_filenames = ['solo_L0_phi-hrt-ilam_0667414905_V202103221851C_0142230602.fits', 'solo_L0_phi-hrt-ilam_0667415054_V202103221851C_0142230603.fits', 'solo_L0_phi-hrt-ilam_0667415205_V202103221851C_0142230604.fits', 'solo_L0_phi-hrt-ilam_0667415354_V202103221851C_0142230605.fits', 'solo_L0_phi-hrt-ilam_0667415505_V202103221851C_0142230606.fits', 'solo_L0_phi-hrt-ilam_0667415654_V202103221851C_0142230607.fits', 'solo_L0_phi-hrt-ilam_0667415805_V202103221851C_0142230608.fits']
 
-flatfield_fits_filename = '../fits_files/solo_L0_phi-hrt-flat_0667134081_V202103221851C_0162201100.fits'
+flatfield_fits_filename = '/data/slam/home/sinjan/fits_files/solo_L1_phi-hrt-ilam_20200417T174538_V202106111549C_0024150020.fits'#solo_L0_phi-hrt-flat_0667134081_V202103221851C_0162201100.fits'
 darkfield_fits_filename = '../fits_files/solo_L0_phi-fdt-ilam_20200228T155100_V202002281636_0022210004_000.fits'
 
-sciencedata_fits_filenames = ['../fits_files/' + i for i in sciencedata_fits_filenames]
+sciencedata_fits_filenames = ['/data/slam/home/sinjan/fits_files/' + i for i in sciencedata_fits_filenames]
 
 prefilter_f = '../fits_files/fitted_prefilter.fits'
 
@@ -27,13 +29,13 @@ prefilter_f = '../fits_files/fitted_prefilter.fits'
 
 c_talk_params = np.zeros((2,3))
 
-q_slope = 0.0038
-u_slope = -0.0077
-v_slope = -0.0009
+q_slope = -0.0098#0.0038
+u_slope = -0.0003#-0.0077
+v_slope = -0.0070#-0.0009
 
-q_int = -0.0056 #the offset, normalised to I_c
-u_int = 0.0031
-v_int = -0.0002 
+q_int = -0.0015#-0.0056 #the offset, normalised to I_c
+u_int = 0.0007#0.0031
+v_int = 0.0006#-0.0002 
 
 c_talk_params[0,0] = q_slope
 c_talk_params[0,1] = u_slope
@@ -43,10 +45,10 @@ c_talk_params[1,0] = q_int
 c_talk_params[1,1] = u_int
 c_talk_params[1,2] = v_int
 
-phihrt_pipe(sciencedata_fits_filenames, flat_f = flatfield_fits_filename, dark_f = darkfield_fits_filename, scale_data = True, norm_f = True, clean_f = True, 
-            sigma = 59, flat_states = 24, norm_stokes = True, prefilter_f = prefilter_f, dark_c = True, flat_c = True, 
+phihrt_pipe(sciencedata_fits_filenames, flat_f = flatfield_fits_filename, dark_f = darkfield_fits_filename, scale_data = True, bit_flat = False, norm_f = True, clean_f = False, 
+            sigma = 59, flat_states = 24, norm_stokes = True, prefilter_f = None, dark_c = True, flat_c = True, 
             fs_c = True, demod = True, ctalk_params = c_talk_params, ItoQUV = True, out_demod_file = True, 
-            out_dir = '/path/to/files', rte = 'RTE', out_rte_filename='your_desired_rte_filename') 
+            out_dir = '/data/slam/home/sinjan/hrt_pipe_results/may_2020_ctalk/', rte = 'False', out_rte_filename='') 
 """
  Input Parameters:
 ----------

utils.py

@@ -62,30 +62,32 @@ def fits_get_sampling(file,verbose = False):
             tunning_constant = 0.0
             ref_wavelength = 0.0
             for v in header:
+                #print(v)
                 if (j < 6):
                     if tunning_constant == 0:
                         tunning_constant = float(v[4])/1e9
                     if ref_wavelength == 0:
                         ref_wavelength = float(v[5])/1e3
-                    if np.abs(np.abs(v[2]) - np.abs(dummy)) > 5:
+                    if np.abs(np.abs(v[2]) - np.abs(dummy)) > 5: #check that the next voltage is more than 5 from the previous, as voltages change slightly
+                        #print(dummy, v[2])
                         voltagesData[j] = float(v[2])
                         dummy = voltagesData[j] 
                         j += 1
-        
-        d1 = voltagesData[0] - voltagesData[1]
-        d2 = voltagesData[4] - voltagesData[5]
-        if np.abs(d1) > np.abs(d2):
-            cpos = 0
-        else:
-            cpos = 5
-        if verbose:
-            print('Continuum position at wave: ', cpos)
-        wave_axis = voltagesData*tunning_constant + ref_wavelength  #6173.3356
-
-        return wave_axis,voltagesData,tunning_constant,cpos
 
     except Exception:
-        print("Unable to open fits file: {}",file)     
+       print("Unable to open fits file: {}",file)     
+    #print(voltagesData)
+    d1 = voltagesData[0] - voltagesData[1]
+    d2 = voltagesData[4] - voltagesData[5]
+    if np.abs(d1) > np.abs(d2):
+        cpos = 0
+    else:
+        cpos = 5
+    if verbose:
+        print('Continuum position at wave: ', cpos)
+    wave_axis = voltagesData*tunning_constant + ref_wavelength  #6173.3356
+    #print(wave_axis)
+    return wave_axis,voltagesData,tunning_constant,cpos
 
 def fix_path(path,dir='forward',verbose=False):
     path = repr(path)