# # This file is part of Herschel Common Science System (HCSS). # Copyright 2001-2011 Herschel Science Ground Segment Consortium # # HCSS is free software: you can redistribute it and/or modify # it under the terms of the GNU Lesser General Public License as # published by the Free Software Foundation, either version 3 of # the License, or (at your option) any later version. # # HCSS is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU Lesser General Public License for more details. # # You should have received a copy of the GNU Lesser General # Public License along with HCSS. # If not, see . # # This script performs a loop over a list of obsid. The loop calls the script provided to reduce # individual obsid for the deep survey/ mini scan map case (scanmap_Deep_Survey_miniscan_Pointsource.py). # This script processes scan map and mini scan map observations from Level 0 (raw data) up to Level 2 # (final map). There are two flavors for the map projection: the highpass filter way and the MadMap # algorithm. This script is using the highpass filter way to remove the 1/f noise.In order to # avoid mistakes in the parameter settings, we suggest to put all the settings in this script. # For deep survey and mini scam map observation of faint sources, we recommend to create the # mask for the masked High-pass filtering after combining all the obsid maps of the observation. # Indeed, only in the final mosaic you can reach the required S/N to detect and properly mask the # faint sources. ###################### SECTION 0: SETTINGS ####################################################### # If the input_setting variable is set to True, the script will read the settings from this script # and will ingnore any other setting given in the individual obsid script. If it is set to False, # the script will read the settings in the Section 0 of the individual obsid script. # input_setting=True # directory name where to read or save files. direc = "/Users/baltieri/work/HSCworkshopOT1/" # By default, the reduction script for the individual observations will be taken from the scripts directory # in the hipe distribution. If you want to use your edited copy of the reduction script for the individual # observations, comment the next line and define scriptsDir with the full path name to the directory # where you have saved your edited copy. #scriptsDir = Configuration.getProperty("var.hcss.dir") + "/scripts/pacs/spg/pipeline/ipipe/phot/" scriptsDir = "/Users/baltieri/work/HSCworkshopOT1/" # the user can provide directly a list of obsid number in this way # obsidall=[1342189191,1342189192] #obsidall=[1342195466,1342195467] # #or the list of obsid can be read from a ascii file in this way # #splitfilesname='your_obsid_list' #infile=open(splitfilesname,'r') #obsidall=infile.readlines() #infile.close() #for i in range(len(obsidall)): # obsidall[i]=obsidall[i].strip() # here we set the channel to be reduced camera = "red" #camera = "blue" # setting of output map pixel and drop size: the choice of the output pixel and drop size depends on the # redundancy of your observation. See the PACS Data Reduction Guide or the PACS User Manual for more # details about the choice of those parameters. Several tips are provided also in the individual obsid # scripts. Here we set the output pixel size to 2 arcsec for the blue channel and 3 arcsec for the red # channel. These parameter values are found to work generally well with a drop size (pixfrac) of 1/10 the # input pixel size. Anyhow, we invite the user to play with these values to find the optimal combination # for the specific scientific case. # if camera=='blue': outpixsz=2. elif camera=='red': outpixsz=2. pixfrac=1.0 # A TIP: you might want to use a specific wcs for your final map. Here it is shown how to set the wcs # parameters. The user have to set the sky coordinates of the reference pixel (ra_reference and dec_reference) # which is supposed to be at the center of the map, and the half dimension of the map in unit of pixels # (rad1 and rad2) for the given output pixel size (outpixsz set above). To set these parameters the user # must have an idea of the map dimensions for the given output pixel size. You can get the required values # by looking at the Level 2 map provided by the automatic pipeline. Setting a fixed wcs could be very # useful if your observation comprises many obsid and you want to produce maps all with the same wcs. # Indeed, in this case you can combine the individual obsid maps easily to get the final map by just # averaging or by taking a weighted mean (e.g. weighted by coverage) of the individual output pixels. # In addition the standard deviation of the output pixel signal can provide a good estimate of the noise. # In this way you can obtain a meaningful noise map associated to your final map, which is still not # available as output of the current pipeline version. # #ra_reference=0.000 #provide here the map center ra in deg #dec_reference=0.000 #provide here the map center dec in deg ## here we build the wcs # #pixsize=outpixsz #rad1=250 ##(half the dimension of the map in the x direction in pixel units) #rad2=250 ##(half the dimension of the map in the y direction in pixel units) # #my_wcs=Wcs(cunit1="Degrees",cunit2="Degrees",cdelt1=-pixsize/3600.,\ # cdelt2=pixsize/3600.,crota2=0.,crpix1=rad1,crpix2=rad2,\ # crval1=ra_reference,crval2=dec_reference,ctype1="RA---TAN",ctype2="DEC--TAN",\ # equinox=2000.0) #my_wcs.setParameter("naxis1",2*rad1,"naxis1") #my_wcs.setParameter("naxis2",2*rad2,"naxis2") # #### HP RADIUS AND MASK SETTINGS FOR DEEP SURVEY/MINI SCAN MAP CASE ######################################### # # setting of the highpass filter radius in number of readouts. The following numbers are suited for point # sources. Larger radii should be set for very bright or extended sources with the caveat that the 1/f # noise will not be properly removed. As a general rule, the smaller the high pass filter radius, the # better you remove the 1/f noise. The suggested values allow to remove as much as possible the 1/f without # damaging the PSF. Alternative values have to be chosen on the basis of the source dimension. As a rule # of thumb the hpfradius should be as large as the source (mask) size. However, if the interpolation method # is used in the high-pass filter task, hpfradius can be smaller then the source size. In general, values # larger that 100 readouts are not reccomended. Several tips about the choise of this parameter can be # found in the scanmap_Deep_Survey_miniscan_Pointsource.py script. # if camera=='blue': hpfradius=15 elif camera=='red': hpfradius=25 # # The individual obsid script provides several options for masking the sources in the high pass filtering # task. # If your scientific case is a deep survey of a black field, your sources will not be detectable on the # single obsid map due to the low S/N. We suggest to perform a preliminary reduction without masking # any source by using the option =1 . You need to reduce all obsid, combine the individual maps and base # the mask on a first preliminary mosaic (see below for the details how to make the mosaic). # Although the preliminary map will show strong high-pass filtering residuals, you can use it to obtain a good # mask of the sources. You can, then, reduce once again the individual obsid by using the mask based on the much # deeper mosaic by using option = 2. # # If your observation containes relatively bright point sources and comprises only two or a small number of obsid, # you probably have enough S/N to detect the sources in the individual obsid maps. So you can create the mask # directly from there. In this case we suggest to choose the option = 3 (see details in the # scanmap_Deep_Survey_miniscan_Pointsource.py script). # # We set by default the option =1, which is providing the first preliminary reduction. Edit the script if you prefer # to use another option. # option = 1 ############### END OF SECTION 0 ################################################################################# ################################################################################################################## # # # The loop is starting. for i in range(len(obsidall)): obsid= obsidall[i] print "Reducing OBSID:", obsid, "camera:", camera execfile(scriptsDir + 'scanmap_Deep_Survey_miniscan_Pointsource.py') # # At this point, a preliminary map still with high-pass filtering residuals is available for any # obsid. If you did not provide a common wcs for all maps, you can simply combine them by using the # MakeMosaicTask to obtain the final map as in the following few lines of code. The script assumes that you # used the same name convention adopted in the individual obsid script for saving the maps in fits files. # If you used a different name convention, please, edit the script and change the map_file variable that # contains the map file names. from java.util import ArrayList from herschel.ia.toolbox.image import MosaicTask images=ArrayList() for i in range(len(obsidall)): obsid=obsidall[i] map_file = direc+ "map_"+"_" + str(obsid) + "_" + camera + ".fits" map=simpleFitsReader(map_file) map.exposure=map.coverage images.add(map) # the variable "images" contains all the individual obsid maps you want to combine. The mosaic is done in # the following way: # mosaic=MosaicTask()(images=images,oversample=0) Display(mosaic) # AN IMPORTANT TIP: The MakeMosaic task relies on the error maps provided by the pipeline to estimate the # error map of the final mosaic. We stress here that the error propagation of the PACS pipeline is still # under construction and that the current available error maps produced by the pipeline are not reliable. # To create a reliable error map in the final mosaic we suggest to project all the individual obsid maps # with the same wcs. The way to set the wcs is suggested in the SECTION 0 of this script and how to force # PhotProject to use the desired wcs is described in the individual obsid script. Once the maps are created # with the same wcs, the final map can be obtained by taking a simple average or weighted mean # (weighted by coverage) of the individual output pixels. The standard deviation of the mean or weighted mean # can provide a meaningful final error map. If the number of obsid is not high enough to let you estimate # the mean (weighted mean) and standard deviation, you can slice the single obsid frames per repetition # or odd and even scan legs to produce a sufficient number of individual maps to accurately estimate # mean (weighted mean) and standard deviation. # If you chose the option = 3 and you are now happy with the result, you can simply save the final # mosaic and skip the rest of the script. #outfile_name = direc+ "final_map"+ "_" + camera + ".fits" #print "Saving file: " + outfile_name #simpleFitsWriter(mosaic,outfile_name) # If you chose the option =1, you just did a quicklook data reduction. Thus, you need to proceed to improve # the data reduction with the following scheme. You can use the preliminary and deep mosaic to detect also # the faint sources and create the mask. You need, then, to re-reduce the data, starting from Level 1, and # used the new mask to perform the masked high-pass filtering and projection. # We propose here the same method used in the single obsid script for masking the sources. Be aware that # here we use the exposure and not the coverage because the MakeMosaicTask produce only the exposure map # with the same meaning of the coverage map. This method follows the following steps: # # 1) identifying the region of the map with high coverage: this is done by taking the region where the # coverage # is above 1sigma # med=STDDEV(mosaic.exposure[mosaic.exposure.where(mosaic.exposure > 0.)]) # # 2) using this region to estimate the signal standard deviation (stdev) # signal_stdev=STDDEV(mosaic.image[mosaic.image.where(mosaic.exposure > med)]) # # 3) defining the threshold as cutlevel*stdev: we set the cutlevel equal to 3.0 by default. # cutlevel=3.0 threshold=cutlevel*signal_stdev # # 4) masking everything above the threshold # mask=mosaic.copy() mask.image[mask.image.where(mosaic.image > threshold)] = 1.0 mask.image[mask.image.where(mosaic.image < threshold)] = 0.0 Display(mask) # you can save the mask in a fits file to be used for any further re-processing maskfile = direc+ "mask"+ "_" + camera + ".fits" print "Saving file: " + maskfile simpleFitsWriter(mask,maskfile) # if you commented properly the SECTION 0 of the individual obsid script, the maskfile name will be passed # to the script for the masked high-pass filtering # A TIP: the cut level to use in order to identify the sources depends on the scientific case. Generally # a 3sigma level is sufficient for the deep survey/mini scan map case. You can also chose to refine the # method by recalculating iteratively the noise and the threshold by avoiding the masked pixels. Indeed, # if you discard the sources in the map the standard deviation should be less affected by the source flux # and should be smaller. It can also be useful to smooth the image (e.g. with the GaussianSmooth Task) # to avoid masking individual noisy pixels. # If you are satisfied by the mask, then you can re-run the final processing by providing to the # individual obsid script the new source mask by choosing option =2. We suggest to comment the whole # section 0, 1, and 2 in the individual obsid script and restart the data processing from the Level 1. option = 2 for i in range(len(obsidall)): obsid= obsidall[i] print "re-reducing OBSID:", obsid, "camera:", camera execfile(scriptsDir + 'scanmap_Deep_Survey_miniscan_Pointsource.py') # If you did not provide a common wcs for all maps, you can simply combine them by using the MakeMosaicTask # to obtain the final map as in the following few lines of code. The script assumes that you used the same # name convention adopted in the individual obsid script for saving the maps in fits files. If you used a # different name convention, please, edit the script and change the map_file variable that contains the # map file name. from java.util import ArrayList from herschel.ia.toolbox.image import MosaicTask images=ArrayList() for i in range(len(obsidall)): obsid=obsidall[i] map_file = direc+ "map_"+"_" + str(obsid) + "_" + camera + ".fits" map=simpleFitsReader(map_file) map.exposure=map.coverage images.add(map) # mosaic=MosaicTask()(images=images,oversample=0) Display(mosaic) # Now you can finally save your final map in a fits file: outfile_name = direc+ "final_map"+ "_" + camera + ".fits" print "Saving file: " + outfile_name simpleFitsWriter(mosaic,outfile_name)