3.6) The PIA Mapping Software

3.6.1) Data Description

Measurement:

PIA mapping software is specially suited for processing measurements which have been performed in one of the raster mode AOTs (for PHT-P, PHT-C or PHT-S) or with AOT P32 (oversampling with PHT-C). This is not a limitation for producing a map with any kind of PHT observations (eg. staring observations).

Note that it processes one AAP data set, also called an AAP measurement, at a time. In other words: an AAP measurement contains in the normal case data corresponding of one measurement, performed with one filter/ aperture for all sky positions of the raster or chopper directions, pre-processed to the AAP level. Nevertheless, PIA allows for concatenating different measurements creating a "pseudo" AAP measurement, making possible eg. coherent co-addition of different maps, but also of different staring observations, like in the case of sparse mapping.

Format:

The PIA mapping software reads data in the AAP (Auto Analysis Product) format, which can be processed to a map. This data format contains the measurements in units of flux densities (Jy) and brightnesses (MJy/sr). The PIA mapping software offers the possibility of using any of them (including the different computations of averages or medians of signals per chopper plateau) with a default choice of mean brightness.

AAP measurements can exist in two different forms, which are:

(a) PIA internal AAP structure saved on hard-disk
(b) PIA internal AAP structure in the buffer

The PIA mapping software expects AAP data, ready for processing, in the form (b), i.e. as PIA internal AAP structure in the buffer. There are essentially two methods for providing such internal data structures:

1. Raw data, which have been processed by earlier PIA processing steps to the AAP level, are automatically present inside PIA as data type (b). Individual AAP measurements can be selected directly for processing without further action.
2. File types (a) can be read in from disk into the buffer (converted to data type (b)) by activating the proper menu entry in the PIA Top Window: Files - Read External Files.

Hint: Technically, the internal AAP structure is stored in the IDL common block PHTAAP.

Produced Output:

PIA_MAP produces three kinds of maps of each AAP measurement, which are:

1. A brightness/flux map: The map of either brightness/flux mean values, median values, 1st quartile, or 3rd quartile values, depending on the user's selection. The unit of map values is: MJy/sr or Jy/pixel.
2. An uncertainty map: The map of uncertainty values, processed from the uncertainty array in the AAP data. The unit of map values is: MJy/sr or Jy/pixel.
3. An exposure map: The map of exposure times for each map point. The unit of map values is: s (seconds). The exposure time is the total of all exposure times from the contributing detector signals, corrected for the map pixel size.

The maps are visualized on the screen or can be printed out onto paper output. It is possible to write maps to files and to reload them, as well as to perform operations, superpose contours, etc.

3.6.2) The PIA main Mapping Window and associated GUIs

FIGURE 1: The PIA Mapping Window 

The window offers the following menu buttons:

File

Select Measurement:

This opens a pop-up window ( Please Select a Measurement as described in section 3.2.4) listing AAP measurements stored in the buffer. Click with the mouse onto one of them and press the OK button.

The pop-up window should now disappear, and the PIA Mapping Window should fill with some informational items of the selected measurement.

IRAF's Drizzle Batch:

To send a batch job to the IRAF package for co-adding the map using the drizzling algorithm, as explained under section 3.6.4.
Read Fits: PHT Type

This reads PHT maps produced and saved by either PIA or the PHT-pipeline. Also for reading every default FITS map (stored as a primary data array).

Read Fits: CAM Type

This reads CAM images as produced by the CAM-pipeline (CMAP*.FITS). They are stored as binary extension FITS files.

Read Internal files:

This reads maps which have already been produced and saved in internal (IDL/XDR) format (Filename *.IMAP).

Quit:

Leaves the PIA Mapping Window. It will close all the displayed maps and re-initialize the color table to the PIA color table. Of course, internally stored AAP buffers (=the input data) will be kept in memory.

Process

Define Map Parameters

This opens a window as showed in Figure 2. 

---

FIGURE 2: The Define Map Parameters Window 

In this window the following selections can be made:

1. Control section:
• OK: stores the chosen parameters in the I/F, returning to the calling main map window for further action.
• Cancel: does not change previous chosen parameters by returning to the calling window.
• Do the map: produces immediately the map according to all the selections given in this window, closing it. The process performed is the same as started by the button Produce Map in the PIA Mapping window, which is described below. At the same time storesthe parameters and return to the calling window.
  
  
2. Mapping method: Gives the user the possibility of choosing among different mapping techniques. They are described under the next sub-section: Produce Map.
   
   
3. Pixel size of map in arcsec: Lets you select the size of the resulting map points in X and Y directions. X points on the screen to the right, Y points upwards. X and Y values can be set to different values. The map-pixel size has a major influence on the processing time. Different values for X and Y map-pixel sizes will be transformed correctly to screen-pixel sizes, so that the resulting view of the map is "aspect-free"


Note: The entries `Pixel size of map' are only meaningful for the regular map processing and not for the "simple"  maps (see below).



4. Rotation angle of map relative to RA,DEC: This angle is the orientation of the resulting map X- and Y-coordinates relative to sky coordinates RA and DEC, counter-clockwise. An angle of 0 deg means that the +X coordinate is aligned with the -RA coordinate and the +Y coordinate is aligned with the +DEC coordinate.


Note: the rotation angle is not user configurable, thus giving only the information. The angle is determined depending on the planning of the observation as follows:



A) if NO precise angle was required by the observer at the proposal time (keyword ATTRORIE=1): in this case the rotation angle is simply given by the spacecraft's roll angle, which is taken as the average of the roll angles of all data points during the raster, since in this case the raster Z-axis is coincident with the Y-spacecraft axis.


B) if the observer had required a determined orientation of the raster wrt RA-DEC (keyword ATTRORIE=0): in this case the rotation angle quoted is the required one plus the average of the roll angles of all data points.


5. Set default values: Pressing this button sets the following default values:
• Pixels sizes in X and Y direction are set to default values. They are calculated trying to match the resolution achieved taking into account the raster steps in both directions (and the chopper increment for the P32 case).
• The input data selection is reset to: mean values
  
  
6. Select detector pixel: Lets you select or de-select individual detector pixels for PHT-C. All data from pixels which are de-selected are treated as invalid and do not contribute to the map. For test purposes it might be useful to map only the data from one pixel.
   
   
7. Flat Fielding: There are four flat-fielding procedures which can be enabled,
• Central overlapping Flat Fielding: in this case individual C-detector pixel maps will be created, and their central parts integrated to obtain flat fielding parameters. This is only possible if the image pixels taken have a complete overlap, e.g. no empty image pixels for any one of the individual maps. It also requires at least nine image pixels in each direction.
• Median filter + smooth: this includes a flat fielding method explained below. It is embedded in a menu procedure "PIA mapping median filtering, smoothing and flat-fielding", which is used for selecting/deselecting data, processing it, looking at the results, etc.
• Chosen Background normalization: by this method the user can choose portions of data for each pixel to be used for normalization among the pixels. Since the background portion should be normally flat for all pixels it is better to use only those portions. (If nothing is deselected the integrals of the pixel flux distributions are used for normalization). It uses practically the same GUI as the Median filter + smooth method
• First Quartile normalization: a simple flat-fielding method, by which the computed 1st quartiles corresponding to each detector pixel QI (as representative for a measured background) are used for normalization. Flux/surface brightness of a detector pixel FCORR-I = FI *  <Q> /  QI

Produce Map:

This produces a map using the full pointing information: RA, DEC and, in particular, the roll angle for each raster point.

It processes the AAP data according to the pointing information and the definition of map parameters (above). It produces the brightness map, the uncertainty map and the exposure map. The brightness map will automatically be shown on the screen. The uncertainty and exposure maps can be brought to the screen manually with the menu entries "Show - Map -..." as explained below.

The procedure uses full pointing information as provided in the AAP data structure. It expects three values RA, DEC, and ROLL (right ascension, declination, and the roll angle) for each raster point. Since the AAP data contains one line of data for each chopper step, all values of RA (DEC, ROLL, respectively) must be the same for all chopper plateaux of the same raster point. The procedure calculates from this the positions of the individual detector pixels at different chopper positions.

After calculation of all positions the values of the detector signals are binned into map pixels. For this task PIA offers the user three different internal algorithms as well as the possibility to define a batch job for IRAF to make use of their implemented drizzling method:

A) Full coverage: A simple gridding function is applied, which is a trapezoidal function, i.e. the geometric overlap of detector pixel and map pixel is used as the contributing part of the detector brightness to the map pixel. The algorithm also assumes that the detector pixel is aligned with the map pixel. I.e., for large roll angles and small map pixels, compared to the detector pixel size, the edges of the individual detector pixels are not properly binned into the right map pixels. (Don't be confused by the orientation of the whole detector: the position of the centre of each individual detector pixel is properly calculated from the roll angle.)

B) Distance weighting:  By this method all the different signals overlapping on an image pixel contribute with a weight. The weight is determined using the distance between center of detector pixel contribution and the center of the image pixel DIST, given in units of detector pixel size (or aperture in the case of PHT-P).  The weight for each signal is computed as w = (1-DIST)**3. This method avoids the smoothing produced by the Full Coverage method, at the price of increasing the noise in the map.

C) Trigrid Interpolation: By this method the TRIANGULATE procedure offered by IDL for a Delaunay triangulation of a planar set of irregularly gridded points is used and a subsequent interpolation (IDL's TRIGRID) is performed. A choice among two different interpolations, linear or quintic polynomial can be taken. This method uses the closest neighbours to each point for interpolation. Since in general the PHT maps have a (more or less) large oversampling, previously to the use of this method, a re-sampling of the data is required for a good use of the full information. PIA groups points together which were observed with their centers up to a distance of a quarter of a detector pixel (or aperture size in the case of PHT-P maps), weighting them with their uncertainties.

D) IRAF's Drizzle Method: I/F to run under IRAF the implementation of the so-called "drizzling" method for map co-addition (the formal name for this technique is variable-pixel linear reconstruction). IRAF has to be available and the PIA system variables for the IRAF login directory (!IRAF_LOGIN) and for an IRAF data handling directory (!IRAF_DATA) defined. The original I/F routines have been written by Manfred Stickel (MPIA). Details are given under section 3.6.4. on IRAF's Drizzling I/F.
Produce Map Ignoring Pointing:

This calls simplified mapping routines which ignore true pointing information. The roll angle, in particular, is not accounted for. This means that the map axes are not parallel to RA and DEC.

The map is calculated under the "assumption" that the spacecraft raster builds a perfectly regular grid and the full coverage method is used for binning the map.

To reconstruct the raster, the procedure reads from the header file the number of raster points and raster lines, and the distance of raster points and raster lines. In addition, the central coordinate of the map is read. (Note that the coordinates displayed in the map plots are not correct, because the roll angle is not accounted for!)

The simplified mapping routines can be called for various cases: P32 or PHT-P simple map (staring/chopped): Since the programme knows the measurement parameters from the AAP header, it will call the proper routine for the cases.

In the case of P32 (PHT-C) a further assumption is that the chopper steps, as well as the spacecraft raster points, coincide perfectly with multiples of 1, or 1/2, or 1/3, of the detector pixel size. This assumption is only approximately true: e.g. in a certain oversampling mode the true raster step size in Y direction is 92 arcsec and not 90 arcsec.

As for the regular mapping case, only the brightness map will be shown automatically on the screen. The uncertainty and exposure maps must be brought to the screen manually.

Show

Header: This shows the contents of the AAP header.

Data Arrays This lets you view the data arrays of the AAP structure as line plots. Most of the features are self-explanatory. You can select individual arrays for the X- or the Y-axis. For a two-dimensional array, the index of the first dimension can be selected individually. You also have the option to restrict the plot to certain raster points or raster legs.

Note: The implementation presently expects that you press the RETURN key in the raster legs or raster points field after you have typed in a number there. Otherwise the programme will not detect your changes.

Map:

Brightness, Uncertainty, Exposure(,Undefined) draw the respective maps onto the screen. Of course, the maps must have been produced already. It is possible to keep a window with the map view on the screen, while processing another measurement or processing the same measurement with changed map parameters. This might be useful for comparison purposes. After processing a new map the `Show - Map -...' entries will show the newly processed maps; old maps are lost.  Undefined refers to maps which are read by the PIA mapping GUI, but have no defined units. In this case they can be displayed and saved giving "Arbitrary" as the unit used.

Options

Help

The PIA mapping median filtering, smoothing and flat-fielding window and the PIA flat-fielding by background choice

• baseline drift correction per detector pixel through median filtering and smoothing of the brightness surface distributions +
• the corresponding obtainment of flat-fielding parameters,

OR
• the obtainment of flat-fielding parameters by computing the factors from integrated parts of the different detector pixels' brightness surface distributions.

The algorithms:

In what cases should/could it be used? This method assumes a flat and constant background. If a long term signal drift affects the map, it can correct for it. If gradients are present in the background distribution, it will falsify it.
Point sources are not affected by the correction, while extended structures can be "washed out" if their extension is not much smaller than the width of the median filtering applied. In any case, a careful deselection of the extended source regions for applying median filtering is recommended.
Special care is always adviced for chopped maps, where extended regions and/or neighbouring point sources can be a problem for the recognition of the baseline drift by median filtering.

How does it work? a modelling of the baseline drift of the signal distribution is obtained via a median filtering of variable width fw. This is performed on all the distributions of the different detector pixels, taking only the portions of data which are not de-selected (for the de-selected parts a simple linear interpolation of the median filtered elements is taken). Special care is taken on the edges of the distributions, by applying succesive median filterings using smaller widths fw_i=[fw/2, fw/4 and fw/8] for the first and last fw_i points of the distribution respectively.

The median filtered distributions can then be smoothed via a smoothing algorithm using a variable width sw.

A normalized distribution for every pixel is then calculated by obtaining the quotient of the original distribution divided by the modelled baseline. Re-normalization to the original integral flux of all distributions is then performed, and the corresponding average flat field factors are calculated.

In this case a portion of the flux distributions per detector pixel is taken, normalized to the number of points and ratios between integrations obtained. The ratios obtained can be used as flat field factors, if they represent for every pixel a portion of the background with the same flux and if the distributions are sufficiently flat and not affected by large long term drifts.

The graphical I/Fs:

• Pixel to be used ...: a text window in which the detector pixel can be chosen, which should be used for data selection/de-selection and for showing the results of the processing performed,
• Graphical Selections: from here special PIA_Xplot displays can be used for selecting data portions to be used for the data processing. It includes a reset button (selecting everything) and the possibility of extending the selection made for one pixel to all the others. Please note that the extension is made taking the same points in time and not in sky coordinates, thus care has to be taken that they represent similar features for the different detector pixels.
• Data Processing: both data processing methods can be started from here, labelled Perform baseline correction or Flatfield on background, which are explained below. The parameters used for driving the median filtering and smoothing procedure can be changed by a simple menu started by the button labelled Change params.
• Plotting: original and modelled distributions can be plotted/overplotted for judging the quality of the results. A text window containing the derived flat field parameters can be obtained under Plot/Print Results -> Print derived flatfields.
• Leaving: Pressing the Accept & Quit button returns with the data corrected, the Cancel button leaves the data unchanged.


Please note that median filtering can be performed renouncing to flat fielding by un-choosing the Button labeld  Use derived flatfield for correcting the data.

FIGURE 3: PIA mapping median filtering, smoothing and flat-fielding window.

3.6.3) The PIA Map Display Window and associated GUIs

FIGURE 3: The Brightness Map 

This window has the following options to modify the plot:

Colors:

This allows for several colour settings via standard IDL routines.

Printer:

This calls the PIA Printer Definition Window.

Save Map:

This allows the calculated map to be saved as a FITS file or an internal PIA file.

Custom:

This brings up the customisation window (Set Image Parameters) as shown in Figure 4 to modify the following options:

• Changing the title of the image and character size.
• Enabling/disabling the drawing of coordinate grid lines.
• The maximum number of grid lines per coordinate direction. The programme will optimize this number so that labels appear at nice values.
• Interpolation of map pixels switch on/off. Take care: interpolation might shift brightness peaks for about one map pixel.
• Continuous zoom: enable/disable continuous update in the zooming application.
• Change the aspect ratio. PIA uses by default the image pixel ratio, unless this would give a strange appearance.
• Set cut values: to optimize colour distribution.
• Set number and thickness of contours with different options.


The plot is refreshed after leaving the customisation window (with ok). 

---

FIGURE 4: The Map Customization Window 

---

Redraw:

This refreshes the plot explicitly.

Smooth map:

Applies the IDL procedure SMOOTH to the image and draw the smoothed image into a new mapping GUI.

Deconvolve:

Deconvolution techniques can be applied to the obtained maps, by using maximum Likelihood or maximum Entropy methods in conjunction with the footprint matrices established for all the different PHT filters (except for the P1 detector filters). PIA has a GUI for deconvolution (Fig. 5a), from which the different methods can be chosen, the footprint matrix to be used displayed (Fig. 5b), parameters for the deconvolution procedure changed (Fig.5c) and the results visualized in a normal PIA mapping window (like in Fig. 3).

The footprint matrix is obtained by folding the point spread function with the optical aperture. Footprint matrices are kept in the Cal G files ##FOOTPR.FITS, with ## being the identifier per detector.

The deconvolution parameters are:
* Intrinsic noise: necessary for calculating the Chi-square of the approximation to the image with the footprint matrix. It is computed from the image pixel to pixel variation.
* Chi-square limit: a limit in the minimization of Chi-square, which ends the iterations when reached.
* Max. iterations: a number of iterations which should be maximally applied.
* Chi-square difference limit: for the iterations to be stopped by a minimal change of the chi-squares of successive iterations.
* Subtraction of constant: a constant background can be subtracted to the full image, for optimization of the algorithm.
* Use background constant: the background to be used, if to be subtracted. Calculated from the image.

FIGURE 5a: The Map Deconvolution Menu

---

FIGURE 5b: The footprint matrix as obtainable from the deconvolution menu.

---

FIGURE 5c: The Map Deconvolution Parameters Menu

---

Convolve:

This option starts a convolution menu for converting the map (raster line) to a given resolution. The menu given in Fig. 6 offers the user the choice of converting to every PHT-filter's resolution, given by the Theta convolution parameter (which can also be written to every value). If it is started with the option Scan, a row of the map can be convolved. In the first case a 2-dimensional gaussian approximation to the PSF is used, while in the second a 1-dimensional gauss function is taken.

---

FIGURE 6: The PIA map/scan convolution menu 

---

Rotate:

A rotation of the map can be obtained via a menu which asks for a desired angle wrt the RA-DEC plane. A new map window is produced, in which all the options are available, except for the rotation itself.

Zoom:

With this option a zoomed portion of the image can be obtained in a separate window. The image can be zoomed by constant factors or continuously, according to the option chosen in the Customize window. Switch the continuous zooming on only on fast computers!
The mouse pointer on the map window gives the center of a new zoomed image, by clicking the lefthand mouse buttton. If the continuous zooming option is on, it is continuously updated by moving the mouse on the image.
With the central mouse button the zooming factor can be changed.
Clicking the righthand mouse button exits the zooming application.

Surface:

This displays the map as a 3-dimensional surface in the PIA Xsurface Window.

Surface (lego):

This displays the map as a 3-dimensional surface 'legoplot' in the PIA Xsurface Window.

Find Point Sources:

Gives the position and intensity of point sources.

Profiles:

This allows a continous profile to be viewed cut along the x and the y axis. The position of the axes is given by the mouse pointer. Clicking the lefthand mouse button interchanges the axes. Clicking the righthand mouse button exits the procedure.

Show Header:

This shows the contents of the MAP header.

Map parameters:

This shows the parameters used for deriving the map.

NaN's -> 0.:

By operations with maps some image pixels can get NaN's as a result (eg. division by a map containing zero points). The problems to handle the display of an image contatining NaN's can be easily avoided by converting the NaN's to 0's.

Extract Flux:

By C100 and C200 maps which have been produced with an image pixel size = detector pixel size (46"x46"  or  92"x92" respectively) point source fluxes can be determined. See section 3.6.5 for details.

Quit:

This leaves the map display window.
 

Row below the image display:

Display of positions and contents:

A mouse click on any point within the displayed map will give RA and DEC position together with the value of the image pixel.

Plot profile:

A profile of every map row can be obtained by clicking with the mouse on any point of the map, after activating one of the options under Plot Profile xxx button. A profile of the corresponding map row is then plotted/overplotted on a PIA_Xplot window. If the option -background is chosen, then the median of the profile edges is subtracted from the full profile before plotting it.

3.6.4) IRAF's Drizzling I/F

FIGURE 7: The PIA_IRAF Drizzle I/F menu

IRAF will work on the background for a while to produce the image (it is a time consuming task!). The final data will be put into the !IRAF_DATA defined directory (re-definition through the Customize -> Paths option under the PIA Top Window) as an IRAF file or as a FITS file, depending on the user's choice. If no name is given to the map, the conventional PIA nomenclature will be used for the name of the output file.
The central parameter for the drizzling task is given by the Drizzle drop factor. By drizzling the pixels in the original input images are mapped into pixels in a subsampled output image, but shrinked using the drop factor.

We do not describe further here the parameters which can be chosen for they are sufficiently self-explaining.

3.6.5) Flux extraction from (mini-)maps

The central part of a mini-map is defined by the 3x3 image pixels covering it. Given the flux FLUX  of a point source the observed fluxes FC , FL and FE for the central pixel, the lateral pixels and the corner pixels respectively are given by

        FC = FLUX * FPSF + BCK
            FL = FLUX * F'PSF + BCK
            FE = FLUX * F"PSF + BCK

Here the FPSFs are accounting for the point spread function factors corrections for central, one pixel distance and pixel distance * square root of two and BCK is assumed as a flat background. Solving these equations is easy and delivers background and source flux at the same time.

This procedure is implemented by PIA not limited to mini-maps but extended to all C100 and C200 maps with an image pixel size corresponding to the detector pixel size (indluding the gap between them, eg. 46"x46" for C100 and 92"x92" by C200). So after initiating the flux extraction the user is asked for giving a mouse click on the point source position in the image displayed. After the computation PIA prints the results on a text window, which can be printed or saved.

In order to make a proper error propagation PIA uses the computed statistical uncertainties for each point, if they are present. If the individual uncertainties are not available (eg. because the data was stored as a map without storing the corresponding uncertainty map), an estimation of the mean uncertainty per point is performed, by taking the median of the flux differences between neighbour image pixels in all directions.

This calculation is complemented by a crude estimation of the source flux, using only the central pixel. The background is estimated using the average of all the 8 neighbour pixels, thus neglecting the point spread function factors for the lateral and corner pixels.

Please note that due to the influence of the gaps between the detector pixels, this procedure is only reliable for observations which have been performed using the mini-map strategy (with the proper distances between raster legs and raster points) and well centering on the point source. However, for other cases it can be a good estimation.
 

3.6.6) The 3D Case - Mapping with ISOPHOT-S

If a PHT-S measurement is loaded using the main PIA mapping window, the user will be prompted to select either a detector pixel, or a wavelength region as shown in Figure 8. The corresponding data will be then extracted from the AAP structure chosen and a map can be built, as for any other detector, using all I/Fs described in the former sections.

Once the map is displayed by the PIA  Map Display Window  several options are active under the same menu button used for plotting the profiles, in the lower right corner, as given in Fig. 9. All of them except for the first one, refers to the extraction of  PHT-SL and PHT_SS spectra corresponding to the position chosen in the map by mouse click. The spectra are then plotted into PIA Xplot windows.
 
 

Figure 8: The PHT-S region choice for mapping

The first option ("Plot underlying signals") extracts the signals on the SRD level corresponding to a position in the map chosen by mouse click (if the corresponding phtsrd data is present in the PIA buffers). In this case the data is plotted on the special version of the PIA_Xplot window, which allows the user to de-select data flagging it as bad for a further processing from the SRD level. This facillitates the fast recognition and elimination of glitch residuals, which might be present in the data.
 

Figure 9: The special options for PHT-S data in the PIA Map Display Window