From dshepher@aoc.nrao.edu Wed Dec 12 10:50:41 2001 Date: Wed, 12 Dec 2001 12:26:33 -0500 (EST) From: Debra Shepherd To: aips2-naug@zia.aoc.NRAO.EDU Cc: 'Debra Shepherd ' Subject: aips++ testing notes AIPS++ Testing This script reduces GBT Orion continuum images, averages the 2 big Orion images togther, regrids the Orion images to match a VLA dataset in the inner region, combines the regridded GBT and VLA data using an image feather technique. bugs() found during testing: --------------------------------- 1. image.feather algorithm doesn't regrid the GBT data properly when combining GBT & VLA images. A grid artifact remains in final image. ==> GBT image must be gridded to VLA image specs before input to feather. 2. imagecalc and spaste documentation needs to be clarified 3. regridding the gbt image to the vla size and shape changes the telescope value of the gbt image to VLA and changes the frequency coord val to that of the VLA image. 4. imagecalc destroys beam information. 5. using imager.feather destroys unit & beam information, sets telescope wrong 6. Something that cannot be done before the im.feather function: The scaled GBT image MUST have a VLA telescope set. I was setting it back to GBT and this caused feather to do strange things. When you regridded the GBT image to the size and shape of the VLA image, the telescope in the regridded image was set to VLA. Since you never changed it back to GBT, feather worked correctly for you. Note, the bug that set the telescope in the GBT regridded image to VLA, also set the frequency to the VLA frequency. This corrected a header problem in the GBT data which had the frequency set to 0. Once the regridding bug is fixed and we can combine GBT and VLA data straight, feather won't work unless the scaled GBT image header is changed to have a VLA telescope. Summary: After much testing and developing workarounds, the images were combined sucessfully. Once the bugs are fixed, additional testing must be done to optimize the imaging technique and make sure the flux scale is correct. I will develop a useable script in the near future and put it in the script library (the one below has too many testing twists and notes to follow easily). Debra Detailed Script: ------------------------------------- gbtcalutils.g is the updated script to reduce GBT continuum data that Joe created. Put this in the directory in which you are working. source /aips++/stable/aipsinit.csh /export/data_1/aips/dss/app/orion.gb AIPS++ Version 1.6 Build 363 glish -l gbtcalutils.g NORMAL: gbtcalutils (gc) is ready for use # Look at the available tools in gc = GBT continuum imaging: - field_names(gc) checkms setdata contcal makeimage covercorr plotsource debug # Make the 1st map: gc.setdata('oria1_DCR'); NORMAL: successful readonly open of default-locked table oria1_DCR/NRAO_GBT_GLISH: 28 columns, 41 rows NORMAL: successful readonly open of default-locked table oria1_DCR/NRAO_GBT_IF: 20 columns, 82 rows NORMAL: successful read/write open of default-locked table oria1_DCR/FEED: 12 columns, 1 rows NORMAL: successful read/write open of default-locked table oria1_DCR/POLARIZATION: 4 columns, 1 rows NORMAL: successful read/write open of default-locked table oria1_DCR/SPECTRAL_WINDOW: 14 columns, 1 rows NORMAL: successful readonly open of default-locked table oria1_DCR/POINTING: 15 columns, 4874 rows # This appears to attach the data to the gc tool # do the basic calibration with "receptor 1" data to get T_ANT # from measured voltages with the noise-diode ON and OFF. gc.contcal(); NORMAL: Starting server imager NORMAL: Server started: /aips++/stable/linux/bin/imager (AIPS++ version: 1.6 (build #363)) NORMAL: Opening MeasurementSet oria1_DCR NORMAL: Successfully closed empty server: imager NORMAL: successful read/write open of default-locked table oria1_DCR: 28 columns, 19516 rows # imager is opened and closed ==> CORRECTED & MODEL_DATA columns are created now. # T_CAL default = 1, receptor default = 1, T_ANT now calculated assuming T_CAL=1 # Select, regrid, and construct an "unnormalized" image gc.makeimage('orion1im'); NORMAL: Starting server imager NORMAL: Server started: /aips++/stable/linux/bin/imager (AIPS++ version: 1.6 (build #363)) NORMAL: Opening MeasurementSet oria1_DCR rm: cannot remove `orion1im': No such file or directory rm: cannot remove `orion1im_weight': No such file or directory NORMAL: Starting imager::setdata NORMAL: Selecting data NORMAL: Performing selection on MeasurementSet NORMAL: Selecting on field and spectral window ids NORMAL: By selection 19516 rows are reduced to 4879 NORMAL: Finished imager::setdata 4.08 real 0.69 user 0.37 system NORMAL: successful readonly open of default-locked table oria1_DCR/NRAO_GBT_GLISH: 28 columns, 41 rows NORMAL: Defining image properties WARN: nx = 66 is not composite; nx = 72 or 64 will be more efficient WARN: ny = 66 is not composite; ny = 72 or 64 will be more efficient NORMAL: Setting processing options NORMAL: Starting imager::weight NORMAL: Weighting MS: IMAGING_WEIGHT column will be changed NORMAL: Natural weighting NORMAL: Sum of weights = 4879 NORMAL: Finished imager::weight 9.38 real 6.02 user 1.12 system NORMAL: Starting imager::makeimage NORMAL: Calculating image (without full skyequation) NORMAL: Making single dish image from corrected data NORMAL: Image is : orion1im NORMAL: Frequency = 8.4351 GHz, synthesized continuum bandwidth = 0.05 GHz NORMAL: Image polarization = Stokes I NORMAL: Frequency = 8.4351 GHz, synthesized continuum bandwidth = 0.05 GHz NORMAL: Image polarization = Stokes I NORMAL: Performing Single Dish gridding with convolution function SF Gridding will use specified common tangent point: 05:35:17.47 -05.23.07.50 J2000 NORMAL: Using prolate spheriodal wave function as the convolution function NORMAL: Transforming I only NORMAL: Finished imager::makeimage 13.06 real 2.06 user 0.52 system NORMAL: Starting imager::makeimage NORMAL: Calculating image (without full skyequation) NORMAL: Making single dish coverage function NORMAL: Image is : orion1im_weight NORMAL: Frequency = 8.4351 GHz, synthesized continuum bandwidth = 0.05 GHz NORMAL: Image polarization = Stokes I NORMAL: Frequency = 8.4351 GHz, synthesized continuum bandwidth = 0.05 GHz NORMAL: Image polarization = Stokes I NORMAL: Using prolate spheriodal wave function as the convolution function NORMAL: Transforming I only NORMAL: Finished imager::makeimage 2.12 real 1.54 user 0.15 system NORMAL: Successfully closed empty server: imager # Created 2 files: the uncorrected image: orion1im/ # & the coverage map: orion1im_weight/ # Create a corrected image by dividing the uncorrected image by the coverage map: gc.covercorr(); NORMAL: Starting server app_image NORMAL: Server started: /aips++/stable/linux/bin/app_image (AIPS++ version: 1.6 (build #363)) NORMAL: Selected bounding box : [1, 1, 1, 1] to [66, 66, 1, 1] (05:35:57.980, -05.33.12.210, I, 8.435100e+09Hz to 05:34:38.209, -05.13.20.955, I, 8.435100e+09Hz) NORMAL: Creating new statistics storage lattice of shape [9] NORMAL: NORMAL: Number points = 4.356000e+03 Sum = 4.355996e+03 NORMAL: Mean = 9.999990e-01 NORMAL: Variance = 1.923394e-01 Std dev = 4.385652e-01 NORMAL: Rms = 1.091922e+00 NORMAL: NORMAL: Minimum value 0.000000e+00 at [1, 1, 1, 1] (05:35:57.980, -05.33.12.210, I, 8.435100e+09Hz) NORMAL: Maximum value 1.426234e+00 at [7, 16, 1, 1] (05:35:50.610, -05.28.37.331, I, 8.435100e+09Hz) NORMAL: NORMAL: Creating image `orion1im_corr' of shape [66, 66, 1, 1] NORMAL: Created mask `mask0' # created a "corrected image" orion1im_corr/ # view the image: include 'image.g' im := image('orion1im_corr'); im.view() # works beautifully # image is what I expect # Try plotsource to fit a source: gc.plotsource('orion1im_corr',nsource=1); rm: cannot remove `tst.tbl': No such file or directory NORMAL: successful creation of default-locked table tst.tbl: 4 columns, 0 rows NORMAL: created table tst.tbl with zero rows NORMAL: Starting image::findsources NORMAL: Selected bounding box : [1, 1, 1, 1] to [66, 66, 1, 1] (05:35:57.980, -05.33.12.210, I, 8.435100e+09Hz to 05:34:38.209, -05.13.20.955, I, 8.435100e+09Hz) NORMAL: Found 1 sources SEVERE: Cannot convert units of brightness to Jy - will assume Jy NORMAL: Finished image::findsources 1.36 real 0.02 user 0.05 system NORMAL: Starting server componentlist NORMAL: Server started: /aips++/stable/linux/bin/componentlist (AIPS++ version: 1.6 (build #363)) NORMAL: defaultentryparser (dep) ready F # Only one source was found (that is all I asked for). OK # there is low level striping in the background in the direction of the otf scanning. # stats function now in the tools option. # max = 84.82 min = 17.29 rms = 22.16 (over total region) # matches Joe's # Trying starting with nsource=0: gc.plotsource('orion1im_corr',nsource=0); # dach, looks like I can't get rid of the "source 1" label NORMAL: successful creation of default-locked table tst.tbl: 4 columns, 0 rows NORMAL: created table tst.tbl with zero rows NORMAL: Starting image::findsources NORMAL: Selected bounding box : [1, 1, 1, 1] to [66, 66, 1, 1] (05:35:57.980, -05.33.12.210, I, 8.435100e+09Hz to 05:34:38.209, -05.13.20.955, I, 8.435100e+09Hz) LocalExec::SetStatus: abnormal child termination for /aips++/stable/linux/bin/app_image Server 'app_image' has failed unexpectedly! You will need to create the relevant tool again. If that causes unexpected behavior, please restart AIPS++ Please submit a bug-report using bug() if you can reproduce the problem. : no sources found File: ./gbtcalutils.g, Line 320 Stack: .() aips++ is hung, kill and restart. # do receptor 2 for image 1 glish -l gbtcalutils.g field_names(gc) gc.setdata('oria1_DCR'); gc.contcal(receptor=2); # NORMAL: Server started: /aips++/stable/linux/bin/imager (AIPS++ version: 1.7.00) # NOTE: this is a new version of AIPS++ # table oria1_DCR: 28 columns, 19516 rows gc.makeimage('orion1im2',receptor=2); # NORMAL: Natural weighting # NORMAL: Sum of weights = 4879 # NORMAL: Image is : orion1im2 # NORMAL: Frequency = 8.4351 GHz, synthesized continuum bandwidth = 0.05 GHz # NORMAL: Performing Single Dish gridding with convolution function SF # Gridding will use specified common tangent point: # 05:35:17.47 -05.23.07.50 J2000 # NORMAL: Using prolate spheriodal wave function as the convolution function # NORMAL: Making single dish coverage function # NORMAL: Image is : orion1im2_weight # NORMAL: Frequency = 8.4351 GHz, synthesized continuum bandwidth = 0.05 GHz # NORMAL: Using prolate spheriodal wave function as the convolution function gc.covercorr(); # NORMAL: Selected bounding box : # [1, 1, 1, 1] to [66, 66, 1, 1] # (05:35:57.980, -05.33.12.210, I, 8.435100e+09Hz to 05:34:38.209, -05.13.20.955, I, 8.435100e+09Hz) # NORMAL: Creating new statistics storage lattice of shape [9] # NORMAL: # NORMAL: Number points = 4.356000e+03 Sum = 4.355996e+03 # NORMAL: Mean = 9.999990e-01 # NORMAL: Variance = 1.923394e-01 Std dev = 4.385652e-01 # NORMAL: Rms = 1.091922e+00 # NORMAL: # NORMAL: Minimum value 0.000000e+00 at [1, 1, 1, 1] (05:35:57.980, -05.33.12.210, I, 8.435100e+09Hz) # NORMAL: Maximum value 1.426234e+00 at [7, 16, 1, 1] (05:35:50.610, -05.28.37.331, I, 8.435100e+09Hz) # NORMAL: # NORMAL: Creating image `orion1im2_corr' of shape [66, 66, 1, 1] # NORMAL: Created mask `mask0' include 'image.g' im := image('orion1im2_corr'); im.view() # max = 88.28 min= 18.67 rms = 18.87 (no emission region) rms = 23.64 (over total image) # matches Joe's # do image 2, receptors 1 & 2: gc.setdata('oria2_DCR'); gc.contcal(receptor=1); # table oria2_DCR: 28 columns, 19516 rows gc.makeimage('orion2im1',receptor=1); # NORMAL: Natural weighting # NORMAL: Sum of weights = 4879 # NORMAL: Image is : orion2im1 # NORMAL: Frequency = 8.4351 GHz, synthesized continuum bandwidth = 0.05 GHz # Gridding will use specified common tangent point: # 05:35:17.47 -05.23.07.50 J2000 # NORMAL: Using prolate spheriodal wave function as the convolution function # NORMAL: Making single dish coverage function # NORMAL: Image is : orion2im1_weight # NORMAL: Frequency = 8.4351 GHz, synthesized continuum bandwidth = 0.05 GHz gc.covercorr(); # NORMAL: Selected bounding box : # [1, 1, 1, 1] to [66, 66, 1, 1] # (05:35:57.980, -05.33.12.210, I, 8.435100e+09Hz to 05:34:38.209, -05.13.20.955, I, 8.435100e+09Hz) # NORMAL: Creating new statistics storage lattice of shape [9] # NORMAL: # NORMAL: Number points = 4.356000e+03 Sum = 4.355996e+03 # NORMAL: Mean = 9.999991e-01 # NORMAL: Variance = 1.924753e-01 Std dev = 4.387201e-01 # NORMAL: Rms = 1.091984e+00 # NORMAL: # NORMAL: Minimum value 0.000000e+00 at [1, 1, 1, 1] (05:35:57.980, -05.33.12.210, I, 8.435100e+09Hz) # NORMAL: Maximum value 1.447141e+00 at [8, 55, 1, 1] (05:35:49.372, -05.16.42.580, I, 8.435100e+09Hz) # NORMAL: # NORMAL: Creating image `orion2im1_corr' of shape [66, 66, 1, 1] include 'image.g' im := image('orion2im1_corr'); im.view() # striping is worse in this smaller image. Clearly modifies source shape. # worst stripe located at the northern edge of the main source # max = 83.89 min = 17.50 rms = 17.97 (no imission region) gc.contcal(receptor=2); # table oria2_DCR: 28 columns, 19516 rows gc.makeimage('orion2im2',receptor=2); # NORMAL: Natural weighting # NORMAL: Sum of weights = 4879 # NORMAL: Image is : orion2im2 # NORMAL: Frequency = 8.4351 GHz, synthesized continuum bandwidth = 0.05 GHz # Gridding will use specified common tangent point: # 05:35:17.47 -05.23.07.50 J2000 # NORMAL: Using prolate spheriodal wave function as the convolution function # NORMAL: Making single dish coverage function # NORMAL: Image is : orion2im2_weight # NORMAL: Frequency = 8.4351 GHz, synthesized continuum bandwidth = 0.05 GHz gc.covercorr(); # NORMAL: Selected bounding box : # [1, 1, 1, 1] to [66, 66, 1, 1] # (05:35:57.980, -05.33.12.210, I, 8.435100e+09Hz to 05:34:38.209, -05.13.20.955, I, 8.435100e+09Hz) # NORMAL: Creating new statistics storage lattice of shape [9] # NORMAL: # NORMAL: Number points = 4.356000e+03 Sum = 4.355996e+03 # NORMAL: Mean = 9.999991e-01 # NORMAL: Variance = 1.924753e-01 Std dev = 4.387201e-01 # NORMAL: Rms = 1.091984e+00 # NORMAL: # NORMAL: Minimum value 0.000000e+00 at [1, 1, 1, 1] (05:35:57.980, -05.33.12.210, I, 8.435100e+09Hz) # NORMAL: Maximum value 1.447141e+00 at [8, 55, 1, 1] (05:35:49.372, -05.16.42.580, I, 8.435100e+09Hz) # NORMAL: # NORMAL: Creating image `orion2im2_corr' of shape [66, 66, 1, 1] include 'image.g' im := image('orion2im2_corr'); im.view() # max = 88.21 min= 18.86 rms = 19.14 (no emission region) # same stripes as in pol1 (as expected) # do image 3, receptors 1 & 2: gc.setdata('oria3_DCR'); gc.contcal(receptor=1); # table oria3_DCR: 28 columns, 19516 rows gc.makeimage('orion3im1',receptor=1); # NORMAL: Natural weighting # NORMAL: Sum of weights = 4879 # NORMAL: Image is : orion3im1 # NORMAL: Frequency = 8.4351 GHz, synthesized continuum bandwidth = 0.05 GHz # Gridding will use specified common tangent point: # 05:35:17.47 -05.23.07.50 J2000 # NORMAL: Using prolate spheriodal wave function as the convolution function # NORMAL: Image is : orion3im1_weight # NORMAL: Frequency = 8.4351 GHz, synthesized continuum bandwidth = 0.05 GHz gc.covercorr(); # NORMAL: Selected bounding box : # [1, 1, 1, 1] to [66, 66, 1, 1] # (05:35:57.980, -05.33.12.210, I, 8.435100e+09Hz to 05:34:38.209, -05.13.20.955, I, 8.435100e+09Hz) # NORMAL: Creating new statistics storage lattice of shape [9] # NORMAL: # NORMAL: Number points = 4.356000e+03 Sum = 4.356001e+03 # NORMAL: Mean = 1.000000e+00 # NORMAL: Variance = 1.997561e-01 Std dev = 4.469409e-01 # NORMAL: Rms = 1.095313e+00 # NORMAL: # NORMAL: Minimum value 0.000000e+00 at [1, 1, 1, 1] (05:35:57.980, -05.33.12.210, I, 8.435100e+09Hz) # NORMAL: Maximum value 1.396221e+00 at [62, 47, 1, 1] (05:34:43.112, -05.19.09.188, I, 8.435100e+09Hz) # NORMAL: # NORMAL: Creating image `orion3im1_corr' of shape [66, 66, 1, 1] include 'image.g' im := image('orion3im1_corr'); im.view() # max = 82.87 min = 18.03 rms = 18.44 (no emission region) # major bad stripe right through the source center gc.contcal(receptor=2); # table oria3_DCR: 28 columns, 19516 rows gc.makeimage('orion3im2',receptor=2); # NORMAL: Natural weighting # NORMAL: Sum of weights = 4879 # NORMAL: Image is : orion3im2 # NORMAL: Frequency = 8.4351 GHz, synthesized continuum bandwidth = 0.05 GHz # Gridding will use specified common tangent point: # 05:35:17.47 -05.23.07.50 J2000 # NORMAL: Image is : orion3im2_weight # NORMAL: Frequency = 8.4351 GHz, synthesized continuum bandwidth = 0.05 GHz gc.covercorr(); # NORMAL: Selected bounding box : # [1, 1, 1, 1] to [66, 66, 1, 1] # (05:35:57.980, -05.33.12.210, I, 8.435100e+09Hz to 05:34:38.209, -05.13.20.955, I, 8.435100e+09Hz) # NORMAL: Creating new statistics storage lattice of shape [9] # NORMAL: # NORMAL: Number points = 4.356000e+03 Sum = 4.356001e+03 # NORMAL: Mean = 1.000000e+00 # NORMAL: Variance = 1.997561e-01 Std dev = 4.469409e-01 # NORMAL: Rms = 1.095313e+00 # NORMAL: # NORMAL: Minimum value 0.000000e+00 at [1, 1, 1, 1] (05:35:57.980, -05.33.12.210, I, 8.435100e+09Hz) # NORMAL: Maximum value 1.396221e+00 at [62, 47, 1, 1] (05:34:43.112, -05.19.09.188, I, 8.435100e+09Hz) # NORMAL: # NORMAL: Creating image `orion3im2_corr' of shape [66, 66, 1, 1] include 'image.g' im := image('orion3im2_corr'); im.view() # max 87.00 min = 19.39 rms = 19.84 (no emission region) # same rotten stripe right through the source... # do image 4, receptors 1 & 2: gc.setdata('oria4_DCR'); gc.contcal(receptor=1); # oria4_DCR: 28 columns, 144720 rows gc.makeimage('orion4im1',receptor=1); # NORMAL: Natural weighting # NORMAL: Sum of weights = 36180 MONSTER SUM HERE, MUCH HIGHER THAN OTHER MAPS. # NORMAL: Image is : orion4im1 # NORMAL: Frequency = 0 GHz, synthesized continuum bandwidth = 0.05 GHz # Gridding will use specified common tangent point: # 05:35:17.47 -05.23.07.50 J2000 # NORMAL: Image is : orion4im1_weight # NORMAL: Frequency = 0 GHz, synthesized continuum bandwidth = 0.05 GHz #### # NOTE: FREQUENCY = 0 GHz... gc.covercorr(); # NORMAL: Selected bounding box : # [1, 1, 1, 1] to [296, 328, 1, 1] # (05:38:19.373, -06.13.11.558, I, 0.000000e+00Hz to 05:32:17.291, -04.33.18.430, I, 0.000000e+00Hz) # NORMAL: Creating new statistics storage lattice of shape [9] # NORMAL: # NORMAL: Number points = 9.708800e+04 Sum = 9.708784e+04 # NORMAL: Mean = 9.999983e-01 # NORMAL: Variance = 4.999809e-02 Std dev = 2.236025e-01 # NORMAL: Rms = 1.024692e+00 # NORMAL: # NORMAL: Minimum value 0.000000e+00 at [1, 1, 1, 1] (05:38:19.373, -06.13.11.558, I, 0.000000e+00Hz) # NORMAL: Maximum value 2.533797e+00 at [8, 8, 1, 1] (05:38:10.757, -06.11.03.411, I, 0.000000e+00Hz) # NORMAL: # NORMAL: Creating image `orion4im1_corr' of shape [296, 328, 1, 1] include 'image.g' im := image('orion4im1_corr'); im.view() # THIS IS A HUGE IMAGE # max = 83.95 min = 17.15 rms = 17.43 (no emission region) # stripes obvious at a low level, significant real low-level emission # If I can get the stripes out, and fix frequency, # then this is the one for the poster. # created ps file 1.orion4pol1.ps # gc.contcal(receptor=2); # table oria4_DCR: 28 columns, 144720 rows gc.makeimage('orion4im2',receptor=2); # WARN: nx = 296 is not composite; nx = 300 will be more efficient # WARN: ny = 328 is not composite; ny = 360 or 324 will be more efficient # NORMAL: Natural weighting # NORMAL: Sum of weights = 36180 # NORMAL: Image is : orion4im2 # NORMAL: Frequency = 0 GHz, synthesized continuum bandwidth = 0.05 GHz # Gridding will use specified common tangent point: # 05:35:17.47 -05.23.07.50 J2000 # NORMAL: Using prolate spheriodal wave function as the convolution function # NORMAL: Image is : orion4im2_weight # NORMAL: Frequency = 0 GHz, synthesized continuum bandwidth = 0.05 GHz gc.covercorr(); # NORMAL: Selected bounding box : # [1, 1, 1, 1] to [296, 328, 1, 1] # (05:38:19.373, -06.13.11.558, I, 0.000000e+00Hz to 05:32:17.291, -04.33.18.430, I, 0.000000e+00Hz) # NORMAL: Creating new statistics storage lattice of shape [9] # NORMAL: # NORMAL: Number points = 9.708800e+04 Sum = 9.708784e+04 # NORMAL: Mean = 9.999983e-01 # NORMAL: Variance = 4.999809e-02 Std dev = 2.236025e-01 # NORMAL: Rms = 1.024692e+00 # NORMAL: # NORMAL: Minimum value 0.000000e+00 at [1, 1, 1, 1] (05:38:19.373, -06.13.11.558, I, 0.000000e+00Hz) # NORMAL: Maximum value 2.533797e+00 at [8, 8, 1, 1] (05:38:10.757, -06.11.03.411, I, 0.000000e+00Hz) # NORMAL: # NORMAL: Creating image `orion4im2_corr' of shape [296, 328, 1, 1] # include 'image.g' im := image('orion4im2_corr'); im.view() # max = 88.12, min = 18.45 rms = 18.75 (no emission region) im.summary() Need to resolve stripe problem, frequency in map 4, calibration in large image. In the mean time, try to average some of the images. include 'image.g' im1 := imagecalc(outfile='orion1ave', pixels='(orion1im_corr + orion1im2_corr)/2'); # NORMAL: Creating image `orion1ave' of shape [66, 66, 1, 1] im := image('orion1ave'); im.summary() Axis Coord Type Name Proj Shape Tile Coord value at pixel Coord incr Units ------------------------------------------------------------------------------------------------ 1 1 Direction Right Ascension SIN 66 66 05:35:17.470 34.00 -1.832705e+01 arcsec 2 1 Direction Declination SIN 66 66 -05.23.07.500 34.00 1.832705e+01 arcsec 3 2 Stokes Stokes 1 1 I 4 3 Spectral Frequency 1 1 8.435100e+09 1.00 5.000000e+07 Hz Velocity -inf 1.00 nan km/s Why is the velocity infinite? im.view() # max = 86.82 min = 17.98 rms = 18.20 (no emission region) # looks like what I expect, stripes are still there at low level include 'image.g' im4 := imagecalc(outfile='orion4ave', pixels='(orion4im1_corr + orion4im2_corr)/2'); # NORMAL: Creating image `orion4ave' of shape [296, 328, 1, 1] im := image('orion4ave'); im.summary() Axis Coord Type Name Proj Shape Tile Coord value at pixel Coord incr Units ------------------------------------------------------------------------------------------------ 1 1 Direction Right Ascension SIN 296 296 05:35:17.470 149.00 -1.832705e+01 arcsec 2 1 Direction Declination SIN 328 328 -05.23.07.500 165.00 1.832705e+01 arcsec 3 2 Stokes Stokes 1 1 I 4 3 Spectral Frequency 1 1 0.000000e+00 1.00 5.000000e+07 Hz Velocity 2.997925e+05 1.00 0.000000e+00 km/s Humm, velocity looks right but the frequency is zero. im.view() # max = 86.03 min = 17.80 rms = 18.09 (no emission region) # looks like what I expect try averaging the big image 4 and small image 1 to see if the software can figure out how to grid them together: im5 := imagecalc(outfile='orion1.4ave', pixels='(orion1ave + orion4ave)/2'); SEVERE: Caught an exception! Event type=create exception=LatticeExprNode - coordinates of operands mismatch Scanned so far: (orion1ave + orion4ave) : Caught an exception! Event type=create exception=LatticeExprNode - coordinates of operands mismatch Scanned so far: (orion1ave + orion4ave) File: servers.g, Line 1134 Stack: .(), image.g line 2413 imagecalc() Nope, software can't figure it out. This should be in the documentation. (submitted a documentation enhancement request) try something else... im := image('orion4ave'); im.shape(); [296 328 1 1] im := image('orion1ave'); im.shape(); [66 66 1 1] im.coordsys(); -- this doesn't work...may not be meant to work. Lets try a straight regrid of orion4ave to the orion1ave coordinate frame & size: im01 := image('orion1ave'); im04 := image('orion4ave'); cs01 := im01.coordsys(); im05 := im04.regrid(outfile='orion4regrid', csys=cs01, shape=im01.shape()); # NORMAL: Starting image::regrid # NORMAL: Selected bounding box : # [1, 1, 1, 1] to [296, 328, 1, 1] # (05:38:19.373, -06.13.11.558, I, 0.000000e+00Hz to 05:32:17.291, -04.33.18.430, I, 0.000000e+00Hz) # NORMAL: Creating image 'orion4regrid' of shape [66, 66, 1, 1] # NORMAL: Created and initialized mask `mask0' # WARN: The Stokes axis cannot be regridded - removing from list # NORMAL: Regridding output axes [1, 2] which are of Coordinate type Direction # NORMAL: Regridding output axis 4 which is of Coordinate type Spectral # NORMAL: Finished image::regrid # E-gads! it seems to have worked... im := image('orion4regrid'); im.view() # it looks like the right size, now try to average the two. im06 := imagecalc(outfile='orion1.4ave', pixels='(orion1ave + orion4regrid)/2'); # NORMAL: Creating image `orion1.4ave' of shape [66, 66, 1, 1] im := image('orion1.4ave'); im.view() # max = 86.43 min = 18.01 rms 18.17 # max and min are about right but this is now the average of 4 images and the # rms has not decreased significantly. im.summary() Image name : orion1.4ave Axis Coord Type Name Proj Shape Tile Coord value at pixel Coord incr Units ------------------------------------------------------------------------------------------------ 1 1 Direction Right Ascension SIN 66 66 05:35:17.470 34.00 -1.832705e+01 arcsec 2 1 Direction Declination SIN 66 66 -05.23.07.500 34.00 1.832705e+01 arcsec 3 2 Stokes Stokes 1 1 I 4 3 Spectral Frequency 1 1 8.435100e+09 1.00 5.000000e+07 Hz Velocity -inf 1.00 nan km/s frequency of the averaged image now matches images 1, 2, & 3... Compare with input images into this average... im10 := image('orion1im2_corr') im10.view() im11 := image('orion1ave') im11.view() im12 := image('orion4regrid') im12.view() the stripes have smoothed out after combining orion1 and orion4 images. Good. But, they remain high. It appears that because the stripes are not random, I can't decrease the noise when I average the images... ===================================================================== ===================================================================== ===================================================================== Ron says that I have to do spatial smoothing on the big map to remove the stripes. Then I have to do the FFT trick on the smaller maps. uv-based cutting removes selected spatial scales ==> it may mess up the GBT-VLA combo. Compare the big and small maps carefully and check the final combo results very closely. GBT Baseline documentation refers only to SSA, which is no longer valid. ------------------------------------------------------------------------ try out new spatial baseline capability in gbtcalutils.g: baseline region on either side of src hardwired into gbtcalutils.g script. Change script if I want to change region. source /aips++/stable/aipsinit.csh /export/data_1/aips/dss/app/orion.gb AIPS++ Version 1.7 Build 000 glish -l gbtcalutils.g # do image 4, receptors 1 & 2: gc.setdata('oria4_DCR'); # do the basic calibration with "receptor 1" data to get T_ANT # from measured voltages with the noise-diode ON and OFF. gc.contcal(receptor=1,baselinerm=T); # oria4_DCR: 28 columns, 144720 rows # lots of baseline messages # Select, regrid, and construct an "unnormalized" image with SF grid fctn: gc.makeimage('orion4im1',receptor=1); # WARN: nx = 296 is not composite; nx = 300 will be more efficient # WARN: ny = 328 is not composite; ny = 360 or 324 will be more efficient # NORMAL: Natural weighting # Gridding will use specified common tangent point: # 05:35:17.47 -05.23.07.50 J2000 # Create a corrected image by dividing the uncorrected image by the coverage map: gc.covercorr(); # orion4im1_corr created of shape [296, 328, 1, 1] include 'image.g' im1 := image('orion4im1_corr'); im1.view() # max = 66.37 min = -1.69 rms = 7e-3 im1.summary() # app crash and dump the core when im.view is active and # gc.contcal is run below... (tell Joe, no bug report -- still too new) # app crash when trying to view 2 images at once... # do the basic calibration with "receptor 1" data to get T_ANT # from measured voltages with the noise-diode ON and OFF. gc.contcal(receptor=2,baselinerm=T); # table oria4_DCR: 28 columns, 144720 rows # Select, regrid, and construct an "unnormalized" image with SF grid fctn: gc.makeimage('orion4im2',receptor=2); # WARN: nx = 296 is not composite; nx = 300 will be more efficient # WARN: ny = 328 is not composite; ny = 360 or 324 will be more efficient # NORMAL: Natural weighting # Gridding will use specified common tangent point: # 05:35:17.47 -05.23.07.50 J2000 # Create a corrected image by dividing the uncorrected image by the coverage map: gc.covercorr(); # NORMAL: Creating image `orion4im2_corr' of shape [296, 328, 1, 1] include 'image.g' im2 := image('orion4im2_corr'); im2.view() # max = 69.2, min = -0.18 rms = 7e-3 (no emission region) im2.summary() # shut down aips++ for the day. # # restart: glish -l gbtcalutils.g # # roughly calibrate each image to Jy scale: # # multiply each polarization by Tcal values from Ron to get Ta: include 'image.g' im1 := imagecalc(outfile='orion4im1_ta', pixels='(orion4im1_corr) * 3.17'); im2 := imagecalc(outfile='orion4im2_ta', pixels='(orion4im2_corr) * 3.03'); # # not sure how to calibrate out atmospheric attenuation # or aperature efficiency as a function of elevation (row)... # not a big effect (< 2%), skip for now. Ask Joe. # # Scale to Jy: 1.76 +/- 0.11 K/Jy so divide the maps by 1.76 to get Jy: im3 := imagecalc(outfile='orion4im1_jy', pixels='(orion4im1_ta) / 1.76'); im4 := imagecalc(outfile='orion4im2_jy', pixels='(orion4im2_ta) / 1.76'); im3.view() # # Flux doesn't show, only sum... assuming sum = flux density? # max = 119.5 Jy/bm min = -0.30 Jy/bm RMS = 1.7e-2 Jy/bm = 17 mJy/bm im4.view() max = 119.2 Jy/bm min = -0.31 Jy/bm RMS = 1.5e-2 Jy/bm = 15 mJy/bm # # average the 2 together: im5 := imagecalc(outfile='orion4ave', pixels='(orion4im1_jy + orion4im2_jy)/2'); # im5 := image('orion4ave'); im5.summary() Axis Coord Type Name Proj Shape Tile Coord value at pixel Coord incr Units ------------------------------------------------------------------------------------------------ 1 1 Direction Right Ascension SIN 296 296 05:35:17.470 149.00 -1.832705e+01 arcsec 2 1 Direction Declination SIN 328 328 -05.23.07.500 165.00 1.832705e+01 arcsec 3 2 Stokes Stokes 1 1 I 4 3 Spectral Frequency 1 1 0.000000e+00 1.00 5.000000e+07 Hz Velocity 2.997925e+05 1.00 0.000000e+00 km/s Humm, velocity looks right but the frequency is zero. im5.view() xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx # convert map to Jy/bm for total flux measurement: rb5 := im5.restoringbeam(); im5.restoringbeam(); # [=] # GBT beam at 3.55 cm 8.435 GHz: 1.46' = 87.6" # GBT beam at 3.55 cm (8.435 GHz) = 1.46' = 4.25e-4 radians, # but, since the GBT data are convolved onto the image grid using an # antenna primary beam, the effective beam is sqrt(2) bigger. # so set GBT beam = 2.06' = 124" rb5.major.value := 124. rb5.major.unit := 'arcsec' rb5.minor.value := 124. rb5.minor.unit := 'arcsec' rb5.positionangle.value := 0 rb5.positionangle.unit := 'deg' im5.setrestoringbeam(beam=rb5); # NORMAL: Set restoring beam # Major : 124 arcsec # Minor : 124 arcsec # Position Angle : 0 deg im5.restoringbeam(); im5.summary(); im5.setbrightnessunit('Jy/beam'); im5.brightnessunit(); # Jy/beam im5.view(); # max = 119.4 Jy/bm min = -0.31 Jy/bm # region 1: RMS = 1.3e-2 Jy/bm = 13 mJy/bm \ # region 2: RMS = 9.1e-3 Jy/bm = 9 mJy/bm | average = 10 mJy/bm # region 3: RMS = 6.8e-3 Jy/bm = 7 mJy/bm / # decreased rms -- good # dynamic range: S(max)/N = 119.4/1e-2 = 11,940 -- monster good. # create a nice image called 1.orion4ave.ps # actual GBT beam: total flux = 861.7 Jy (flux in big orion part = 854.1 Jy) # GBt beam sqrt(2) bigger: flux = 430.0 Jy (flux in big orion part = 424.5 Jy) # compare with VLA: max = 0.89 Jy/bm, flux = 189 Jy xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx im5.done(); include 'image.g'; imgbt := image('orion4ave'); field_names(imgbt); # busy result *agent* adddegaxes ada addnoise arrconvolve boundingbox # bb brightnessunit bu calc calcmask close convolve2d c2d coordmeasures # coordsys deconvolvecomponentlist dcl delete done fft findsources fs # fitsky fitprofile getchunk getregion hanning haslock histograms histo # history id insert isopen ispersistent lock maskhandler mh # maskhandlergui mhgui maxfit miscinfo modify moments momentsgui name # open pixelvalue putchunk putregion regrid rename replacemaskedpixels # rmp restoringbeam rb sepconvolve sc sepconvolvegui scgui set # setbrightnessunit sbu setcoordsys sethistory setmiscinfo # setrestoringbeam srb shape statistics stats subimage subim summary # toascii tofits topixel toworld unlock view convertsm type imgbt.shape(); # [296 328 1 1] imgbt.summary(); ################################################################ ################################################################ # The app vla image is called vla.orion.mem7/ # combine the 2 images... include 'table.g' tabledelete('orion4rgrd') tabledelete('orion4rgrd.scaled') include 'image.g'; imgbt := image('orion4ave'); csgbt := imgbt.coordsys(); csgbt.settelescope('GBT'); imgbt.setcoordsys(csgbt); imgbt.summary(); # freq, beam, & units unknown, telescope correct imvla := image('vla.orion.mem7'); imvla.summary() csvla := imvla.coordsys(); csvla.settelescope('VLA'); imvla.setcoordsys(csvla); imvla.summary(); imregrid := imgbt.regrid(outfile='orion4rgrd', csys=csvla, shape=imvla.shape(), axes=[1,2]); # NORMAL: Selected bounding box : # [1, 1, 1, 1] to [296, 328, 1, 1] # (05:38:19.373, -06.13.11.558, I, 0.000000e+00Hz to 05:32:17.291, -04.33.18.430, I, 0.000000e+00Hz) # NORMAL: Creating image 'orion4rgrd' of shape [300, 300, 1, 1] imregrid.done(); imgbt.done(); imvla.done(); imrgrd := image('orion4rgrd'); imrgrd.summary(); # now says the GBT image has a VLA telescope... ==> bug() # frequency is correct now ==> bug() freq should not have changed. . xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx x csrgrd := imrgrd.coordsys(); x x csrgrd.settelescope('GBT'); x x imrgrd.setcoordsys(csrgrd); x -- this cannot be done or image is destroyed. x imrgrd.summary(); x x # now we have a GBT telescope back x xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx rbrgrd := imrgrd.restoringbeam(); imrgrd.restoringbeam(); rbrgrd.major.value := 124. rbrgrd.major.unit := 'arcsec' rbrgrd.minor.value := 124. rbrgrd.minor.unit := 'arcsec' rbrgrd.positionangle.value := 0 rbrgrd.positionangle.unit := 'deg' imrgrd.setrestoringbeam(beam=rbrgrd); # NORMAL: Set restoring beam # Major : 124 arcsec # Minor : 124 arcsec # Position Angle : 0 deg imrgrd.setbrightnessunit('Jy/beam'); imrgrd.brightnessunit(); # Jy/beam imrgrd.view(); # max = 119.3 Jy/bm, sum = 1.536e6 - no beam correction # max = 119.3 Jy/bm, flux = 352.7 Jy - 124" restoring beam # compare with VLA: max = 0.89 Jy/bm, flux = 189 Jy imrgrd.done(); ################################################################ # # Scale the GBT image by the ratio of # beam solid angles # # First, find the VLA beam solid angle # include 'image.g'; imvla := image('vla.orion.mem7'); beam := imvla.restoringbeam(); beam # [major=[value=8.3126812, unit=arcsec], # minor=[value=7.41258144, unit=arcsec], # positionangle=[value=-32.9935684, unit=deg]]include 'quanta.g'; vlasa := dq.canonical(dq.mul(dq.mul(beam.major, beam.minor), 1.1331)); vlasa # [value=1.64107555e-09, unit=rad2] # # GBT beam at 3.55 cm (8.435 GHz) = 1.46' = 4.25e-4 radians, # but, since the GBT data are convolved onto the image grid using an # antenna primary beam, the effective beam is sqrt(2) bigger. # so set GBT beam = 2.06' = 5.99e-4 radians # scaling estimate: # Tim: This means that the solid angle is # off by a factor of two and so according to the original script that I # sent an extra factor of two increase in the GBT data must be multiplied # in. # OR, equivalently, the beam size can be set to sqrt(2) higher (roughly 2arcmin). gbtbeam := dq.quantity(sqrt(2)*0.000425, 'rad'); gbtbeam # [value=0.000601040764, unit=rad] gbtsa := dq.canon(dq.mul(gbtbeam, gbtbeam)); gbtsa # [value=3.6125e-07, unit=rad2] # gbtsa/vlasa = ratio: ratio := dq.div(gbtsa, vlasa).value; ratio # gbtsa/vlasa = ratio: 220.130023 # # Now use the image calculator to rescale the GBT image # imscale := imagecalc('orion4rgrd.scaled', spaste(ratio, '*', 'orion4rgrd')); imscale.done(); imvla.done(); imsc := image('orion4rgrd.scaled'); imsc.summary(); # header telescope=VLA and freq OK, no beam xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx # change to Jy/bm before the image feathering ... # setrestoringbeam to be the vla beam imvla := image('vla.orion.mem7'); imvla.summary(); # Restoring Beam : 8.31268 arcsec, 7.41258 arcsec, -32.9936 deg vlabeam := imvla.restoringbeam(); imvla.restoringbeam(); # [major=[value=8.3126812, unit=arcsec], minor=[value=7.41258144, unit=arcsec], # positionangle=[value=-32.9935684, unit=deg]] rb := imsc.restoringbeam(); imsc.restoringbeam(); # [=] imsc.setrestoringbeam(beam=vlabeam); NORMAL: Set restoring beam # Major : 8.31268 arcsec # Minor : 7.41258 arcsec # Position Angle : -32.9936 deg imsc.restoringbeam(); imsc.summary(); xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx imsc.setbrightnessunit('Jy/beam'); imsc.brightnessunit(); # Jy/beam imsc.view(); # scaled image stats: max = 2.6e4 sum = 3.4e8 -- no beams, no nothing # scaled image stats: max = 2.6e4 Jy/bm, flux = 2e7 Jy -- vlabeam CRAP # compare with GBT: # max = 119.3 Jy/bm, flux = 352.7 Jy 2.06' beam # # max = 119.3 Jy/bm, flux = 707 Jy 1.46' beam # compare with VLA: max = 0.89 Jy/bm, flux = 189 Jy imsc.done(); # # ################################################################ # Use imager to do the feathering # # first: tell imager to use the voltage pattern (primary beam) of the GBT. # If you don't, the image will be horribly confused. This is because we must # account for the voltage pattern when imaging the same location on the image but # from different pointings. Get the voltage pattern from the original # GBT data set: # include 'imager.g'; im := imager('oria4_DCR'); im.setvp(dovp=T); # # feather info from : # Synopsis: feather(output_image, highres, lowres, usedefaultvp, vptable, async) # Description: # Basically the "imerg" algorithm of AIPS and SDE, or the "feather" # algorithm of MIRIAD, we regrid the total power (or low resolution) # image onto the interferometer (or high resolution) image, Fourier # transform both the interferometer and single dish images, down weight # the Fourier transform of the interferometer image by 1.0 - FT(low res # beam), add the weighted interferometer Fourier plane to the single # dish Fourier plane, and transform back into the image plane. # ##### The scaled orion image MUST have the telescope set to VLA ##### No restoring beams can be set in the scaled-GBT image # im.feather('orion.gbt+vla.image', 'vla.orion.mem7', 'orion4rgrd.scaled'); # NORMAL: Starting imager::feather # NORMAL: feathering together high and low resolution images # NORMAL: Regridding output axes [1, 2] which are of Coordinate type Direction # NORMAL: Shutting down server : /aips++/stable/linux/bin/imager im.done(); # # Look at the result # include 'image.g'; imcombo := image('orion.gbt+vla.image'); imcombo.view() imcombo.summary() # telescope is now VLA ... # restoring beam is gone... imvla := image('vla.orion.mem7'); imvla.summary(); # Restoring Beam : 8.31268 arcsec, 7.41258 arcsec, -32.9936 deg vlabeam := imvla.restoringbeam(); imvla.restoringbeam(); # [major=[value=8.3126812, unit=arcsec], minor=[value=7.41258144, unit=arcsec], # positionangle=[value=-32.9935684, unit=deg]] rbcombo := imcombo.restoringbeam(); imcombo.restoringbeam(); # [=] vlabeam.major.value := sqrt(2)*8.31268 vlabeam.minor.value := sqrt(2)*7.41258144 imvla.setrestoringbeam(beam=vlabeam); imcombo.setrestoringbeam(beam=vlabeam); # NORMAL: Set restoring beam # Major : 11.7559 arcsec --- Not sure which beam to use... # Minor : 10.483 arcsec # Position Angle : -32.9936 deg NORMAL: Set restoring beam # Major : 8.31268 arcsec # Minor : 7.41258 arcsec # Position Angle : -32.9936 deg imcombo.restoringbeam(); imcombo.setbrightnessunit('Jy/beam'); imcombo.brightnessunit(); # Jy/beam imcombo.summary(); imcombo.view(); # combo: max = 1.82 Jy/bm, flux = 832.3 Jy -too high # compare with GBT: max = 119 Jy/bm, flux = Jy # compare with VLA: max = 0.89 Jy/bm, flux = 189 Jy # display settings saved in 'displaydata defaults' -- it can be restored imcombo.done(); include 'image.g'; imvla := image('vla.orion.mem7'); imvla.view() imvla.done() created ps files of gbt, vla, & gbt+vla images: 1.gbt+vla.ps 1.only-gbt.ps 1.only-vla.ps