Difference between revisions of "Einstein-Face (CASA 3.2)"

From CASA Guides
Jump to navigationJump to search
Line 74: Line 74:
 
'''Step 4''' Add FITS header keywords and change the format
 
'''Step 4''' Add FITS header keywords and change the format
  
At this stage, we need to perform some manipulations on the FITS file to get it readable by simdata (an 8bit to 16bit conversion) and trim it down to 300x300 pixels.  
+
At this stage, we need to perform some manipulations on the FITS file to get it readable by simdata (an 8bit to 16bit conversion) and trim it down to 300x300 pixels. For convenience, we are also adding a World Coordinate System (WCS) at this point (though this is not strictly necessary as the WCS parameters can be supplied directly to the simdata task). As with most operations in CASA, we may either use high level tasks, or the low level toolkit to perform these operations. Users who are more familiar with IDL will probably prefer the toolkit, as it provides similar functionality to the IDL astronomy library. Users who are more familiar with AIPS or IRAF, or who are novices, will probably prefer the tasks approach. For completeness, we give both here.
For convenience, we are also adding a WCS at this point. We need to use some of the low level toolkit
 
tasks to perform these conversions:
 
  
Read the FITS file into CASA, trim it, and write it out as a 16-bit integer
+
Toolkit version:
 +
 
 +
First, read in the FITS file. Then define region "box" to trim. The subimage tool is then used to create "testimage2" which has the box applied, and the object becomes im2. im2 then has WCS added and is written to a 16bit integer fits file.
  
 
<source lang="python">
 
<source lang="python">
Line 84: Line 84:
 
box = rg.box([0,0],[299,299])
 
box = rg.box([0,0],[299,299])
 
im2 = ia.subimage('testimage2',box,overwrite=T)
 
im2 = ia.subimage('testimage2',box,overwrite=T)
 +
csys = im2.coordsys()
 +
csys.setdirection(refcode='EQUATORIAL',proj='SIN',projpar=[0,0],refpix=[150,150], refval="52.5deg -28.5deg", incr="-0.043arcsec,0.043arcsec,1.0MHz")
 +
ep = 2000.0
 +
csys.setepoch(ep)
 
ok = im2.tofits('einstein16.fits',bitpix=16,overwrite=true)
 
ok = im2.tofits('einstein16.fits',bitpix=16,overwrite=true)
 +
im2.done()
 
ia.close()
 
ia.close()
im2.close()
+
 
 
</source>
 
</source>
  

Revision as of 13:09, 7 March 2011

Simulations using non-science images: the face of Einstein

Simdata can be used to simulate any digitized image. These toy models can be particularly useful for examining the effects of varying uv-coverage on image fidelity if the "truth" model is a familiar object or image. In this example (which is on page 13 of the ALMA Early Science Primer)we use the face of Albert Einstein.

Step 1: obtain your image. Typically from the internet.

Einstein.jpg

In this case, it is a jpg file.

Step 2: Convert your image to FITS

Various software programs have conversion to FITS enabled. The (GIMP) was used in this case. A handy list of FITS conversion programs is maintained by GSFC here

For the GIMP, start up the software

>gimp &

and in the main window select "Open" from the "File" menu.

The image will open up in a new window, you can use the GIMP to modify the image (adjust contrast, colormap etc).

Then, select "Save as" from the "File" menu in the window containing the image, and hit "Select File Type" in the dialog box to bring up the file type options, and select "Flexible Image Transport System". Pick a name for your file ending in .fits, e.g. einstein.fits

Gimp save.jpg

Step 3: check your file

Read your file into e.g. ds9to check that a valid FITS file has been produced. You can also examine the pixel values to examine the scaling of the image. In this case, Einstein's forehead has pixel values around 230 and the background around 40, so there is plenty of contrast. You can also examine the image header in ds9. Under the "File" menu, select "Display FITS Header" and examine the output. Make sure that SIMPLE = T, NAXIS=2 and check BITPIX. In this case, BITPIX=8, which is not valid for reading into CASA, so we need to change that at the next step.

FITS header produced by GIMP:

SIMPLE = T

BITPIX = 8

NAXIS = 2

NAXIS1 = 300

NAXIS2 = 327

BZERO = 0.000000

BSCALE = 1.000000

DATAMIN = 0.000000

DATAMAX = 255.000000

HISTORY THIS FITS FILE WAS GENERATED BY GIMP USING FITSRW

COMMENT FitsRW is (C) Peter Kirchgessner (peter@kirchgessner.net), but available

COMMENT under the GNU general public licence.

COMMENT For sources see http://www.kirchgessner.net

COMMENT Image type within GIMP: GIMP_GRAY_IMAGE

END

Step 4 Add FITS header keywords and change the format

At this stage, we need to perform some manipulations on the FITS file to get it readable by simdata (an 8bit to 16bit conversion) and trim it down to 300x300 pixels. For convenience, we are also adding a World Coordinate System (WCS) at this point (though this is not strictly necessary as the WCS parameters can be supplied directly to the simdata task). As with most operations in CASA, we may either use high level tasks, or the low level toolkit to perform these operations. Users who are more familiar with IDL will probably prefer the toolkit, as it provides similar functionality to the IDL astronomy library. Users who are more familiar with AIPS or IRAF, or who are novices, will probably prefer the tasks approach. For completeness, we give both here.

Toolkit version:

First, read in the FITS file. Then define region "box" to trim. The subimage tool is then used to create "testimage2" which has the box applied, and the object becomes im2. im2 then has WCS added and is written to a 16bit integer fits file.

ia.fromfits(outfile='testimage',infile='einstein.fits',overwrite=T)
box = rg.box([0,0],[299,299])
im2 = ia.subimage('testimage2',box,overwrite=T)
csys = im2.coordsys()
csys.setdirection(refcode='EQUATORIAL',proj='SIN',projpar=[0,0],refpix=[150,150], refval="52.5deg -28.5deg", incr="-0.043arcsec,0.043arcsec,1.0MHz")
ep = 2000.0
csys.setepoch(ep)
ok = im2.tofits('einstein16.fits',bitpix=16,overwrite=true)
im2.done()
ia.close()

Below is the IDL version. This routine is written in IDL, using the IDL astronomy library, but similar manipulations can be carried out in IRAF, or using the python PyWCS and PyFITS libraries, available from the astropython project.

The IDL script is in File:Make 2dimage.pro.txt (remove the .txt from the filename before using).

IDL>make_2dimage,'einstein.fits',0,299,27,326

The IDL code performs the following manipulations:

1) Reads in the FITS file as a 2D array, trims it to 300x300 pixels and converts it to real, 300x300x1 array (the third dimension is added for generality to allow the construction of an image cube, it is not actually necessary in this particular case).

2) Creates header keywords corresponding to the axis types (CTYPE1,2,3) values at the reference pixels (CRVAL1,2,3), the reference pixel positions (CRPIX1,2,3) and the axis increments (CDELT1,2,3), and the epoch (EPOCH).

3) Writes out the modified FITS file as "twodmodel.fits"

If you want to skip the above steps, the fits file is File:Twodmodel.fits.txt. download it and copy it to twodmodel.fits

Step 5 Start CASA and prepare inputs for simdata

Start with the 10min full science observation. Inputs to simdata are given below. The integration time is set much longer than realistic (300s, compared to 1-10s in practice) to speed the computation. The map spacing is set to ensure that only one pointing is observed:

>casapy

default 'simdata'
project = 'fs_cfg8_10m'
modifymodel = F
skymodel = 'twodmodel.fits'
setpointings = T
integration = '300s'
mapsize = ['1arcmin','1arcmin']
maptype = 'hexagonal'
pointingspacing = '1arcmin'
predict = T

Antenna configuration: ALMA antenna configuration files are stored in a directory that depends on your CASA installation. To be sure of finding them, identify the CASAPATH variable using the os.getenv command, and pick the configuration you want. Details on configuration choices are given in the M51 simulation guide [1].


repodir=os.getenv("CASAPATH").split(' ')[0]
antennalist        =  repodir+"/data/alma/simmos/alma.out08.cfg"
totaltime = '600s'
thermalnoise = ""
image = T
vis = '$project.ms'
imsize = [300,300]
cell = '0.043arcsec'
niter = 2000
weighting = 'natural'
analyze=F
simdata

The output image should have a synthesized beam of 0.62"x0.56" and look something like: Einstein fs cfg8 10min.gif

Now we repeat for an 1hr observation:

tget simdata
project = 'fs_cfg8_1hr'
totaltime = '3600s'
simdata

Which should look something like: Einstein fs cfg8 1hr.gif

Finally, two Early Science simulations, using the 250m configuration. One 10min simulation:

tget simdata
project = 'es_cfg250_10m'
antennalist = repodir+"/data/alma/simmos/alma.early.250m.cfg"
totaltime = '600s'
simdata

which looks like this: Einstein es cfg250 10min.gif

and a 4hr simulation:

tget simdata
project = 'es_cfg250_4hr'
totaltime = '14400s'
simdata

which looks like this: Einstein es cfg250 4hr.gif

Further experiments:

Some more things to try:

An 8hr observation shows the improvement obtained by obtaining fuller uv-coverage in the full science array:

tget simdata
antennalist = repodir+"/data/alma/simmos/alma.out08.cfg"
project = 'fs_cfg8_8hr'
totaltime = '28800s'
simdata

Which should look something like: Einstein fs cfg8 8hr.gif

An attempt to make a higher resolution image shows what happens when short spacings are missing in the configuration. Configuration 16 has a 0.17x0.15 beam, still better than Nyquist sampling of the model image (which has 0.043" pixels). However, the lack of short spacings in the configuration leads to poorly sampled structure on large spatial scales. In practice, one would need to combine these observations with a set in a more compact configuration (such as 8) to sample both the large and small spatial structures.

tget simdata
antennalist = repodir+"/data/alma/simmos/alma.out16.cfg"
project = 'fs_cfg16_1hr'
totaltime = '3600s'
simdata

The result is: Einstein fs cfg16 1hr.gif