Coverage for /wheeldirectory/casa-6.7.0-12-py3.10.el8/lib/py/lib/python3.10/site-packages/casatasks/private/task_imval.py: 83%
220 statements
« prev ^ index » next coverage.py v7.6.4, created at 2024-11-01 07:19 +0000
1import os
2import numpy
4from casatools import image, regionmanager, coordsys
5from casatasks import casalog
7_ia = image( )
8_rg = regionmanager( )
11# AUTHOR: S. Jaeger
12#
13# NAME: getimaxes
14#
15# DESCRIPTION:
16# This function uses the coordinate information associated
17# with an image to find where the directional (sky) axes are,
18# the spectral axes, and the stokes axes.
19#
20# INPUT:
21# imagename string path to a file on disk.
22#
23# RETURN
24# list of four lists, [list1, list2, list3, list4 ], as follows :
25# list1: ['axis num of 1st sky axis', 'Name of axis' ]
26# list2: ['axis num of 2nd sky axis', 'Name of axis' ]
27# list3: ['axis num of spectral axis', 'Spectral' ]
28# list4: ['axis num of stokes axis', 'Stokes' ]
30def getimaxes(imagename):
31 """
32 Open an image file, looking at its coordinate system information
33 to determine which axes are directional, linear, spectral, and
34 the stokes axies.
36 The return list or lists contains the axis numbers and names in
37 the following order:
38 1. Directional or Linear
39 2. Directional or Linear
40 3. Spectral
41 4. Stokes
43 Note that if an axis type is not found an empty list is returned
44 for that axis.
45 """
47 # Get the images coord. sys.
48 csys=None
49 _ia.open( imagename )
50 csys=_ia.coordsys()
52 # Find where the directional and channel axies are
53 # Save the internal placement of the axies in a list
54 # (which will be in the following order:
55 # direction1: RA, Longitude, Linear, el, ..
56 # direction2: DEC, Lattitude, Linear, az, ..
57 # spectral:
58 # stokes: I or V
59 axisTypes=csys.axiscoordinatetypes()
60 axisNames=csys.names()
62 # Note that we make a potentially dangerous assumption here
63 # that the first directional access is always RA, but it
64 # may not be. The names given to the axies are completely
65 # arbitrary and can not be used to determine one axis from
66 # another.
67 # TODO throw exception??? if we find extra axies or
68 # unrecognized axies.
69 retValue=[['',''],['',''],['',''],['','']]
70 foundFirstDir=False
71 for i in range( len( axisTypes ) ):
72 if ( axisTypes[i]=='Direction' or axisTypes[i]=='Linear' ):
73 if ( not foundFirstDir ) :
74 retValue[0]=[i,axisNames[i]]
75 foundFirstDir=True
76 else:
77 retValue[1]=[i,axisNames[i]]
78 elif ( axisTypes[i]=='Spectral' ) :
79 retValue[2]=[i,axisNames[i]]
80 elif ( axisTypes[i]=='Stokes' ) :
81 retValue[3]=[i,axisNames[i]]
83 if ( csys != None ):
84 del csys
85 if ( _ia.isopen() ):
86 _ia.close()
87 return retValue
90# The function that handles the imval task.
91def imval(imagename, region, box, chans, stokes):
92 myia = image()
93 mycsys = coordsys()
94 try:
95 # Blank return value.
96 retValue = { 'blc':[], 'trc':[], 'unit': '', 'data': [], 'mask': []}
97 casalog.origin('imval')
99 try:
100 axes=getimaxes(imagename)
101 except:
102 raise RuntimeError("Unable to determine the axes of image: "+imagename)
105 # Get rid of any white space in the parameters
106 region = region.replace(' ', '' )
107 box = box.replace( ' ', '' )
108 chans = chans.replace( ' ', '' )
109 stokes = stokes.replace( ' ','' )
111 # Default for the chans and stokes parameter is all when the
112 # aren't given. The easy way is to set them to -1, and let
113 # the imageregion.py code handle it.
114 if ( len(chans) < 1 ):
115 chans='-1'
117 if ( len(stokes) < 1 ):
118 stokes='-1'
121 # If the user hasn't specified any region information then
122 # find the very last data point -- what ia.getpixelvalue
123 # defaults too.
125 if ( len(box)<1 and len(chans)<1 and len(stokes)<1 and len(region)<1 ):
126 try:
127 myia.open(imagename)
129 # Get the default pixelvalue() at the referencepixel pos.
130 csys=myia.coordsys()
131 ref_values = csys.referencepixel()['numeric']
132 point=[]
133 for val in ref_values.tolist():
134 point.append( int(round(val) ) )
135 casalog.post( 'Getting data value at point '+str(point), 'NORMAL' )
136 results = myia.pixelvalue(point)
138 retValue = _imval_process_pixel( results, point )
139 retValue['axes']=axes
140 casalog.post( 'imval task complete for point'+str(point), 'NORMAL1' )
141 return retValue
142 finally:
143 myia.done()
145 # If the box parameter only has two value then we copy
146 # them.
147 if ( box.count(',') == 1 ):
148 box = box + ','+ box
150 # If we are finding the value at a single point this
151 # is a special case and we use ia.getpixelvalue()
153 singlePt = _imval_get_single( box, chans, stokes, axes )
154 if ( len( singlePt ) == 4 and singlePt.count( -1 ) < 1 ):
155 try:
156 casalog.post( 'Getting data value at point '+str(singlePt), 'NORMAL' )
157 myia.open( imagename )
158 results = myia.pixelvalue( singlePt )
159 retValue = _imval_process_pixel( results, singlePt )
160 retValue['axes']=axes
161 casalog.post( 'imval task complete for point '+str(singlePt), 'NORMAL1' )
162 return retValue
163 finally:
164 myia.done()
167 # If we've made it here then we are finding the stats
168 # of a region, which could be a set of single points.
169 axes=getimaxes(imagename)
170 statAxes=[]
171 if ( len(box)>0 ):
172 statAxes.append(axes[0][0])
173 statAxes.append(axes[1][0])
174 if ( len(chans)>0 ):
175 statAxes.append(axes[2][0])
177 # If we get to this point and find that nothing was
178 # given for the box parameter we use the reference
179 # pixel values.
180 myia.open(imagename)
181 mycsys = myia.coordsys()
183 if ( len(box) == 0 and len(region) == 0):
184 ctypes = mycsys.axiscoordinatetypes()
185 ndir = 0
186 nlin = 0
187 for ctype in ctypes:
188 if ctype == 'Direction':
189 ndir += 1
190 elif ctype == 'Linear':
191 nlin += 1
192 if ndir == 2 or nlin == 2:
193 try:
194 ref_values = mycsys.referencepixel()['numeric']
195 values = ref_values.tolist()
196 box = str(int(round(values[axes[0][0]])))+','\
197 + str(int(round(values[axes[1][0]])))+',' \
198 + str(int(round(values[axes[0][0]])))+','\
199 +str(int(round(values[axes[1][0]])))
200 except:
201 raise Exception("Unable to find the size of the input image.")
203 # Because the help file says -1 is valid, apparently that's supported functionality, good grief
205 if box.startswith("-1"):
206 box = ""
207 if chans.startswith("-1"):
208 chans = ""
209 if stokes.startswith("-1"):
210 stokes = ""
211 reg = _rg.frombcs(
212 mycsys.torecord(), myia.shape(), box, chans,
213 stokes, "a", region
214 )
217 # Now that we know which axes we are using, and have the region
218 # selected, lets get that stats! NOTE: if you have axes size
219 # greater then 0 then the maxpos and minpos will not be displayed
220 if ( 'regions' in reg ):
221 casalog.post( "Complex region found, only processing the first"\
222 " SIMPLE region found", "WARN" )
223 reg=reg['regions']['*1']
224 retValue = _imval_getregion( imagename, reg )
225 retValue['axes']=axes
227 casalog.post( 'imval task complete for region bound by blc='+str(retValue['blc'])+' and trc='+str(retValue['trc']), 'NORMAL1' )
228 return retValue
230 finally:
231 myia.done()
232 mycsys.done()
234#
235# Take the results from the ia.pixelvalue() function and
236# the position given to the function and turn the results
237# into the desired values; blc, trc, data, and mask
238#
239def _imval_process_pixel( results, point ):
240 retvalue={}
241 # First check that the results are a dictionary type and that
242 # it contains the key/value pairs we expect.
243 if ( not isinstance( results, dict ) ):
244 casalog.post( "ia.pixelvalue() has returned erroneous data, Python dictionary type expectd.", "WARN" )
245 casalog.post( "Value returned is: "+str(results), "SEVERE" )
246 return retvalue
248 if ( 'mask' not in results ):
249 casalog.post( "ia.pixelvalue() has returned unexpected results, no mask value present.", "SEVERE" )
250 return retvalue
252 if ( 'value' not in results or 'unit' not in results['value'] or 'value' not in results['value'] ):
253 casalog.post( "ia.pixelvalue() has returned unexpected results, data value absent or ill-formed.", "SEVERE" )
254 return retvalue
256 retValue={
257 'blc':point, 'trc': point, 'unit':results['value']['unit'],
258 'data': numpy.array([results['value']['value']]),
259 'mask': numpy.array([results['mask']])
260 }
261 return retValue
263#
264# Give the box, channel, and stokes values find out
265# if we are getting the data from a single point in the
266# image, if so then return that point.
267def _imval_get_single( box, chans, stokes, axes ):
268 # If we have more then one point then return an empty list.
269 try:
270 junk=int(chans)
271 junk=int(stokes)
272 except:
273 return []
274 if ( box.count(';')>0 ):
275 return []
277 # If the channel wasn't specified use the first one only.
278 if ( len( chans ) < 1 ):
279 #chans=0
280 return []
282 # If the stokes values weren't specified use the first one only.
283 if ( len( stokes ) < 1 ):
284 #stokes=0
285 return[]
287 # The two x values and two y values need to be the same if its
288 # a single point. There may be only two value if its a single
289 # point too.
290 x=-1
291 y=-1
292 box=box.split(',')
293 if ( len(box) == 2 ):
294 x=int(box[0])
295 y=int(box[1])
296 elif ( len(box) == 4 and box[0]==box[2] and box[1]==box[3]):
297 x=int(box[0])
298 y=int(box[1])
299 else:
300 # We have more then one point, return empty list.
301 return []
303 retvalue=[-1,-1,-1,-1]
305 retvalue[axes[0][0]]=x
306 retvalue[axes[1][0]]=y
307 retvalue[axes[2][0]]=int(chans)
308 retvalue[axes[3][0]]=int(stokes)
309 return retvalue
311#
312# Use the ia.getregion() function to construct the requested data.
313def _imval_getregion( imagename, region):
314 retvalue= {}
315 myia = image()
316 try:
317 # Open the image for processing!
318 myia.open(imagename)
319 # Find the array of data and mask values.
320 data_results=myia.getregion( region=region, dropdeg=True, getmask=False )
321 mask_results=myia.getregion( region=region, dropdeg=True, getmask=True )
323 # Find the bounding box of the region so we can report it back to
324 # the user.
325 bbox = myia.boundingbox( region=region )
327 if ( 'blc' not in bbox ):
328 casalog.post( "ia.boundingbox() has returned unexpected results, blc value absent.", "SEVERE" )
329 myia.done()
330 return retvalue
331 if ( 'trc' not in bbox ):
332 casalog.post( "ia.boundingbox() has returned unexpected results, trc value absent.", "SEVERE" )
333 myia.done()
334 return retvalue
336 # get the pixel coords
337 mycsys = myia.coordsys()
338 myarrays = _imval_iterate(bbox['blc'], bbox['trc'])
339 mycoords = mycsys.toworldmany(myarrays)
340 outcoords = _imval_redo(data_results.shape, mycoords['numeric'])
342 avalue = myia.pixelvalue( bbox['blc'].tolist() )
343 if ( 'value' not in avalue or 'unit' not in avalue['value'] ):
344 casalog.post( "ia.pixelvalue() has returned unexpected results, data value absent or ill-formed.", "SEVERE" )
345 myia.done()
346 return retvalue
348 retvalue={
349 'blc':bbox['blc'].tolist(),'trc':bbox['trc'].tolist(),
350 'unit':avalue['value']['unit'], 'data':data_results,
351 'mask': mask_results, 'coords': outcoords
352 }
353 finally:
354 myia.done()
355 return retvalue
357def _imval_iterate(begins, ends, arrays=None, place=None, depth=0, count=None):
358 if (depth == 0):
359 count = [0]
360 mylist = []
361 diff = numpy.array(ends) - numpy.array(begins) + 1
362 prod = diff.prod()
363 for i in range(len(begins)):
364 mylist.append(numpy.zeros([prod]))
365 arrays = numpy.array(mylist)
366 for i in range(begins[depth], ends[depth] + 1):
367 if (depth == 0):
368 tmpplace = []
369 for j in range(len(begins)):
370 tmpplace.append(0)
371 place = numpy.array(tmpplace)
372 place[depth] = i
373 if (depth == len(begins) - 1):
374 for k in range(depth + 1):
375 arrays[k][count[0]] = place[k]
376 count[0] = count[0] + 1
377 else:
378 mydepth = depth + 1
379 _imval_iterate(begins, ends, arrays, place, mydepth, count)
380 return arrays
382def _imval_redo(shape, arrays):
383 list_of_arrays = []
384 for x in range(arrays[0].size):
385 mylist = []
386 for arr in arrays:
387 mylist.append(arr[x])
388 list_of_arrays.append(numpy.array(mylist))
389 array_of_arrays = numpy.array(list_of_arrays)
390 # because shape is an immutable tuple
391 newshape = list(shape)
392 newshape.append(array_of_arrays.shape[1])
393 return array_of_arrays.reshape(newshape)