debug smooth_data call with 'arcsec' scale and reduce 3C405

src/FOC_reduction.py

@@ -29,10 +29,10 @@ def main():
 #            'x14w0104t_c1f.fits','x14w0105p_c1f.fits','x14w0106t_c1f.fits']
 #    globals()['plots_folder'] = "../plots/NGC1068_x14w010/"
 
-#    globals()['data_folder'] = "../data/3C405_x136060/"
+    globals()['data_folder'] = "../data/3C405_x136060/"
-#    infiles = ['x1360601t_c0f.fits','x1360602t_c0f.fits','x1360603t_c0f.fits']
+    infiles = ['x1360601t_c0f.fits','x1360602t_c0f.fits','x1360603t_c0f.fits']
-#    infiles = ['x1360601t_c1f.fits','x1360602t_c1f.fits','x1360603t_c1f.fits']
+    infiles = ['x1360601t_c1f.fits','x1360602t_c1f.fits','x1360603t_c1f.fits']
-#    globals()['plots_folder'] = "../plots/3C405_x136060/"
+    globals()['plots_folder'] = "../plots/3C405_x136060/"
 
 #    globals()['data_folder'] = "../data/CygnusA_x43w0/"
 #    infiles = ['x43w0101r_c0f.fits', 'x43w0104r_c0f.fits', 'x43w0107r_c0f.fits',
@@ -46,14 +46,6 @@ def main():
 #    infiles = ['x3mc0101m_c0f.fits','x3mc0102m_c0f.fits','x3mc0103m_c0f.fits']
 #    globals()['plots_folder'] = "../plots/3C109_x3mc010/"
 
-    globals()['data_folder'] = "../data/PGN1630377_x39510/"
-    infiles = ['x3990201m_c0f.fits', 'x3990205m_c0f.fits', 'x3995101r_c0f.fits',
-            'x3995105r_c0f.fits', 'x3995109r_c0f.fits', 'x3995201r_c0f.fits',
-            'x3995205r_c0f.fits', 'x3990202m_c0f.fits', 'x3990206m_c0f.fits',
-            'x3995102r_c0f.fits', 'x3995106r_c0f.fits', 'x399510ar_c0f.fits',
-            'x3995202r_c0f.fits','x3995206r_c0f.fits']
-    globals()['plots_folder'] = "../plots/PGN1630377_x39510/"
-
 #    globals()['data_folder'] = "../data/MKN463_x2rp030/"
 #    infiles = ['x2rp0201t_c0f.fits', 'x2rp0203t_c0f.fits', 'x2rp0205t_c0f.fits',
 #            'x2rp0207t_c0f.fits', 'x2rp0302t_c0f.fits', 'x2rp0304t_c0f.fits',
@@ -62,6 +54,14 @@ def main():
 #            'x2rp0305t_c0f.fits', 'x2rp0307t_c0f.fits']
 #    globals()['plots_folder'] = "../plots/MKN463_x2rp030/"
 
+#    globals()['data_folder'] = "../data/PG1630+377_x39510/"
+#    infiles = ['x3990201m_c0f.fits', 'x3990205m_c0f.fits', 'x3995101r_c0f.fits',
+#            'x3995105r_c0f.fits', 'x3995109r_c0f.fits', 'x3995201r_c0f.fits',
+#            'x3995205r_c0f.fits', 'x3990202m_c0f.fits', 'x3990206m_c0f.fits',
+#            'x3995102r_c0f.fits', 'x3995106r_c0f.fits', 'x399510ar_c0f.fits',
+#            'x3995202r_c0f.fits','x3995206r_c0f.fits']
+#    globals()['plots_folder'] = "../plots/PG1630+377_x39510/"
+
     ## Reduction parameters
     # Deconvolution
     deconvolve = False
@@ -70,28 +70,28 @@ def main():
         psf_FWHM = 0.10
         psf_scale = 'arcsec'
         psf_shape=(9,9)
-        iterations = 5
+        iterations = 10
     # Error estimation
-    error_sub_shape = (100,100)
+    error_sub_shape = (200,200)
-    display_error = True
+    display_error = False
     # Data binning
     rebin = True
     if rebin:
-        pxsize = 0.10
+        pxsize = 0.35
         px_scale = 'arcsec'         #pixel or arcsec
         rebin_operation = 'sum'     #sum or average
     # Alignement
     align_center = 'image'        #If None will align image to image center
-    display_data = True
+    display_data = False
     # Smoothing
     smoothing_function = 'combine'  #gaussian or combine
     smoothing_FWHM = None           #If None, no smoothing is done
-    smoothing_scale = 'pixel'       #pixel or arcsec
+    smoothing_scale = 'arcsec'       #pixel or arcsec
     # Rotation
     rotate = False                  #rotation to North convention can give erroneous results
     rotate_library = 'scipy'        #scipy or pillow
     # Polarization map output
-    figname = 'PGN1630+377_FOC'         #target/intrument name
+    figname = '3C405_FOC'         #target/intrument name
     figtype = ''    #additionnal informations
     SNRp_cut = 3    #P measurments with SNR>3
     SNRi_cut = 30   #I measurments with SNR>30, which implies an uncertainty in P of 4.7%.

src/lib/reduction.py

@@ -716,7 +716,7 @@ def align_data(data_array, error_array=None, upsample_factor=1., ref_data=None,
         return rescaled_image, rescaled_error
 
 
-def smooth_data(data_array, error_array, FWHM=1., scale='pixel',
+def smooth_data(data_array, error_array, headers, FWHM=1., scale='pixel',
         smoothing='gaussian'):
     """
     Smooth a data_array using selected function.
@@ -727,6 +727,8 @@ def smooth_data(data_array, error_array, FWHM=1., scale='pixel',
     error_array : numpy.ndarray
         Array of images (2D floats, aligned and of the same shape) containing
         the error in each pixel of the observation images in data_array.
+    headers : header list
+        List of headers corresponding to the images in data_array.
     FWHM : float, optional
         Full Width at Half Maximum for desired smoothing in 'scale' units.
         Defaults to 1.
@@ -887,14 +889,18 @@ def polarizer_avg(data_array, error_array, headers, FWHM=None, scale='pixel',
         err60_array = error_array[is_pol60]
         err120_array = error_array[is_pol120]
 
+        headers0 = [header for header in headers if header['filtnam1']=='POL0']
+        headers60 = [header for header in headers if header['filtnam1']=='POL60']
+        headers120 = [header for header in headers if header['filtnam1']=='POL120']
+
         if not(FWHM is None) and (smoothing.lower() in ['combine','combining']):
             # Smooth by combining each polarizer images
-            pol0, err0 = smooth_data(pol0_array, err0_array, FWHM=FWHM,
+            pol0, err0 = smooth_data(pol0_array, err0_array, headers0,
-                    scale=scale, smoothing=smoothing)
+                    FWHM=FWHM, scale=scale, smoothing=smoothing)
-            pol60, err60 = smooth_data(pol60_array, err60_array, FWHM=FWHM,
+            pol60, err60 = smooth_data(pol60_array, err60_array, headers60,
-                    scale=scale, smoothing=smoothing)
+                    FWHM=FWHM, scale=scale, smoothing=smoothing)
-            pol120, err120 = smooth_data(pol120_array, err120_array, FWHM=FWHM,
+            pol120, err120 = smooth_data(pol120_array, err120_array, headers120,
-                    scale=scale, smoothing=smoothing)
+                    FWHM=FWHM, scale=scale, smoothing=smoothing)
 
         else:
             # Average on each polarization filter.
@@ -908,10 +914,12 @@ def polarizer_avg(data_array, error_array, headers, FWHM=None, scale='pixel',
             err60 = np.mean(err60_array,axis=0)/np.sqrt(err60_array.shape[0])
             err120 = np.mean(err120_array,axis=0)/np.sqrt(err120_array.shape[0])
             polerr_array = np.array([err0, err60, err120])
+
+            headers_array = headers0 + headers60 + headers120
             if not(FWHM is None):
                 # Smooth by convoluting with a gaussian each polX image.
                 pol_array, polerr_array = smooth_data(pol_array, polerr_array,
-                        FWHM=FWHM, scale=scale)
+                        headers_array, FWHM=FWHM, scale=scale)
                 pol0, pol60, pol120 = pol_array
                 err0, err60, err120 = polerr_array