This site requires javascript to be enabled to fuction correctly


How to regrid a CASA image to match another one?
13 April 2011 12:58 PM
Here is a little script "demo_imageregrid.py" to demonstrate how to make the grids of two images agree.

The script first creates its own test data by generating two images with different coordinate systems and angular grids but with a pattern of the same angular size. This part can be ignored by most users but it demonstrates how to modify coordinate systems. The script then regrids one of the images to the grid of the other. It finally tests if a calculation of the difference of the images with the immath task works, first on the original images (it doesn't), and then on the regridded one (there it does).

Look at the individual output images smaller.im, bigger.im, bigger_regridded.im, and difference.im with the imview task in CASA to see what happens.

# Demo of image regridding (as a reply to helpdesk ticket 1530)
# D. Petry, ESO, 14 April 2011
# Requires CASA 3.1, verified to work up to CASA 4.3.0

import numpy as np
# prepare test data
a = np.array([1.]) # fill image with value 1.0
b = np.array([3.]) # fill image with value 3.0
abigger = np.resize(a,(20,20)) # 20 by 20 pixels
bsmaller = np.resize(b,(10,10)) # 10 by 10 pixels
# set a little square in the bigger image to a different value
for i in range(7,13):
    for j in range(7,13):
        abigger[i][j] = 10.
# set a corresponding square in the smaller image to the same value
for i in range(3,7):
    for j in range(3,7):
        bsmaller[i][j] = 10.


os.system('rm -rf bigger.im smaller.im bigger_regridded.im difference.im')

ia.fromarray(outfile='bigger.im',pixels = abigger, overwrite=True)
cs=ia.coordsys()
cs.setepoch({'type': 'epoch', 'm0': {'value': 0.0, 'unit': 'd'}, 'refer': 'LAST'})
cs.setprojection(parameters = [ 0.,  0.], type= 'SIN')
cs.setreferencecode('J2000')
cs.setunits(['deg','deg'])
cs.setreferencepixel(type='direction', value=[9.,9.])
cs.setreferencevalue(type='direction', value=[0.,0.])
cs.setincrement(type='direction', value=[0.001,0.001])
cs.setnames(type='direction', value=['RA','Dec'])
ia.setcoordsys(cs.torecord())
ia.close()

ia.fromarray(outfile='smaller.im',pixels = bsmaller, overwrite=True)
cs=ia.coordsys()
cs.setepoch({'type': 'epoch', 'm0': {'value': 0.0, 'unit': 'd'}, 'refer': 'LAST'})
cs.setprojection(parameters = [ 0.,  0.], type= 'SIN')
cs.setreferencecode('J2000')
cs.setunits(['deg','deg'])
cs.setreferencepixel(type='direction', value=[4.,4.])
cs.setreferencevalue(type='direction', value=[0.,0.])
cs.setincrement(type='direction', value=[0.0015,0.0015]) # different from above
cs.setnames(type='direction', value=['RA','Dec'])
ia.setcoordsys(cs.torecord())
ia.close()

# actual demo
ia.open('smaller.im')
cs1 = ia.coordsys()
s1 = ia.shape()
ia.close()
ia.open('bigger.im')
ia.regrid(outfile='bigger_regridded.im', shape=s1, csys=cs1.torecord(), overwrite=True, decimate=1)
ia.close()

try:
    rval = immath(imagename=['smaller.im', 'bigger.im'], expr='IM0-IM1', outfile='difference.im')
    if not rval:
        raise
except:
    print "Expected failure since grids do not agree."

try:
    rval = immath(imagename=['smaller.im', 'bigger_regridded.im'], expr='IM0-IM1', outfile='difference.im')
    if not rval:
        raise
    else:
        print "Second attempt worked since now the grids agree."
        print "Result image difference.im should contain 2.0 and in the center a 4 x 4 square of value 0."
except:
    print "Unexpected failure."

(0 vote(s))
This article was helpful
This article was not helpful

Comments (0)
Help Desk Software by Kayako Resolve