Changeset 2882 for trunk/python
- Timestamp:
- 12/20/13 17:34:30 (11 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/python/scantable.py
r2877 r2882 301 301 return res 302 302 303 def is_number(s): 304 s = s.strip() 305 res = True 306 try: 307 a = float(s) 308 res = True 309 except: 310 res = False 311 finally: 312 return res 313 314 def is_frequency(s): 315 s = s.strip() 316 return (s[-2:].lower() == "hz") 317 318 def get_freq_by_string(s): 319 if not is_frequency(s): 320 raise RuntimeError("Invalid input string.") 321 322 prefix_list = ["a", "f", "p", "n", "u", "m", ".", "k", "M", "G", "T", "P", "E"] 323 factor_list = [1e-18, 1e-15, 1e-12, 1e-9, 1e-6, 1e-3, 1.0, 1e+3, 1e+6, 1e+9, 1e+12, 1e+15, 1e+18] 324 325 s = s.strip() 326 factor = 1.0 327 328 prefix = s[-3:-2] 329 if is_number(prefix): 330 res = float(s[:-2]) 331 else: 332 res = float(s[:-3]) * factor_list[prefix_list.index(prefix)] 333 334 return res 335 336 def is_velocity(s): 337 s = s.strip() 338 return (s[-3:].lower() == "m/s") 339 340 def get_velocity_by_string(s): 341 if not is_velocity(s): 342 raise RuntimeError("Invalid input string.") 343 344 prefix_list = [".", "k"] 345 factor_list = [1e-3, 1.0] 346 347 s = s.strip() 348 factor = 1.0 349 350 prefix = s[-4:-3] 351 if is_number(prefix): 352 res = float(s[:-3]) * 1e-3 353 else: 354 res = float(s[:-4]) * factor_list[prefix_list.index(prefix)] 355 356 return res # in km/s 357 358 def get_frequency_by_velocity(restfreq, vel): 359 # vel is in unit of km/s 360 361 # speed of light 362 vel_c = 299792.458 363 364 import math 365 r = vel / vel_c 366 return restfreq * math.sqrt((1.0 - r) / (1.0 + r)) 367 368 303 369 class scantable(Scantable): 304 370 """\ … … 1623 1689 return istart,iend 1624 1690 1691 @asaplog_post_dec 1692 def parse_spw_selection(self, selectstring, restfreq=None, frame=None, doppler=None): 1693 """ 1694 Parse MS type spw/channel selection syntax. 1695 1696 Parameters: 1697 selectstring : A string expression of spw and channel selection. 1698 Comma-separated expressions mean different spw - 1699 channel combinations. Spws and channel selections 1700 are partitioned by a colon ':'. In a single 1701 selection expression, you can put multiple values 1702 separated by semicolons ';'. Both for spw and 1703 channel selection, allowed cases include single 1704 value, blank('') or asterisk('*') to specify all 1705 available values, two values connected with a 1706 tilde ('~') to specify an inclusive range. Unit 1707 strings for frequency or velocity can be added to 1708 the tilde-connected values. For channel selection 1709 expression, placing a '<' or a '>' is possible to 1710 specify a semi-infinite interval as well. 1711 1712 examples: 1713 '' or '*' = all spws (all channels) 1714 '<2,4~6,9' = Spws 0,1,4,5,6,9 (all channels) 1715 '3:3~45;60' = channels 3 to 45 and 60 in spw 3 1716 '0~1:2~6,8' = channels 2 to 6 in spws 0,1, and 1717 all channels in spw8 1718 '1.3GHz~1.5GHz' = all spws that fall in or have 1719 at least some overwrap with 1720 frequency range between 1.3GHz 1721 and 1.5GHz. 1722 '1.3GHz~1.5GHz:1.3GHz~1.5GHz' = channels that 1723 fall between the 1724 specified frequency 1725 range in spws that 1726 fall in or have 1727 overwrap with the 1728 specified frequency 1729 range. 1730 '1:-200km/s~250km/s' = channels that fall between 1731 the specified velocity range 1732 in spw 1. 1733 Returns: 1734 A dictionary of selected (valid) spw and masklist pairs, 1735 e.g. {'0': [[50,250],[350,462]], '2': [[100,400],[550,974]]} 1736 """ 1737 if not isinstance(selectstring, str): 1738 asaplog.post() 1739 asaplog.push("Expression of spw/channel selection must be a string.") 1740 asaplog.post("ERROR") 1741 1742 orig_unit = self.get_unit() 1743 self.set_unit('channel') 1744 1745 orig_restfreq = self.get_restfreqs() 1746 orig_restfreq_list = [] 1747 for i in orig_restfreq.keys(): 1748 if len(orig_restfreq[i]) == 1: 1749 orig_restfreq_list.append(orig_restfreq[i][0]) 1750 else: 1751 orig_restfreq_list.append(orig_restfreq[i]) 1752 1753 orig_coord = self._getcoordinfo() 1754 orig_frame = orig_coord[1] 1755 orig_doppler = orig_coord[2] 1756 1757 if restfreq is None: restfreq = orig_restfreq_list 1758 if frame is None: frame = orig_frame 1759 if doppler is None: doppler = orig_doppler 1760 1761 self.set_restfreqs(restfreq) 1762 self.set_freqframe(frame) 1763 self.set_doppler(doppler) 1764 1765 valid_ifs = self.getifnos() 1766 1767 comma_sep = selectstring.split(",") 1768 res = {} 1769 1770 for cms_elem in comma_sep: 1771 if (cms_elem.strip() == ""): continue 1772 1773 colon_sep = cms_elem.split(":") 1774 1775 if (len(colon_sep) > 2): 1776 raise RuntimeError("Invalid selection expression: more than two colons!") 1777 1778 # parse spw expression and store result in spw_list. 1779 # allowed cases include '', '*', 'a', '<a', '>a', 'a~b', 1780 # 'a*Hz~b*Hz' (where * can be '', 'k', 'M', 'G' etc.), 1781 # 'a*m/s~b*m/s' (where * can be '' or 'k') and also 1782 # several of the above expressions connected with ';'. 1783 1784 spw_list = [] 1785 1786 semicolon_sep = colon_sep[0].split(";") 1787 1788 for scs_elem in semicolon_sep: 1789 scs_elem = scs_elem.strip() 1790 1791 lt_sep = scs_elem.split("<") 1792 gt_sep = scs_elem.split(">") 1793 ti_sep = scs_elem.split("~") 1794 1795 lt_sep_length = len(lt_sep) 1796 gt_sep_length = len(gt_sep) 1797 ti_sep_length = len(ti_sep) 1798 1799 len_product = lt_sep_length * gt_sep_length * ti_sep_length 1800 1801 if (len_product > 2): 1802 # '<', '>' and '~' must not coexist in a single spw expression 1803 1804 raise RuntimeError("Invalid spw selection.") 1805 1806 elif (len_product == 1): 1807 # '', '*', or single spw number. 1808 1809 if (scs_elem == "") or (scs_elem == "*"): 1810 spw_list = valid_ifs[:] # deep copy 1811 1812 else: # single number 1813 try: 1814 #checking if the given number is valid for spw ID 1815 idx = valid_ifs.index(int(scs_elem)) 1816 spw_list.append(valid_ifs[idx]) 1817 1818 except: 1819 asaplog.post() 1820 asaplog.push("Wrong spw number (" + scs_elem + ") given. ignored.") 1821 asaplog.post("WARNING") 1822 1823 else: # (len_product == 2) 1824 # namely, one of '<', '>' or '~' appers just once. 1825 1826 if (lt_sep_length == 2): # '<a' 1827 if is_number(lt_sep[1]): 1828 for i in valid_ifs: 1829 if (i < float(lt_sep[1])): 1830 spw_list.append(i) 1831 1832 else: 1833 RuntimeError("Invalid spw selection.") 1834 1835 elif (gt_sep_length == 2): # '>a' 1836 if is_number(gt_sep[1]): 1837 for i in valid_ifs: 1838 if (i > float(gt_sep[1])): 1839 spw_list.append(i) 1840 1841 else: 1842 RuntimeError("Invalid spw selection.") 1843 1844 else: # (ti_sep_length == 2) where both boundaries inclusive 1845 expr0 = ti_sep[0].strip() 1846 expr1 = ti_sep[1].strip() 1847 1848 if is_number(expr0) and is_number(expr1): 1849 # 'a~b' 1850 expr_pmin = min(float(expr0), float(expr1)) 1851 expr_pmax = max(float(expr0), float(expr1)) 1852 for i in valid_ifs: 1853 if (expr_pmin <= i) and (i <= expr_pmax): 1854 spw_list.append(i) 1855 1856 elif is_frequency(expr0) and is_frequency(expr1): 1857 # 'a*Hz~b*Hz' 1858 expr_f0 = get_freq_by_string(expr0) 1859 expr_f1 = get_freq_by_string(expr1) 1860 1861 for coord in self._get_coordinate_list(): 1862 expr_p0 = coord['coord'].to_pixel(expr_f0) 1863 expr_p1 = coord['coord'].to_pixel(expr_f1) 1864 expr_pmin = min(expr_p0, expr_p1) 1865 expr_pmax = max(expr_p0, expr_p1) 1866 1867 spw = coord['if'] 1868 pmin = 0.0 1869 pmax = float(self.nchan(spw) - 1) 1870 1871 if ((expr_pmax - pmin)*(expr_pmin - pmax) <= 0.0): 1872 spw_list.append(spw) 1873 1874 elif is_velocity(expr0) and is_velocity(expr1): 1875 # 'a*m/s~b*m/s' 1876 expr_v0 = get_velocity_by_string(expr0) 1877 expr_v1 = get_velocity_by_string(expr1) 1878 expr_vmin = min(expr_v0, expr_v1) 1879 expr_vmax = max(expr_v0, expr_v1) 1880 1881 for coord in self._get_coordinate_list(): 1882 spw = coord['if'] 1883 1884 pmin = 0.0 1885 pmax = float(self.nchan(spw) - 1) 1886 1887 vel0 = coord['coord'].to_velocity(pmin) 1888 vel1 = coord['coord'].to_velocity(pmax) 1889 1890 vmin = min(vel0, vel1) 1891 vmax = max(vel0, vel1) 1892 1893 if ((expr_vmax - vmin)*(expr_vmin - vmax) <= 0.0): 1894 spw_list.append(spw) 1895 1896 else: 1897 # cases such as 'aGHz~bkm/s' are not allowed now 1898 raise RuntimeError("Invalid spw selection.") 1899 1900 # parse channel expression and store the result in crange_list. 1901 # allowed cases include '', 'a~b', 'a*Hz~b*Hz' (where * can be 1902 # '', 'k', 'M', 'G' etc.), 'a*m/s~b*m/s' (where * can be '' or 'k') 1903 # and also several of the above expressions connected with ';'. 1904 1905 for spw in spw_list: 1906 pmin = 0.0 1907 pmax = float(self.nchan(spw) - 1) 1908 1909 if (len(colon_sep) == 1): 1910 # no expression for channel selection, 1911 # which means all channels are to be selected. 1912 crange_list = [[pmin, pmax]] 1913 1914 else: # (len(colon_sep) == 2) 1915 crange_list = [] 1916 1917 found = False 1918 for i in self._get_coordinate_list(): 1919 if (i['if'] == spw): 1920 coord = i['coord'] 1921 found = True 1922 break 1923 1924 if not found: 1925 raise RuntimeError("Invalid spw value.") 1926 1927 semicolon_sep = colon_sep[1].split(";") 1928 for scs_elem in semicolon_sep: 1929 scs_elem = scs_elem.strip() 1930 1931 ti_sep = scs_elem.split("~") 1932 ti_sep_length = len(ti_sep) 1933 1934 if (ti_sep_length > 2): 1935 raise RuntimeError("Invalid channel selection.") 1936 1937 elif (ti_sep_length == 1): 1938 if (scs_elem == "") or (scs_elem == "*"): 1939 # '' and '*' for all channels 1940 crange_list = [[pmin, pmax]] 1941 break 1942 elif (is_number(scs_elem)): 1943 # single channel given 1944 crange_list.append([float(scs_elem), float(scs_elem)]) 1945 else: 1946 raise RuntimeError("Invalid channel selection.") 1947 1948 else: #(ti_sep_length == 2) 1949 expr0 = ti_sep[0].strip() 1950 expr1 = ti_sep[1].strip() 1951 1952 if is_number(expr0) and is_number(expr1): 1953 # 'a~b' 1954 expr_pmin = min(float(expr0), float(expr1)) 1955 expr_pmax = max(float(expr0), float(expr1)) 1956 1957 elif is_frequency(expr0) and is_frequency(expr1): 1958 # 'a*Hz~b*Hz' 1959 expr_p0 = coord.to_pixel(get_freq_by_string(expr0)) 1960 expr_p1 = coord.to_pixel(get_freq_by_string(expr1)) 1961 expr_pmin = min(expr_p0, expr_p1) 1962 expr_pmax = max(expr_p0, expr_p1) 1963 1964 elif is_velocity(expr0) and is_velocity(expr1): 1965 # 'a*m/s~b*m/s' 1966 restf = self.get_restfreqs().values()[0][0] 1967 expr_f0 = get_frequency_by_velocity(restf, get_velocity_by_string(expr0)) 1968 expr_f1 = get_frequency_by_velocity(restf, get_velocity_by_string(expr1)) 1969 expr_p0 = coord.to_pixel(expr_f0) 1970 expr_p1 = coord.to_pixel(expr_f1) 1971 expr_pmin = min(expr_p0, expr_p1) 1972 expr_pmax = max(expr_p0, expr_p1) 1973 1974 else: 1975 # cases such as 'aGHz~bkm/s' are not allowed now 1976 raise RuntimeError("Invalid channel selection.") 1977 1978 cmin = max(pmin, expr_pmin) 1979 cmax = min(pmax, expr_pmax) 1980 # if the given range of channel selection has overwrap with 1981 # that of current spw, output the overwrap area. 1982 if (cmin <= cmax): 1983 cmin = float(int(cmin + 0.5)) 1984 cmax = float(int(cmax + 0.5)) 1985 crange_list.append([cmin, cmax]) 1986 1987 if (len(crange_list) == 0): 1988 crange_list.append([]) 1989 1990 if res.has_key(spw): 1991 res[spw].extend(crange_list) 1992 else: 1993 res[spw] = crange_list 1994 1995 # restore original values 1996 self.set_restfreqs(orig_restfreq_list) 1997 self.set_freqframe(orig_frame) 1998 self.set_doppler(orig_doppler) 1999 self.set_unit(orig_unit) 2000 2001 return res 2002 2003 @asaplog_post_dec 2004 def get_first_rowno_by_if(self, ifno): 2005 found = False 2006 for irow in xrange(self.nrow()): 2007 if (self.getif(irow) == ifno): 2008 res = irow 2009 found = True 2010 break 2011 2012 if not found: raise RuntimeError("Invalid IF value.") 2013 2014 return res 2015 2016 @asaplog_post_dec 2017 def _get_coordinate_list(self): 2018 res = [] 2019 spws = self.getifnos() 2020 for spw in spws: 2021 elem = {} 2022 elem['if'] = spw 2023 elem['coord'] = self.get_coordinate(self.get_first_rowno_by_if(spw)) 2024 res.append(elem) 2025 2026 return res 2027 2028 ################################## 2029 1625 2030 @asaplog_post_dec 1626 2031 def parse_maskexpr(self, maskstring):
Note:
See TracChangeset
for help on using the changeset viewer.