Commits

David Mehringer authored ea2e9c50578
        New Development: No             JIRA Issue: CAS-5261          Ready for Test: Yes      Interface Changes: No What Interface Changed: Please list interface changes          Test Programs: List test programs   Put in Release Notes: Yes/No              Module(s): Module Names change impacts.            Description: Added exception handling, dimensionality checking, and documentation

No tags

gcwrap/python/scripts/task_boxit.py

Modified
1 1 from taskinit import *
2 2 import numpy
3 3 import sys
4 4 import re
5 5
6 6 # Writes out regions above threshold to regionfile+'.box'
7 7
8 8 def boxit(imagename, regionfile, threshold, maskname, chanrange, polrange, minsize, diag, boxstretch, overwrite):
9 9
10 10 casalog.origin('boxit')
11 -
12 - if not(regionfile):
13 - regionfile = imagename + '.box'
14 - if not regionfile.endswith('.box'):
15 - regionfile = regionfile + '.box'
16 -
17 - if not(overwrite):
18 - if(os.path.exists(regionfile)):
19 - casalog.post('file "' + regionfile + '" already exists.', 'WARN')
20 - return
21 - if(maskname and os.path.exists(maskname)):
22 - casalog.post('output mask "' + maskname + '" already exists.', 'WARN')
11 + myia = iatool()
12 + try:
13 + if not(regionfile):
14 + regionfile = imagename + '.box'
15 + if not regionfile.endswith('.box'):
16 + regionfile = regionfile + '.box'
17 +
18 + if not(overwrite):
19 + if(os.path.exists(regionfile)):
20 + casalog.post('file "' + regionfile + '" already exists.', 'WARN')
21 + return
22 + if(maskname and os.path.exists(maskname)):
23 + casalog.post('output mask "' + maskname + '" already exists.', 'WARN')
24 + return
25 +
26 + # If no units, assume mJy for consistency with auto/clean tasks.
27 + # But convert to Jy, because that's what units the images are in.
28 + threshold = qa.getvalue(qa.convert(qa.quantity(threshold,'mJy'),'Jy'))[0]
29 + casalog.post("Setting threshold to " + str(threshold) + "Jy", "INFO")
30 + newIsland = numpy.zeros(1, dtype=[('box','4i4'),('npix','i4')])
31 + # Find all pixels above the threshold
32 + myia.open(imagename)
33 + if len(myia.shape()) != 4:
34 + raise Exception("Only 4D images with direction, spectral, and stokes coordinates are supported")
35 + # CAS-2059 escape characters in image name that will confuse the lattice expression processor
36 + escaped_imagename = imagename
37 + for escapeme in ['-', '+', '*', '/' ]:
38 + escaped_imagename = re.sub("[" + escapeme + "]", "\\" + escapeme, escaped_imagename)
39 + mask = escaped_imagename+'>'+str(threshold)
40 + fullmask = myia.getregion(mask=mask, getmask=True)
41 + if not(fullmask.max()):
42 + casalog.post('Maximum flux in image is below threshold.', 'WARN')
23 43 return
24 -
25 - # If no units, assume mJy for consistency with auto/clean tasks.
26 - # But convert to Jy, because that's what units the images are in.
27 - threshold = qa.getvalue(qa.convert(qa.quantity(threshold,'mJy'),'Jy'))[0]
28 - casalog.post("Setting threshold to " + str(threshold) + "Jy", "INFO")
29 -
30 - newIsland = numpy.zeros(1, dtype=[('box','4i4'),('npix','i4')])
31 - # Find all pixels above the threshold
32 - ia.open(imagename)
33 - # CAS-2059 escape characters in image name that will confuse the lattice expression processor
34 - escaped_imagename = imagename
35 - for escapeme in ['-', '+', '*', '/' ]:
36 - escaped_imagename = re.sub("[" + escapeme + "]", "\\" + escapeme, escaped_imagename)
37 - mask = escaped_imagename+'>'+str(threshold)
38 - fullmask = ia.getregion(mask=mask, getmask=True)
39 - if not(fullmask.max()):
40 - casalog.post('Maximum flux in image is below threshold.', 'WARN')
41 - return
42 - writemask = bool(maskname)
43 - csys = ia.coordsys()
44 -
45 - shape = fullmask.shape
46 -
47 - nx = int(shape[0])
48 - ny = int(shape[1])
49 - n2 = n3 = 1
50 - if len(shape)==3:
51 - n2 = int(shape[2])
52 - if len(shape)==4:
53 - n2 = int(shape[2])
54 - n3 = int(shape[3])
55 - #print 'nx:', nx, 'ny:', ny, 'n2:', n2, 'n3:', n3
56 -
57 - #casa generated images always 4d and in order of [ra, dec, stokes, freq]
58 - #other images can be in the order of [ra, dec, freq, stokes]
59 - nms = csys.names()
60 - #axes=[]
61 - #for i in xrange(len(nms)):
62 - # if nms[i]=='Right Ascension':
63 - # axes[0]=i
64 - # if nms[i]=='Declination':
65 - # axes[1]=i
66 - # if nms[i]=='Stokes'
67 - # axes[2]=i
68 - # if nms[i]=='Frequency':
69 - # axes[3]=i
70 -
71 - chmax=n3
72 - pomax=n2
73 - chmin=0
74 - pomin=0
75 - if len(nms)==4 and nms[3]=='Stokes':
76 - chmax=n2
77 - pomax=n3
78 -
79 - if chanrange:
80 - try:
81 - if str.count(chanrange, '~') == 1:
82 - ch1=int(str.split(chanrange, '~')[0])
83 - ch2=int(str.split(chanrange, '~')[1])
84 - else:
85 - ch1=int(chanrange)
86 - ch2=ch1
87 - if ch1 > ch2 or ch1 < 0 or ch2 > chmax:
88 - casalog('invalid channel range', "SEVERE")
89 - return
90 - if ch1>chmin:
91 - chmin = ch1
92 - if ch2<chmax:
93 - chmax = ch2
94 - except:
95 - casalog('bad format for chanrange', "SEVERE")
96 - return
97 - #print 'chmin', chmin, 'chmax', chmax
98 -
99 - if polrange:
100 - try:
101 - if str.count(polrange, '~') == 1:
102 - po1=int(str.split(polrange, '~')[0])
103 - po2=int(str.split(polrange, '~')[1])
104 - else:
105 - po1=int(polrange)
106 - po2=po1
107 - if po1 > po2 or po1 < 0 or po2 > pomax:
108 - casalog( 'invalid stokes range', "SEVERE")
109 - return
110 - if po1>pomin:
111 - pomin = po1
112 - if po2<pomax:
113 - pomax = po2
114 - except:
115 - casalog( 'bad format for polrange', "SEVERE")
116 - return
117 - #print 'pomin', pomin, 'pomax', pomax
118 -
119 - n1=pomin
120 - n2=pomax
121 - n3=chmax
122 - n4=chmin
123 - if len(nms)==4 and nms[3]=='Stokes':
124 - n1=chmin
125 - n2=chmax
126 - n3=pomax
127 - n4=pomin
128 - #print 'n1:', n1, 'n2:', n2, 'n4:', n4, 'n3:', n3
44 + writemask = bool(maskname)
45 + csys = myia.coordsys()
129 46
130 - f = open(regionfile, 'w')
131 - totregions = 0
132 - outputmask = []
133 - if writemask:
134 - outputmask = ia.getchunk()
135 - outputmask.fill(False)
136 -
137 - for i3 in xrange(n4, n3):
138 - for i2 in xrange(n1, n2):
139 - regions = {}
140 - boxRecord = {}
141 - if len(shape)==2:
142 - mask = fullmask
143 - if len(shape)==4:
144 - mask = fullmask[:,:,i2,i3].reshape(nx, ny)
145 - islands = []
146 - pos = numpy.unravel_index(mask.argmax(), mask.shape)
147 - #print pos
148 - while(mask[pos]):
149 - # found pixel in new island
150 - island = {}
151 - island['box'] = [pos[0], pos[1], pos[0], pos[1]]
152 - island['npix'] = 1
153 - find_nearby_island_pixels(island, mask, pos, shape[0], shape[1], diag)
154 - islands.append(island)
47 + shape = fullmask.shape
48 + nx = int(shape[0])
49 + ny = int(shape[1])
50 + n2 = n3 = 1
51 + if len(shape)==3:
52 + n2 = int(shape[2])
53 + if len(shape)==4:
54 + n2 = int(shape[2])
55 + n3 = int(shape[3])
56 +
57 + #casa generated images always 4d and in order of [ra, dec, stokes, freq]
58 + #other images can be in the order of [ra, dec, freq, stokes]
59 + nms = csys.names()
60 +
61 + chmax=n3
62 + pomax=n2
63 + chmin=0
64 + pomin=0
65 + if len(nms)==4 and nms[3]=='Stokes':
66 + chmax=n2
67 + pomax=n3
68 + if chanrange:
69 + try:
70 + if str.count(chanrange, '~') == 1:
71 + ch1=int(str.split(chanrange, '~')[0])
72 + ch2=int(str.split(chanrange, '~')[1])
73 + else:
74 + ch1=int(chanrange)
75 + ch2=ch1
76 + if ch1 > ch2 or ch1 < 0 or ch2 > chmax:
77 + casalog('invalid channel range', "SEVERE")
78 + return
79 + if ch1>chmin:
80 + chmin = ch1
81 + if ch2<chmax:
82 + chmax = ch2
83 + except:
84 + casalog('bad format for chanrange', "SEVERE")
85 + return
86 +
87 + if polrange:
88 + try:
89 + if str.count(polrange, '~') == 1:
90 + po1=int(str.split(polrange, '~')[0])
91 + po2=int(str.split(polrange, '~')[1])
92 + else:
93 + po1=int(polrange)
94 + po2=po1
95 + if po1 > po2 or po1 < 0 or po2 > pomax:
96 + casalog( 'invalid stokes range', "SEVERE")
97 + return
98 + if po1>pomin:
99 + pomin = po1
100 + if po2<pomax:
101 + pomax = po2
102 + except:
103 + raise Exception('bad format for polrange')
104 + n1=pomin
105 + n2=pomax
106 + n3=chmax
107 + n4=chmin
108 + if len(nms)==4 and nms[3]=='Stokes':
109 + n1=chmin
110 + n2=chmax
111 + n3=pomax
112 + n4=pomin
113 +
114 + f = open(regionfile, 'w')
115 + totregions = 0
116 + outputmask = []
117 + if writemask:
118 + outputmask = myia.getchunk()
119 + outputmask.fill(False)
120 + for i3 in xrange(n4, n3):
121 + for i2 in xrange(n1, n2):
122 + regions = {}
123 + boxRecord = {}
124 + if len(shape)==2:
125 + mask = fullmask
126 + if len(shape)==4:
127 + mask = fullmask[:,:,i2,i3].reshape(nx, ny)
128 + islands = []
155 129 pos = numpy.unravel_index(mask.argmax(), mask.shape)
156 - for record in islands:
157 - box = record['box']
158 - # check size of box
159 - if record['npix'] < minsize:
160 - continue
161 - totregions += 1
162 - # need pixel corners, not pixel centers
163 - box[0] -= boxstretch
164 - box[1] -= boxstretch
165 - box[2] += boxstretch
166 - box[3] += boxstretch
167 - # in case we used boxstretch < 0 and a one-pixel sized box:
168 - if box[0] > box[2]:
169 - box[0] += boxstretch
170 - box[2] -= boxstretch
171 - if box[1] > box[3]:
172 - box[1] += boxstretch
173 - box[3] -= boxstretch
174 - if writemask:
175 - # avoid pixels in boxes that have been stretched beyond the image limits
176 - for ii in range(max(0, box[0]), min(box[2], outputmask.shape[0] - 1)):
177 - for jj in range(max(0, box[1]), min(box[3], outputmask.shape[1] - 1)):
178 - outputmask[ii][jj][i2][i3] = True
179 - blccoord = [box[0]-0.5, box[1]-0.5, i2, i3]
180 - trccoord = [box[2]+0.5, box[3]+0.5, i2, i3]
181 - # note that the toworld() calls are likely very expensive for many boxes. But then
182 - # again, the box-finding algorithm itself seems pretty inefficient, but resource
183 - # constraints only permit a band aid fix at this time.
184 - blc = ia.toworld(blccoord, 'm')['measure']
185 - trc = ia.toworld(trccoord, 'm')['measure']
186 - # RA/Dec reference frame
187 - outstring = "worldbox " + blc['direction']['refer']
188 - # RA blc/trc
189 - outstring += " [" + quantity_to_string(blc["direction"]["m0"], "rad") + ", "
190 - outstring += quantity_to_string(trc["direction"]["m0"], "rad") + "]"
191 - # Dec blc/trc
192 - outstring += " [" + quantity_to_string(blc["direction"]["m1"], "rad") + ", "
193 - outstring += quantity_to_string(trc["direction"]["m1"], "rad") + "]"
194 - # frequency blc/trc
195 - freq = blc["spectral"]['frequency']['refer'] + " "
196 - freq += quantity_to_string(blc["spectral"]["frequency"]["m0"], "Hz", False)
197 - outstring += " ['" + freq + "', '" + freq + "']"
198 - # Stokes blc/trc
199 - outstring += " ['" + blc["stokes"] + "', '" + trc["stokes"] + "']"
200 - # add the mask flag
201 - outstring += " " + str(1)
202 - f.write(outstring + "\n")
203 - casalog.post("Wrote " + str(totregions) + " regions to file " + regionfile, 'INFO')
204 - if writemask:
205 - myia = iatool()
206 - myia.fromimage(infile=imagename, outfile=maskname, overwrite=True)
207 - myia.done()
208 - myia.open(maskname)
209 - myia.putchunk(outputmask)
210 - myia.done()
211 - f.close()
212 - return True
130 + while(mask[pos]):
131 + # found pixel in new island
132 + island = {}
133 + island['box'] = [pos[0], pos[1], pos[0], pos[1]]
134 + island['npix'] = 1
135 + find_nearby_island_pixels(island, mask, pos, shape[0], shape[1], diag)
136 + islands.append(island)
137 + pos = numpy.unravel_index(mask.argmax(), mask.shape)
138 + for record in islands:
139 + box = record['box']
140 + # check size of box
141 + if record['npix'] < minsize:
142 + continue
143 + totregions += 1
144 + # need pixel corners, not pixel centers
145 + box[0] -= boxstretch
146 + box[1] -= boxstretch
147 + box[2] += boxstretch
148 + box[3] += boxstretch
149 + # in case we used boxstretch < 0 and a one-pixel sized box:
150 + if box[0] > box[2]:
151 + box[0] += boxstretch
152 + box[2] -= boxstretch
153 + if box[1] > box[3]:
154 + box[1] += boxstretch
155 + box[3] -= boxstretch
156 + if writemask:
157 + # avoid pixels in boxes that have been stretched beyond the image limits
158 + for ii in range(max(0, box[0]), min(box[2], outputmask.shape[0] - 1)):
159 + for jj in range(max(0, box[1]), min(box[3], outputmask.shape[1] - 1)):
160 + outputmask[ii][jj][i2][i3] = True
161 + blccoord = [box[0]-0.5, box[1]-0.5, i2, i3]
162 + trccoord = [box[2]+0.5, box[3]+0.5, i2, i3]
163 + # note that the toworld() calls are likely very expensive for many boxes. But then
164 + # again, the box-finding algorithm itself seems pretty inefficient, but resource
165 + # constraints only permit a band aid fix at this time.
166 + blc = myia.toworld(blccoord, 'm')['measure']
167 + trc = myia.toworld(trccoord, 'm')['measure']
168 + # RA/Dec reference frame
169 + outstring = "worldbox " + blc['direction']['refer']
170 + # RA blc/trc
171 + outstring += " [" + quantity_to_string(blc["direction"]["m0"], "rad") + ", "
172 + outstring += quantity_to_string(trc["direction"]["m0"], "rad") + "]"
173 + # Dec blc/trc
174 + outstring += " [" + quantity_to_string(blc["direction"]["m1"], "rad") + ", "
175 + outstring += quantity_to_string(trc["direction"]["m1"], "rad") + "]"
176 + # frequency blc/trc
177 + freq = blc["spectral"]['frequency']['refer'] + " "
178 + freq += quantity_to_string(blc["spectral"]["frequency"]["m0"], "Hz", False)
179 + outstring += " ['" + freq + "', '" + freq + "']"
180 + # Stokes blc/trc
181 + outstring += " ['" + blc["stokes"] + "', '" + trc["stokes"] + "']"
182 + # add the mask flag
183 + outstring += " " + str(1)
184 + f.write(outstring + "\n")
185 + casalog.post("Wrote " + str(totregions) + " regions to file " + regionfile, 'INFO')
186 + if writemask:
187 + myia.fromimage(infile=imagename, outfile=maskname, overwrite=True)
188 + myia.done()
189 + myia.open(maskname)
190 + myia.putchunk(outputmask)
191 + myia.done()
192 + f.close()
193 + return True
194 + except Exception, instance:
195 + casalog.post( str( '*** Error *** ') + str(instance), 'SEVERE')
196 + raise
197 + finally:
198 + if myia:
199 + myia.done()
213 200
214 201 def quantity_to_string(quantity, unit=None, quotes=True):
215 202 if unit != None:
216 203 quantity = qa.convert(quantity, unit)
217 204 string = str(quantity['value']) + quantity['unit']
218 205 if quotes:
219 206 string = "'" + string + "'"
220 207 return string
221 208
222 209 def find_nearby_island_pixels(island, mask, pos, xmax, ymax, diag):

Everything looks good. We'll let you know here if there's anything you should know about.

Add shortcut