Using external cut in a script

brentdutra wrote on Dec 06, 2010

So I've run into a stumbling block with the syntax that the iraf cl scripting language uses. I've created a script which takes in raw images and coordinates as input and generates a series of output images per my specifications. The whole idea is such: I'm looking for tidal distortions in elliptical galaxies which are members of compact groups. I have the coordinate data for such groups. When I use my known coordinates (from SDSS) I get a slew of images output which are then sorted (using a python script I have) into individual group folders, however each group may contain 1 or more images and/or 1 or more objects of interest. Using my cl script I then would like to "catalog" (using ccfind against the images and source coordinates) and fit isophotal ellipses to the galaxies, then create a model and subtract it from my image data. Ok, well my script below essentially does that; however one line does not function correctly....

match( s1,"*.match") | cut("-f 1 -d:"," > objlist")

This line works fine under iraf cl calling it as an individual "task", but when run in the script it just halts and the script stalls. I have added cut to my external cl packages, so I know it obviously works in some sense. Am I using an incorrect syntax? Is there a simpler way to do what I'm doing? And I certainly know there are more efficient ways of doing certain parts of my script but they work for now.

Now I have included the script below but the reason I am resorting to using cut is that I generate a file which has file names prepended to coordinates; I need to separate the file names from the coordinates
(eg.
isophote> match(s1,"*.match")
fpC-005598-r3-0180.fit.match:00:38:54.69 07:03:21.78 87.036 445.801
fpC-005641-r2-0061.fit.match:00:38:54.69 07:03:21.78 1994.828 635.513
fpC-005656-r3-0080.fit.match:00:38:54.69 07:03:21.78 96.205 644.781

--->
fpC-005598-r3-0180.fit.match
fpC-005641-r2-0061.fit.match
fpC-005656-r3-0080.fit.match
)
Perhaps there's a better way to do this, that I'm missing?

#real x,y
int howMany
int lineNumber = 1
int naxis1
int naxis2
int objects, currentObj
real centerDistance
real dist1,dist2
string filename
string someReferenceImage
string currentFile
#string currentLine
string ra,dec
struct *dlist
struct *objList

#package isophote


###Run over objects from HCG 1 to 100-----------------------------------

for (j=1; j <= =100; j+=1) {
	print(j)
		if (access("./hcg" // j)) {
			cd ("./hcg" // j) 
			
###Check for all files in current group directory----------------------

			files *.fit > files
			count files | scan(howMany)
			
			list = "files"
			while ( fscan(list, currentFile) != EOF) {

			ccfind ("/Research/scripts/hcgccfind.coords", currentFile //".match", currentFile, lngcolumn=1, latcolumn=2, 				lngunits="",latunits="",insystem="j2000", usewcs=yes, verbose=no)

			match("^[0-9]", currentFile //".match") | count | scan(objects)
			print("Number of Objects Found: ",objects)
			print(currentFile, ,'\t', objects, >> "objects")

			#!cat ./*.match | grep -h ^[0-9] | cut -c -30 | uniq >> uniquecoords
			
			}
			
			!cat ./*.match | grep -h ^[0-9] | cut -c -30 | sort -u > uniquecoords

			#for( ToLoopOver = objects; ToLoopOver >=1; ToLoopOver -=1 ) {
			objList = "uniquecoords"
			while ( fscan (objList, ra, dec) !=EOF ) {

				s1 = ra#dec
				match( s1,"*.match") | cut("-f 1 -d:"," > objlist")
				#!cat *.match | grep -h ^[0-9] | cut -c -30 | sort -u | xargs -i grep {} *.match | cut -f 1 -d: > objlist

				count objlist | scan(howMany)

				if (howMany > 1) {
					list = "objlist"
					while ( fscan(list, currentFile) != EOF) {
						print ("Working...")

						print (currentFile)
						keypar(currentFile,"naxis1")
						naxis1 = int(keypar.value)
						keypar(currentFile,"naxis2")
						naxis2 = int(keypar.value)

						#!grep ^[0-9] *.match >> matchedxycoords
						match("^[0-9]", currentFile //".match") | count | scan(objects)
						print("Number of Objects Found: ",objects)
						print(currentFile, ,'\t', objects, >> "objects")
						fields(currentFile //".match", "3,4",lines=lineNumber+9) | scan(x,y)

						centerDistance = sqrt((abs((naxis1/2.)-x)**2)+(abs((naxis2/2.)-y)**2))

						print(currentFile, '\t',x,'\t',y,'\t',centerDistance, >> "distance")
					
						}

###If more than one image in current group for current object, find the best, set as refereence image for SRegister-------------------
					dist2 = 10000.
					dlist = "distance"
					while( fscan(dlist, currentFile,s1,s2,dist1) !=EOF){
					
						#fields(currentLine,"4") | scan(dist1)

						if (dist1 < dist2){
							dist2 = dist1
							someReferenceImage = currentFile
							x = real(s1)
							y = real(s2)
						}
					}

###Finding most centered image above for sregister below---------------------------------------------------------------------------------



###Generate a single combined image if more than one file for image and group------------------------------------------------------------
					sregister ("@files", someReferenceImage, "sregd_"//"@files", wcs="world", interpolant="linear",boundary="constant", 						constant=-10000., wcsinherit=yes)

					imcombine ("sregd_"//"@files", "hcg" //j //"_combined.fits", combine="median", reject="avsigclip", lthreshold=-100.)
					filename = "hcg"//j//"_combined.fits"
###Output a single combined image from above---------------------------------------------------------------------------------------------------------
				}	

		
###If only one image find the pixel space coordinates of the single object and prepare variables for "Interesting Stuff"--------------------
				else {

				!grep ^[0-9] objlist > matchedxycoords
				fields("matchedxycoords", "3,4", lines = lineNumber) | scan(x,y)
				files *.fit | scan(filename)

				}

###Outputing all necessary variables for "Interesting Stuff"-----------------------------------------------------------------------------------


###Do all the interesting stuff, fit isophotal ellipses to object coords in pixel space, generate a *.fits model, subtract model from source data

				geompar.x0 = x
				geompar.y0 = y
				lpar geompar
			
				ellipse (filename, "hcg" //j //"_ellipse.tbl", verbose=no)

				bmodel ("hcg" //j //"_ellipse.tbl", "hcg" //j //"_model.fits",  minsma=1., maxsma=1., backgr=0., interp="spline")

				imarith (filename, "-", "hcg"//j//"_model.fits", "hcg"//j//"_removed.fits", title="", divzero=0., hparams="", pixtype="", 					calctype="", verbose=no, noact=no)
###Interesting stuff above-------------------------------------------------------------------------------------------------------------------------

			}

			cd("../")

		}

;
}

---

Last post on Dec 06, 2010