-
Notifications
You must be signed in to change notification settings - Fork 9
/
ms2util.py
1344 lines (1191 loc) · 55.1 KB
/
ms2util.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
#
# A collection of utilities for ms2 (lookup tables etc..)
# Rewrite in Python, January 2013, Harro Verkouter
#
#
# $Id: ms2util.py,v 1.20 2017-01-27 13:50:28 jive_cc Exp $
#
# $Log: ms2util.py,v $
# Revision 1.20 2017-01-27 13:50:28 jive_cc
# HV: * jplotter.py: small edits
# - "not refresh(e)" => "refresh(e); if not e.plots ..."
# - "e.rawplots.XXX" i.s.o. "e.plots.XXX"
# * relatively big overhaul: in order to force (old) pyrap to
# re-read tables from disk all table objects must call ".close()"
# when they're done.
# Implemented by patching the pyrap.tables.table object on the fly
# with '__enter__' and '__exit__' methods (see "ms2util.opentable(...)")
# such that all table access can be done in a "with ..." block:
# with ms2util.opentable(...) as tbl:
# tbl.getcol('DATA') # ...
# and then when the with-block is left, tbl gets automagically closed
#
# Revision 1.19 2015-09-23 12:28:36 jive_cc
# HV: * Lorant S. requested sensible requests (ones that were already in the
# back of my mind too):
# - option to specify the data column
# - option to not reorder the spectral windows
# Both options are now supported by the code and are triggered by
# passing options to the "ms" command
#
# Revision 1.18 2015-04-29 14:34:14 jive_cc
# HV: * add support for retrieving the actual frequencies of the channels
#
# Revision 1.17 2015-02-16 12:51:07 jive_cc
# HV: * "getcolslice()" in combination with j2ms2 inefficient tiling scheme
# results in very poor performance. So now we do the "slicing"
# ourselves in Python [might change again when j2ms2 does things more
# efficiently or when the casa
#
# Revision 1.16 2014-04-25 12:22:30 jive_cc
# HV: * deal with lag data correctly (number of channels/number of lags)
# * there was a big problem in the polarization labelling which is now fixed
#
# Revision 1.15 2014-04-24 20:09:19 jive_cc
# HV: * indexr now uses 'SCAN_NUMBER' column for scan determination
#
# Revision 1.14 2014-04-14 22:08:01 jive_cc
# HV: * add support for accessing scan properties in time selection
#
# Revision 1.13 2014-04-14 14:46:05 jive_cc
# HV: * Uses pycasa.so for table data access waiting for pyrap to be fixed
# * added "indexr" + scan-based selection option
#
# Revision 1.12 2014-04-10 21:14:40 jive_cc
# HV: * I fell for the age-old Python trick where a default argument is
# initialized statically - all data sets were integrating into the
# the same arrays! Nice!
# * Fixed other efficiency measures: with time averaging data already
# IS in numarray so no conversion needs to be done
# * more improvements
#
# Revision 1.11 2014-04-08 23:34:12 jive_cc
# HV: * Minor fixes - should be better now
#
# Revision 1.10 2014-04-08 22:41:11 jive_cc
# HV: Finally! This might be release 0.1!
# * python based plot iteration now has tolerable speed
# (need to test on 8M row MS though)
# * added quite a few plot types, simplified plotters
# (plotiterators need a round of moving common functionality
# into base class)
# * added generic X/Y plotter
#
# Revision 1.9 2013-12-12 14:10:16 jive_cc
# HV: * another savegame. Now going with pythonic based plotiterator,
# built around ms2util.reducems
#
# Revision 1.8 2013-08-20 18:23:50 jive_cc
# HV: * Another savegame
# Got plotting to work! We have stuff on Screen (tm)!
# Including the fance standardplot labels.
# Only three plot types supported yet but the bulk of the work should
# have been done I think. Then again, there's still a ton of work
# to do. But good progress!
#
# Revision 1.7 2013-06-19 12:28:44 jive_cc
# HV: * making another savegame
#
# Revision 1.6 2013-03-31 17:17:56 jive_cc
# HV: * another savegame
#
# Revision 1.5 2013-03-09 16:59:07 jive_cc
# HV: * another savegame
#
# Revision 1.4 2013-02-19 16:53:29 jive_cc
# HV: * About time to commit - make sure all edits are safeguarded.
# Making good progress. baselineselection, sourceselection and
# timeselection working
#
# Revision 1.3 2013-02-11 09:40:33 jive_cc
# HV: * saving work done so far
# - almost all mapping-functionality is in place
# - some of the API functions are starting to appear
#
# Revision 1.2 2013-01-29 12:23:45 jive_cc
# HV: * time to commit - added some more basic stuff
#
# Revision 1.1.1.1 2013-01-25 09:53:40 jive_cc
# HV: Initial import of some of the python code for the
# rewrite of standardplots into python
#
# Revision 1.3 2003-10-29 12:35:38 verkout
# HV: Changed the way of how to open sub-tables of a MS; now using the right way rather than the JIVE way..:)
#
# Revision 1.2 2001/04/19 14:59:32 verkout
# HV: Added project detection+added it to the plots properties
#
# Revision 1.1.1.1 2001/04/06 13:34:34 verkout
# Files, new + from jivegui/MS1
from __future__ import print_function
from functional import map_, filter_, zip_, enumerate_
from six import iteritems
from functools import reduce
import itertools, operator, math, re, jenums, hvutil, numpy, pyrap.tables, pyrap.quanta, sys
import datetime, copy, collections
# Return a nicely human readable string representation
# of a MS::TIME column value
as_time = lambda x : pyrap.quanta.quantity(x, "s").formatted("dmy", 9)
# string -> time conversion with a knack
# The system supports dayoffsets. In order for this to work correctly
# it is advisory to fill in the refdata parameter with the zero-date.
# It does not matter what the actual time of day of that refdate will be
# since the date will be automatically truncated to 0h00m on the refdate
def str2time(s, refdate=0.0):
totime = lambda x : pyrap.quanta.quantity(x).totime()
# compute 0h00 on refdate
secday = 60.0*60.0*24.0
time0 = (math.floor(refdate/secday))*secday
# remove all whitespace
s = re.sub(r"\s+", "", s)
# If it was all numeric, assume it's MJD seconds - the internal
# measurement set format
if re.match(r"^\d+(\.\d*)?$", s):
return pyrap.quanta.quantity(float(s), "s")
# rxrelday = RE to see if it was a relative day + timestamp
# no more slashes may occur after the first one
relday = re.match(r"^([+-]?\d+)/([^/]+)$", s)
# IF it was a relative day - take out the relative day nr
# and the actual timestamp
if relday:
time0 += (secday * float(relday.group(1)))
s = relday.group(2)
# now we should be left with a timestamp which casa(core) can
# parse as a timestring
# There is one form of timestring that casacore does not recognize:
# MMmSS.SSSs (just minutes and seconds)
# Let's add our own format for that
minsec = re.match(r"^(\d+)m(\d+(\.\d*)?s)?$", s)
tmq = None
if minsec:
# set the timequantity
tmq = totime(minsec.group(1)+"min")
if minsec.group(2):
tmq += totime(minsec.group(2))
else:
# attempt the conversion
tmq = totime(s)
# if we have < 1 day's worth of seconds, add the refdate
if tmq.get_value("s")<secday:
tmq += pyrap.quanta.quantity(time0, "s").totime()
return tmq
# we like our table opens to be silent
# so if the user did not pass a value for
# 'ack' (the 5th argument for the table
# function) we'll add a default of 'False'
def opentable(*args, **kwargs):
d = {}
if len(args)<5:
d['ack'] = False
d.update(kwargs)
# HV: Before we used to return the pyrap.tables.table() return value unadorned.
# This, however, does not have success written all over it. Specifically,
# either the pyrap, boost.python or C++ Table::... code holds a reference
# from '<name>' to internal object.
# This means that if '<name>' changes on disk, pyrap.tables.table('<name>')
# will keep on returning the _old_ (!!!!) table information AS LONG AS THERE
# ARE pyrap.tables.table() OBJECTS ALIVE THAT HAVE NOT __EXPLICITLY__ (!!!!)
# CALLED ".close()" ON THE TABLE.
#
# This seems to happen based on *name* of the table because the table objects
# you get returned from pyrap.tables.table() are /different/; you really
# get different instances of table object:
# >>> t1 = table('foo.ms')
# >>> t2 = table('foo.ms')
# >>> id(t1) == id(t2)
# False
#
# We fix this by adding '__enter__' and '__exit__' "methods" to the pyrap.tables.table
# object such that it can be used as a context manager; i.e.:
#
# with opentable('<name>') as tbl:
# # blah
# print len(tbl)
#
# When 'tbl' goes out of scope (leaves the with-scope), the table is automagically closed.
realtable = pyrap.tables.table(*args, **d)
# casacore-python (or python-casacore) - the successor of pyrap,
# does seem to have '__enter__'/'__exit__' methods
if not hasattr(realtable, '__enter__'):
setattr(realtable, '__enter__', lambda *args: realtable)
if not hasattr(realtable, '__exit__'):
setattr(realtable, '__exit__' , lambda *args: realtable.close())
return realtable
# Exceptions that can be thrown
class NotAMeasurementSet(Exception):
def __init__(self, ms):
self.ms = ms
def __str__(self):
return "{0} is not a minimal MeasurementSet!".format(self.ms)
class InvalidFreqGroup(Exception):
def __init__(self,fgrp):
self.freqGroup = fgrp
def __str__(self):
return "FREQ_GROUP {0} not found!".format(self.freqGroup)
class InvalidSubband(Exception):
def __init__(self,sb):
self.subbandIdx = sb
def __str__(self):
return "Subband #{0} is out of range".format(self.subbandIdx)
class InvalidSpectralWindow(Exception):
def __init__(self, spw):
self.spWinId = spw
def __str__(self):
return "SPECTRAL_WINDOW_ID {0} is unknown".format(self.spWinId)
class InvalidDataDescId(Exception):
def __init__(self, ddid):
self.dataDescId = ddid
def __str__(self):
return "DATA_DESC_ID {0} is unknown".format(self.dataDescId)
class InvalidPolIdForSubband(Exception):
def __init__(self,pid,*args):
self.polarizationId = pid
self.moreArgs = args
def __str__(self):
return "The polarizationId {0} is not valid for {1}".format(self.polarizationId, self.moreArgs)
class InvalidAntenna(Exception):
def __init__(self, id):
self.antennaId = id
def __str__(self):
return "The antenna ID {0} is unknown".format(self.antennaId)
class InvalidBaselineId(Exception):
def __init__(self, blid):
self.baselineIdx = blid
def __str__(self):
return "The baseline id '{0}' (number or name) is unknown".format(self.baselineIdx)
class InvalidPolarizationCode(Exception):
def __init__(self,polstr):
self.polarizationCode = polstr
def __str__(self):
return "The polarization code {0} is not a recognized polarization".format(self.polarizationCode)
class InvalidFieldId(Exception):
def __init__(self,fld):
self.fieldId = fld
def __str__(self):
return "The field code {0} is not a recognized field".format(self.fieldId)
##############################################################################################
## data structures
##############################################################################################
##
## The polarization mapping
##
class polarizationmap():
## there exists a static mapping between AIPS++ polarization code
## numerical value (as found in the MS) to the human readable equivalent
## so we start with a bit of class-level code to deal with that
class polcode():
def __init__(self, **kwargs):
self.code = kwargs['id']
self.name = kwargs['name']
CorrelationStrings = [
polcode(id=1 ,name='I' ),
polcode(id=2 ,name='Q' ),
polcode(id=3 ,name='U' ),
polcode(id=4 ,name='V' ),
polcode(id=5 ,name='RR'),
polcode(id=6 ,name='RL'),
polcode(id=7 ,name='LR'),
polcode(id=8 ,name='LL'),
polcode(id=9 ,name='XX'),
polcode(id=10,name='XY'),
polcode(id=11,name='YX'),
polcode(id=12,name='YY')#,
#polcode(id=13,name='RX'),
#polcode(id=14,name='RY'),
#polcode(id=15,name='LX'),
#polcode(id=16,name='LY'),
#polcode(id=17,name='XR'),
#polcode(id=18,name='XL'),
#polcode(id=19,name='YR'),
#polcode(id=20,name='YL')
]
@staticmethod
def polcodeById(corrtype):
try:
[pcode] = filter_(lambda x: x.code==corrtype, polarizationmap.CorrelationStrings)
return copy.deepcopy(pcode)
except ValueError:
raise InvalidPolarizationCode(corrtype)
# Return a list of correlationIDs for which the name matches the regex...
@staticmethod
def matchingName2ID(rx):
return [x.code for x in polarizationmap.CorrelationStrings if rx.match(x.name)]
@staticmethod
def correlationId2String(id):
return ",".join( polarizationmap.correlationId2Strings(id) )
@staticmethod
def correlationId2Strings(id):
# Assume 'id' is a list-of-correlation-types
return map_(lambda x: polarizationmap.polcodeById(x).name, id)
# Turn an array of polarizationcombination strings into
# an array of correlation types
# FIXME XXX FIXME
@staticmethod
def string2CorrelationId(s):
lows = s.lower()
cs = map_(lambda x: polarizationmap.polcode(id=x.code, name=x.name.lower()), polarizationmap.CorrelationStrings)
def findid(nm):
try:
[x] = filter_(lambda y: y.name.lower()==nm.lower(), cs)
return x.code
except ValueError:
raise InvalidPolarizationCode(nm)
(result,names) = ([],s.lstrip().rstrip().split(','))
map(result.append, map(findid, names))
return result
# This is where the instance methods go - these are
# for mapping what was actually contained in the MS
# polmap = [(polid, poltypes)]
def __init__(self, polmap):
self.polarizationMap = polmap
# Note: technically speaking the polarization id is
# a direct index into the list of polarizations
# so *technically* we could do self.polarizationMap[id]
# but then we - because of pythonicity - also can
# address polarization id "-1", which would be the
# last one in the list. Which we don't want because
# then the "id" isn't actually the id anymore.
# return the polarizations as a string with comma-separated
# polarization strings
def getPolarizationsString(self, id):
return ",".join(self.getPolarizations(id))
# return the polarization strings as list of strings
def getPolarizations(self, id):
try:
# enforce unique search result otherwise bomb out
[x] = filter_(lambda y: y[0]==id, self.polarizationMap)
return polarizationmap.correlationId2Strings(x[1])
except ValueError:
raise InvalidPolarizationCode(id)
def getCorrelationTypes(self, id):
try:
# enforce unique search result otherwise bomb out
[x] = filter_(lambda y: y[0]==id, self.polarizationMap)
return x[1]
except ValueError:
raise InvalidPolarizationCode(id)
def polarizationIDs(self):
return [x[0] for x in self.polarizationMap]
def nrPolarizationIDs(self):
return len(self.polarizationMap)
def makePolarizationMap(nm, **args):
errf = hvutil.mkerrf("makePolarizationMap({0})".format(nm))
with opentable(nm) as tbl:
with opentable(tbl.getkeyword('POLARIZATION')) as poltab:
if len(poltab)==0:
return errf("No rows in POLARIZATION table?!")
# read all polarization types
get_type = lambda x: x['CORR_TYPE']
return polarizationmap(zip_(itertools.count(), map(get_type, poltab)))
ct2str = polarizationmap.correlationId2String
class subband:
# ddmap = Map POLARIZATION_ID => DATA_DESCRIPTION_ID,
# all polarization IDs this subband was correlated with
#
# TODO FIXME
# it is possible that >1 DATA_DESCRIPTION_ID maps to
# the same (SPECTRAL_WINDOW, POLID) tuple. So basically
# the ddmap should be
# Map POLARIZATION_ID => [DATA_DESCRIPTION_ID, ...]
def __init__(self, spwid, f0, nchan, bw, ddmap, **kwargs):
self.spWinId = spwid
self.frequency = f0
self.numChannels = nchan
self.bandWidth = bw
self.datadescMap = ddmap
self.formatStr = "{frequency[0]:5.4f}MHz/{bw:<.1f}MHz nch={nch} {polarizations}"
self.domain = kwargs.get('domain', jenums.Type.Unknown)
# in lag domain the frequency axis becomes delay in units of s
if self.domain==jenums.Type.Lag:
freqs = self.frequency
self.frequency = numpy.arange(-self.numChannels, self.numChannels) * 1./(2*self.bandWidth)
self.numChannels = len(self.frequency)
self.formatStr = "{0:5.4f}MHz {1:.2}us - {2:.2}us".format(freqs[0]/1.e6, self.frequency[0]/1e-6, self.frequency[-1]/1e-6) + " nlag={nch} {polarizations}"
def __str__(self, **kwargs):
pmap = kwargs.get('polarizationMap', None)
if pmap is None:
polarizations = " ".join(map("P{0}".format, self.datadescMap.keys()))
else:
polarizations = " ".join(map(lambda x: "P{0}={1}".format(x, ct2str(pmap.getCorrelationTypes(x))), self.datadescMap.keys()))
return self.formatStr.format(frequency=(self.frequency/1.0e6), bw=self.bandWidth/1.e6, spw=self.spWinId, ddMap=self.datadescMap, nch=self.numChannels, polarizations=polarizations)
# (attempt to) map the indicated polarization id to a DATA_DESCRIPTION_ID
# None if we don't have this polid
# FIXME TODO
# would become a list-of-data_desc_id's
def mapPolId(self, polid):
try:
return self.datadescMap[polid]
except KeyError:
raise InvalidPolIdForSubband(polid, str(self))
# unmap a DATA_DESC_ID to a POLARIZATION_ID, if we have the DATA_DESC_ID,
# [] otherwise
def unmapDDId(self, ddid):
def foldfn(x, acc):
# Note: 'x' is always a (int, int) tuple but by analyzing it like below
# we can perform the extending of the accumulator only if the
# datadescription id we're looking for is in this subband's
# data description map in a one liner
###
### FIXME TODO
### the data desc map would be a list-of-data_desc_id's so we must
### look for 'ddid' in the list of datadescids!
### (use "if ddid in v")
acc.extend([k for (k,v) in [x] if v==ddid])
return acc
return hvutil.dictfold(foldfn, [], self.datadescMap)
def getDataDescIDs(self):
return list(self.datadescMap.values())
## the spectral map class
## holds a mapping and defines lookup functions on those mappings
class spectralmap:
# spm = Map FQGROUP -> [subband()]
# but we want each subband tagged with a numberical index;
# our spectralmap should look like:
# spectralMap = Map FQGROUP -> [(sbidx, subband())]
# thus we tag the incoming list of subbands with a
# logical subband number, starting from 0
#
# Update: assume FQGROUP == (FREQ_GROUP :: int, FREQ_GROUP_NAME :: string)
def __init__(self, spm):
self.spectralMap = hvutil.dictmap(lambda k_v: enumerate_(k_v[1]), spm)
# Simple API functions
def nFreqId(self):
return len(self.spectralMap)
def freqIds(self):
return map_(operator.itemgetter(0), self.spectralMap.keys())
def freqGroupName(self, fq):
try:
[key] = filter_(lambda num_name: fq==num_name[0] or fq==num_name[1], self.spectralMap.keys())
return key[1]
except ValueError:
raise InvalidFreqGroup(fq)
# our lookup key may now be either int or string.
# may raise KeyError if not found
def _findFQ(self, fqkey):
try:
[key] = filter_(lambda num_name: fqkey==num_name[0] or fqkey==num_name[1], self.spectralMap.keys())
return self.spectralMap[key]
except ValueError:
raise KeyError
def subbandsOfFREQ(self, fq):
try:
return self._findFQ(fq)
except KeyError:
raise InvalidFreqGroup(fq)
# FREQ0 indexed by spectral-window id and via "FREQGRP/SUBBAND"
# typically we first unmap the SPECTRAL_WINDOW_ID to a
# FREQGROUP/SUBBAND and then do the (easy) lookup
def frequencyOfSPW(self, spwid):
try:
rv = self.unmap(spwid)
return self.frequenciesOfFREQ_SB(rv.FREQID, rv.SUBBAND)[0]
except AttributeError:
raise InvalidSpectralWindow(spwid)
def frequenciesOfSPW(self, spwid):
try:
rv = self.unmap(spwid)
return self.frequenciesOfFREQ_SB(rv.FREQID, rv.SUBBAND)
except AttributeError:
raise InvalidSpectralWindow(spwid)
def frequencyOfFREQ_SB(self, fq, sb):
try:
fqref = self._findFQ(fq)
return fqref[sb][1].frequency[0]
except KeyError:
raise InvalidFreqGroup(fq)
except IndexError:
raise InvalidSubband(sb)
def frequenciesOfFREQ_SB(self, fq, sb):
try:
fqref = self._findFQ(fq)
return fqref[sb][1].frequency
except KeyError:
raise InvalidFreqGroup(fq)
except IndexError:
raise InvalidSubband(sb)
# Id. for NUMCHAN
def numchanOfSPW(self, spwid):
try:
rv = self.unmap(spwid)
return self.numchanOfFREQ_SB(rv.FREQID, rv.SUBBAND)
except AttributeError:
raise InvalidSpectralWindow(spwid)
def numchanOfFREQ_SB(self, fq, sb):
try:
fqref = self._findFQ(fq)
return fqref[sb][1].numChannels
except KeyError:
raise InvalidFreqGroup(fq)
except IndexError:
raise InvalidSubband(sb)
# Id. for the polarization IDs - NOTE! These are lists!
def polarizationIdsOfSPW(self, spwid):
try:
rv = self.unmap(spwid)
return self.polarizationIdsOfFREQ_SB(rv.FREQID, rv.SUBBAND)
except AttributeError:
raise InvalidSpectralWindow(spwid)
def polarizationIdsOfFREQ_SB(self, fq, sb):
try:
# the polarization Ids are the keys in the datadescription map,
# the values are the associated data description id's
fqref = self._findFQ(fq)
return list(fqref[sb][1].datadescMap.keys())
except KeyError:
raise InvalidFreqGroup(fq)
except IndexError:
raise InvalidSubband(sb)
# Id. for the bandwidth
def bandwidtOfSPW(self, spwid):
try:
rv = self.unmap(spwid)
return self.bandwidthOfFREQ_SB(rv.FREQID, rv.SUBBAND)
except AttributeError:
raise InvalidSpectralWindow(spwid)
def bandwidthOfFREQ_SB(self, fq, sb):
try:
fqref = self._findFQ(fq)
return fqref[sb][1].bandWidth
except KeyError:
raise InvalidFreqGroup(fq)
except IndexError:
raise InvalidSubband(sb)
def typeOfFREQ_SB(self, fq, sb):
try:
fqref = self._findFQ(fq)
return fqref[sb][1].numStr
except KeyError:
raise InvalidFreqGroup(fq)
except IndexError:
raise InvalidSubband(sb)
# Id. for the DATA_DESCRIPTION_IDs - these can be used in TaQL directly
# we go from spwid => (FQGROUP, SUBBAND), then look inside that SUBBAND
# for POLID
def datadescriptionIdOfSPW(self, spwid, polid):
try:
rv = self.unmap(spwid)
return self.datadescriptionIdOfFREQ_SB_POL(rv.FREQID, rv.SUBBAND, polid)
except AttributeError:
raise InvalidSpectralWindow(spwid)
def datadescriptionIdOfFREQ_SB_POL(self, fq, sb, polid):
try:
# the polarization Ids are the keys in the datadescription map,
# the values are the associated data description id's
fqref = self._findFQ(fq)
return fqref[sb][1].mapPolId(polid)
except KeyError:
raise InvalidFreqGroup(fq)
except IndexError:
raise InvalidSubband(sb)
# return all the spectral window Id's for the given FREQGROUP
def spectralWindows(self, fq):
try:
fqref = self._findFQ(fq)
return [x[1].spWinId for x in fqref]
except KeyError:
raise InvalidFreqGroup(fq)
# How many subbands has the indicated FREQGROUP?
def nSubbandOfFREQ(self, fq):
try:
return len(self._findFQ(fq))
except KeyError:
raise InvalidFreqGroup(fq)
# Return all DATA_DESCRIPTION_IDs (mostly used to be able
# to see if all datadescription ids are selected: if that's true,
# the TaQL can be dropped...)
def datadescriptionIDs(self):
# All the datadescription IDs are contained in the subband objects,
# which are contained in lists-per-freqgroup
# so we 'fold' over all our freq-groups and collect all data-desc-ids
def foldfn(x, acc):
(fgrp, sblist) = x
for (idx,sb) in sblist:
acc.extend(sb.getDataDescIDs())
return acc
return sorted(list(hvutil.dictfold(foldfn, [], self.spectralMap)))
#return sorted(list(set(dictfold(foldfn, [], self.spectralMap))))
# Map FREQID (integer) + zero-based subbandnr into
# a spectral window id
#
# return None on error, SEPCTRAL_WINDOW_ID (int) otherwise
def map(self, freqid, sbid):
try:
if sbid<0:
raise InvalidSubband(sbid)
fqref = self._findFQ(freqid)
return fqref[sbid][1].spWinId
except KeyError:
raise InvalidFreqGroup(freqid)
except IndexError:
raise InvalidSubband(sbid)
# Unmap a given spectral window id (note: this is the *zero*based rownumber!!)
#
# retval = None on error, record with FREQID=xxx and SUBBAND=yyyy where
# FREQID is the freqid (you'd never guessed that eh?) and SUBBAND is the *zero*
# based subband nr in this freqid that the requested spectral window represents..
def unmap(self, spwid):
for (fq, sbs) in iteritems(self.spectralMap):
for (idx, sb) in sbs:
if sb.spWinId==spwid:
o = type('',(),{})()
o.FREQID = fq[0]
o.SUBBAND = idx
return o
raise InvalidSpectralWindow(spwid)
# Unmap a DATA_DESC_ID into FREQID/SUBBAND/POLID
def unmapDDId(self, ddid):
# look in all FREQGROUPS, in all SUBBANDS for the given DATA_DESC_ID
for (k,v) in iteritems(hvutil.dictmap(lambda k_v: [sb for sb in k_v[1] if sb[1].unmapDDId(ddid)], self.spectralMap)):
if len(v)>1:
raise RuntimeError("Non-unique search result for DATA_DESC_ID={0}".format(ddid))
if len(v)==1:
class sbres:
def __init__(self,fid,sb,pid):
self.FREQID = fid
self.SUBBAND = sb
self.POLID = pid
def fsbpol(self):
return (self.FREQID, self.SUBBAND, self.POLID)
[pol] = v[0][1].unmapDDId(ddid)
return sbres(k[0], v[0][0], pol)
else:
raise InvalidDataDescId(ddid)
return None
# Print ourselves in readable format
def __str__(self):
r = "*** SPWIN <-> FREQID/SUBBAND MAP ***\n"
for (fgrp,subbands) in iteritems(self.spectralMap):
r = r+"FQ={0} ({1})\n".format(fgrp[0], fgrp[1])
for (idx,sb) in subbands:
r = r + " {0:2d}: {1}".format(idx,sb)
if len(sb.datadescMap)>0:
r = r + " // {0}\n".format( \
",".join(["POLID #{0} (DDI={1})".format(pol,ddi) \
for (pol,ddi) in iteritems(sb.datadescMap)]) \
)
return r
def assertMinimalMS(nm):
ti = pyrap.tables.tableinfo(nm)
if ti['type']!='Measurement Set':
raise NotAMeasurementSet(nm)
def makeSpectralMap(nm, **kwargs):
errf = hvutil.mkerrf("makeSpectralMap({0})".format(nm))
with opentable(nm) as tbl:
datadom = getDataDomain(nm, **kwargs)
# verify that we have a non-empty spectral window table
with opentable( tbl.getkeyword('SPECTRAL_WINDOW') ) as spwint:
if len(spwint)==0:
return errf("No rows in spectral window table")
if 'FREQ_GROUP' not in spwint.colnames():
return errf("No 'FREQ_GROUP' column found?")
# id. for the datadescription table
with opentable(tbl.getkeyword('DATA_DESCRIPTION')) as ddt:
if len(ddt)==0:
return errf("No rows in data description table")
# spmap:
# Map [FREQ_GROUP] -> { Map [SubbandIdx] = > {FREQ0, NCHAN, DATADESCMAP} }
spmap = collections.defaultdict(list)
# columns in the spectral window table
freqgroup = lambda x : x['FREQ_GROUP']
fgrpname = lambda x : x['FREQ_GROUP_NAME']
chanfreq = lambda x : x['CHAN_FREQ']
numchan = lambda x : x['NUM_CHAN']
totbw = lambda x : x['TOTAL_BANDWIDTH']
# columns in the datadescription table
spwid = lambda x : x['SPECTRAL_WINDOW_ID']
polid = lambda x : x['POLARIZATION_ID']
# Do different things depending on wether to load the
# full spectral window table or only the ones that
# are actually used in the MS
# function to generate the spectral window's "key" function
# from a row in the SPECTRAL_WINDOW table
key_f = lambda row: (freqgroup(row), fgrpname(row))
# function to create a fully fledged subband() object
# from a row in the SPECTRAL_WINDOW table + some extra
# info (the row# in the SPECTRAL_WINDOW table and
# the 'polarization-id' => 'data_description_id' mapping)
sb_f = lambda n, row, pmap: subband(n, chanfreq(row), numchan(row), totbw(row), pmap, domain=datadom.domain)
if not kwargs.get('unique', False):
## Read the spectral window table and find all possible
## correlations with it
for (idx,row) in enumerate(spwint):
# each spectral window may appear multiple times in the datadescription
# table - with different polarization IDs. So we build a mapping
# of polid -> datadescid for all the datadescriptions that match the
# current spectral window id
ddmap = dict((polid(row), ddid) for (ddid,row) in enumerate(ddt) if spwid(row)==idx)
# build an entry for this subband and add to current mapping
spmap[ key_f(row) ].append( sb_f(idx, row, ddmap) )
else:
## Input: list-of-DATA_DESC_IDs
## unmap to (key, subband) pairs
## group by (key)
## join subbands [could be same subband but different polarization id]
## Extract the unique DATA_DESC_IDs from the main table, then we
## simply unmap those
def reductor(acc, ddid):
# get the spectral-window id and polid for this DDID
(spw, pol) = (spwid(ddt[ddid]), polid(ddt[ddid]))
# create the subband
row = spwint[spw]
sb = sb_f(spw, row, {pol:ddid})
key = key_f(row)
try:
# check if this sb is already present in the current
# freqgroup, we may have hit another POLID
processed = False
for k in spmap[key]:
if k.spWinId==sb.spWinId:
# irrespective of the outcome of this operation we can
# flag this particular data desc id as processed, such that
# it don't get added
processed = True
# now update the subband object's data desc map with this one's
# let's verify that existing (POLID -> DATA_DESC_ID) are not violated!
for (p,d) in iteritems(sb.datadescMap):
if p in k.datadescMap:
if k.datadescMap[p]!=d:
#### FIXME TODO
#### if we update the subband() object to support >1 DATA_DESC_ID
#### mapping to the same (SPECTRAL_WINDOW, POLID) tuple this can go
raise RuntimeError("Inconsistent DATA_DESCRIPTION/SPECTRAL_WINDOW table. \
Spectral window {0} appears with POLID={1} -> DATA_DESC_ID={2} \
but also as POLID{3} -> DATA_DESC_ID{4}".format( \
sb.spWinId, p, d, p, k.datadescMap[p] ))
else:
# ok, no entry for this polarization yet
k.datadescMap[p] = d
# if we didn't find the subband yet, add it to the current key!
if not processed:
spmap[key] = spmap[key]+[sb]
except KeyError:
spmap[key] = [sb]
return spmap
# On an 8.7Mrow (that's 8.7e6 rows) MS this statement takes 54.96(!) seconds
# Un-effing-believable!
# pyrap.tables.taql("select unique DATA_DESC_ID from $tbl", locals={"tbl":tbl}).getcol("DATA_DESC_ID"), \
#
# Let's see if we can do better!
# Thus:
#
# pyrap.tables.taql("select unique DATA_DESC_ID ...") takes 54.96s (8.7Mrow)
# set(sorted(tbl.getcol('DATA_DESC_ID'))) takes 4.38s ( .. )
# set(tbl.getcol('DATA_DESC_ID')) takes 1.35s ( .. )
# tbl.getcol('DATA_DESC_ID') takes 0.36s ( .. )
# numpy.unique( tbl.getcol('DATA_DESC_ID') ) takes 0.64s ( .. )
#
# Looks like we have a winner!
spmap = reduce(reductor, numpy.unique( tbl.getcol('DATA_DESC_ID') ), {})
# do not forget to sort all subbands by frequency
sort_order = kwargs.get('spw_order', 'by_frequency').lower()
if sort_order=='by_frequency':
sortfn = lambda x: x.frequency[0] #lambda x,y: cmp(x.frequency[0], y.frequency[0])
elif sort_order=='by_id':
sortfn = lambda x: x.spWinId #lambda x,y: cmp(x.spWinId, y.spWinId)
else:
raise RuntimeError("The spectral window ordering function {0} is unknown".format( kwargs.get('spw_order') ))
return spectralmap( hvutil.dictmap( lambda kvpair : sorted(kvpair[1], key=sortfn), spmap) )
##
## The baseline mapping
##
class baselinemap:
# antennalist is [(antenna, id), ...]
# baselineleList = [(baselineId, baselineName),...]
def __init__(self, antennalist, **kwargs):
# keep a sorted list ourselves
self.antennaDict = dict( (al[1], al[0]) for al in antennalist )
# Check if the 'baselines' were explicitly given. If not, we form
# the baselines ourselves out of all antenna pairs
# The baselines can be passed in as keyword-arg:
# ..., baselines=[(x,y), ...], ...
# a list of antenna pairs
bls = kwargs.get('baselines', None)
if bls is None:
# Make local list to make sure we can iterate over it > once (on Py3 it's a view)
aidx = list(enumerate(self.antennaDict.keys()))
bls = [(x, y[1]) for (idx, x) in aidx for y in aidx[idx:]]
# Now transform the list of indices into a list of codes + names
ANAME = self.antennaDict.get
self.baselineList = map_(lambda x_y: (x_y, "{0}{1}".format(ANAME(x_y[0]), ANAME(x_y[1]))), bls)
def baselineNames(self):
return map_(operator.itemgetter(1), self.baselineList)
def baselineIndices(self):
return map_(operator.itemgetter(0), self.baselineList)
def baselineName(self, blidx):
try:
[(x,y)] = filter_(lambda idx_nm: idx_nm[0]==blidx, self.baselineList)
return y
except ValueError:
raise InvalidBaselineId(blidx)
def baselineIndex(self, blname):
try:
[(x,y)] = filter_(lambda idx_nm: idx_nm[1]==blname, self.baselineList)
return x
except ValueError:
raise InvalidBaselineId(blname)
def antennaNames(self):
return list( self.antennaDict.values() )
# return the list of (antenna, id) tuples, sorted by id
def antennas(self):
return sorted(self.antennaDict.items(), key=operator.itemgetter(0))
def antennaName(self, aid):
try:
return self.antennaDict[aid]
except KeyError:
raise InvalidAntenna(aid)
def antennaId(self, name):
try:
namelower = name.lower()
[(nm, aid)] = filter_(lambda ant_antid: ant_antid[1].lower()==namelower, self.antennaDict.items())
return aid
except ValueError:
raise InvalidAntenna(name)
def makeBaselineMap(nm, **kwargs):
errf = hvutil.mkerrf("makeBaselineMap({0})".format(nm))
with opentable(nm) as tbl:
with opentable( tbl.getkeyword('ANTENNA') ) as antab:
if len(antab)==0:
return errf("No rows in the ANTENNA table")
# If we want to know only the antenna's that are
# actually *used* ....
filter_f = None
baselines = None
if 'unique' not in kwargs or not kwargs['unique']:
filter_f = lambda x : True
else:
# compile the list of unique antenna-id's that are
# actually used in the MS
#
# Again - we hit the ***** 'performance' of
# pyrap.tables.taql.
#
# uniqry = lambda col:
# set(
# pyrap.tables.taql("SELECT {0} AS FOO from $tbl".format(col), locals={"tbl":tbl}).getcol("FOO")
# )
#
# All measurements taken on a 8.7Mrow (8.7e6 rows) MS
#
# uniqry( "ANTENNA1" ) takes 1.51s
# tbl.getcol('ANTENNA1') takes 0.36s
# x = tbl.getcol('ANTENNA1')
# numpy.unique( x ) takes 0.23s
# set(numpy.unique(x)) takes 0.24s
#
# Thus for getting the unique antennas it's fastests to
# get the columns for ANTENNA1, ANTENNA2, create sets and union them
#
#
#uniqry = lambda col: set(pyrap.tables.taql("SELECT {0} AS FOO from $tbl".format(col), locals={"tbl":tbl}).getcol("FOO"))
#ants = uniqry("ANTENNA1") | uniqry("ANTENNA2")
a1 = tbl.getcol('ANTENNA1')
a2 = tbl.getcol('ANTENNA2')
ants = set(numpy.unique(a1)) | set(numpy.unique(a2))
filter_f = lambda x : x in ants
# retrieve the uniqe baselines. This also takes a LOOOONG time
# so we gonna do it mighty different.
# since we already have ANTENNA1, ANTENNA2 columns we're going
# to play a neat trick.