Running getRM in python3

[AstroRipples]

Hi folks. Here's hoping some of y'all still monitor this repo...

I've just built RMextract inside a py3 virtualenv, and I'm trying to run getRM to estimate the ionospheric RM for some of my MeerKAT data. The virtualenv has been built okay, but I cannot get anything out of getRM.

I'm trying to run the following :

targ = SkyCoord('13:31:08.290000', '+30:30:33.00000', frame='fk5', equinox='J2000', unit=(u.hourangle, u.degree))
rd = np.array([ targ.ra.to(u.rad).value, targ.dec.to(u.rad).value ]) #3C286  Ra, Dec in radians
a=tab.taql('calc MJD("2021/04/25/19:33:12")')[0]*3600*24
b=tab.taql('calc MJD("2021/04/25/22:29:16")')[0]*3600*24

ms472='msdir/1619378472-cal.MS'
RMdict=grm.getRM(MS=ms472,
                 server="ftp://gnss.oma.be/gnss/products/IONEX/",
                 prefix='ROBR',
                 ionexPath="IONEXdata/",
                 radec=rd,
                 timestep=600,
                 timerange = [a, b],
                 ha_limit=np.pi)

and when I do, the MS is successfully read, and RMextract prints out numbers from 001 to 250 (twice over...):

Successful readonly open of default-locked table 1619378472-cal.MS: 25 columns, 797880 rows
Successful readonly open of default-locked table 1619378472-cal.MS/FIELD: 9 columns, 3 rows
Successful readonly open of default-locked table 1619378472-cal.MS/ANTENNA: 8 columns, 61 rows
Successful readonly open of default-locked table 1619378472-cal.MS/ANTENNA: 8 columns, 61 rows
001
002
003
004
005
006
007
008
...

but then I get a UnicodeDecodeError:

---------------------------------------------------------------------------
UnicodeDecodeError                        Traceback (most recent call last)
<ipython-input-8-7e8e34bf8feb> in <module>
      6                  timestep=600,
      7                  timerange = [a, b],
----> 8                  ha_limit=np.pi)

~/rmextract/lib/python3.6/site-packages/RMextract/getRM.py in getRM(MS, server, prefix, ionexPath, earth_rot, timerange, use_azel, ha_limit, use_filter, **kwargs)
    157         if not use_proxy:
    158             if not "http" in server: #ftp server use ftplib
--> 159                 ionexf=ionex.getIONEXfile(time=date_parms,server=server,prefix=prefix,outpath=ionexPath)
    160             else:
    161                 ionexf=ionex.get_urllib_IONEXfile(time=date_parms,server=server,prefix=prefix,outpath=ionexPath)

~/rmextract/lib/python3.6/site-packages/RMextract/getIONEX.py in getIONEXfile(time, server, prefix, outpath, overwrite)
    692                  overwrite=False):
    693     getIONEXfile.__doc__ = _get_IONEX_file.__doc__
--> 694     return _get_IONEX_file(time, server, prefix, outpath, overwrite)
    695 
    696 def get_TEC_data(times, lonlatpp, server, prefix, outpath, use_filter=None,earth_rot=0.):

~/rmextract/lib/python3.6/site-packages/RMextract/getIONEX.py in _get_IONEX_file(time, server, prefix, outpath, overwrite, backupserver)
    559         filenames += nfilenames
    560         _combine_ionex(outpath, filenames,
--> 561                        prefix + "%03d0.%sI" % (dayofyear, yy))
    562         ftp.quit()
    563         return os.path.join(outpath, prefix + "%03d0.%sI" % (dayofyear, yy))

~/rmextract/lib/python3.6/site-packages/RMextract/getIONEX.py in _combine_ionex(outpath, filenames, newfilename)
    340     firstfile = open(filenames[0])
    341     lastfile = open(filenames[-1])
--> 342     for line in lastfile:
    343         if "EPOCH OF LAST MAP" in line:
    344             lastepoch = line

~/rmextract/lib64/python3.6/encodings/ascii.py in decode(self, input, final)
     24 class IncrementalDecoder(codecs.IncrementalDecoder):
     25     def decode(self, input, final=False):
---> 26         return codecs.ascii_decode(input, self.errors)[0]
     27 
     28 class StreamWriter(Codec,codecs.StreamWriter):

UnicodeDecodeError: 'ascii' codec can't decode byte 0xc2 in position 1589: ordinal not in range(128)

Then if I try and re-run the same task again, it looks like it's not parsing the files correctly as I get an UnboundLocalError:

...
~/rmextract/lib/python3.6/site-packages/RMextract/getIONEX.py in _read_ionex_header(filep)
     59             ntimes = int(stripped.split()[0])
     60 
---> 61     lonarray = np.arange(start_lon, end_lon + step_lon, step_lon)
     62     latarray = np.arange(start_lat, end_lat + step_lat, step_lat)
     63     dtime = endtime - starttime

UnboundLocalError: local variable 'start_lon' referenced before assignment

So what's the deal here? Is something screwed with the format of the IONEX files downloaded from this particular server (ftp://gnss.oma.be/gnss/products/IONEX/)? I tried with the other examples that I've used before (CDDIS, AIUB) but both of those gave me server timeouts...

Suggestions would be much appreciated!

[gmloose]

UnicodeDecodeError errors are a real PITA. There may be a bug somewhere in the code, but a quick possible workaround would be to do the following:

export LC_ALL=C.UTF-8

[AstroRipples]

Okay, this might be a real basic question, but where should I try that? As you'd expect, that syntax doesn't play well with python.

[maaijke]

Hi as for the time outs, did you try

ionex_server = "ftp://cddis.nasa.gov/gnss/products/ionex/" or ionex_server = "ftp://gssc.esa.int/gnss/products/ionex/" or ftp://ftp.aiub.unibe.ch/CODE/ (CODG only)

I will look into the two other problems. It seems there is a corruption in the ROBR file, they are anyway not the best to use (if you want high time resolution go for uqrg).

[maaijke]

So it is indeed an issue with the IONEX files. For ROBR we need to combine several files of one day into one, to get the standard IONEX format. The first time that fails, but it still creates the combined IONEX files (probably empty) The second time it tries to read from the empty file, thus causing the (not so informative error on uninitialized start_lon) error.

I'll check if the general format of ROBR has changed and if the code needs updating, but otherwise I strongly suggest to find a different server prefix combination that woks for you

[AstroRipples]

Hi Maaijke!

Thanks for the details! Okay, now I'm not getting server timeouts -- must have been some slight formatting differences (I'd ripped the server links from an older script). Those updates links now work okay.

However, it seems that getRM is now unable to find the leap second table... I'm not crash-hot at working with virtualenvs, but as far as I can tell I set up the virtualenv correctly (I ran pip3 install python-casacore before running pip3 install rmextract for example) so as far as I can tell it should work, but I get the following:

2021-09-07 11:31:14 WARN    MeasIERS::findTab (file /code/measures/Measures/MeasIERS.cc, line 387)  Requested data table TAI_UTC cannot be found in the searched directories:
2021-09-07 11:31:14 WARN    MeasIERS::findTab (file /code/measures/Measures/MeasIERS.cc, line 387)+ /usr/share/casacore/data/ephemerides/
2021-09-07 11:31:14 WARN    MeasIERS::findTab (file /code/measures/Measures/MeasIERS.cc, line 387)+ /usr/share/casacore/data/geodetic/
2021-09-07 11:31:14 SEVERE  MeasTable::dUTC(Double) (file /code/measures/Measures/MeasTable.cc, line 4281)  Cannot read leap second table TAI_UTC

Do I need to run some additional command(s) to get rmextract to play nice with casacore?

[maaijke]

This is really a casacore warning, leapseconds are of no importance to RMextract at all. To get it corrected I think you need to update your measure tables or the path to it, but as said, not important.

I did find a bug when you have http in the server address, maybe this was what you ran into. Will be solved soon

[AstroRipples]

Thanks Maijke... Okay, hm, if not finding the leap second table is a warning rather than an error, that's okay. However, I still get a RunTimeError when trying to run getRM:

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
<ipython-input-2-1ab25ef87fdf> in <module>
      6                  timestep=600,
      7                  timerange = [c, d],
----> 8                  ha_limit=np.pi)

~/rmextract/lib/python3.6/site-packages/RMextract/getRM.py in getRM(MS, server, prefix, ionexPath, earth_rot, timerange, use_azel, ha_limit, use_filter, **kwargs)
    185               flags[station].append(1)
    186             else:
--> 187                az,el = PosTools.getAzEl(pointing,time,position,ha_limit)
    188                if az==-1 and el==-1:
    189                   return

~/rmextract/lib/python3.6/site-packages/RMextract/PosTools.py in getAzEl(pointing, time, position, ha_limit)
    586             me.doframe(p)
    587             me.doframe(t)
--> 588             hadec=me.measure(phasedir,"HADEC")
    589             if abs(hadec['m0']['value'])>ha_limit:
    590                 print ("below horizon",tab.taql('calc ctod($time s)')[0],degrees(hadec['m0']['value']),degrees(hadec['m1']['value']))

~/rmextract/lib/python3.6/site-packages/casacore/measures/__init__.py in measure(self, v, rf, off)
    118                 if dq.is_quantity(v[key]):
    119                     v[key] = v[key].to_dict()
--> 120         return _measures.measure(self, v, rf, off)
    121 
    122     def direction(self, rf='', v0='0..', v1='90..', off=None):

RuntimeError: Cannot convert due to missing frame information

for reference, I'm trying to run:

a  = tab.taql('calc MJD("2021/04/25/19:33:12")')[0]*3600*24
b  = tab.taql('calc MJD("2021/04/25/22:29:16")')[0]*3600*24
rd = np.array([3.53925792, 0.53248541])
ms472='msdir/1619378472-cal.MS'
RMdict=grm.getRM(MS=ms472,
                 server="ftp://ftp.aiub.unibe.ch/CODE/",
                 prefix='CODG',
                 ionexPath="IONEXdata/",
                 radec=rd,
                 timestep=600,
                 timerange = [a, b],
                 ha_limit=np.pi)

I did find a bug when you have http in the server address, maybe this was what you ran into. Will be solved soon

Thanks, this was indeed exactly the bug I ran into with the server timeouts. 👍

[maaijke]

Ok, this seems to be an error of your measures data installation. One day I will remove all dependencies on casacore measures, until then either make sure all necessary measures data is in: /usr/share/casacore/data/ or change the path in your .casarc : measures.directory: ~/opt/casacore/data/ (for example). I am not an expert on casacore installation, however.

It might run smoothly if you remove the ha_limit option, but since I use casacore.measures elsewhere I cannot guarantee it does.

[maaijke]

I could not reproduce the issue with the ROBR files. I managed to download and use them (of the same date), so maybe somehow the download got interrupted. The best you can do in such case is to remove all recently downloaded IONEX data. Or use the overwrite option in getRM that I will add within a day.

[twillis449]

I just got IONEX files from AUIB (CODE) last week with no problems. I thought that to use CDDIS these days you had to have an account / password.

Also I don't have a whole lot of faith in URLLIB. I use PyCurl for data retrieval.

To avoid the CASA complaints about leap seconds you have to keep the CASA data directory up-to-data but as Maaijke saya this error is totally insignificant for ionosphere rotation measure analysis.

[gmloose]

I guess this issue can be closed. Remaining problems were with casacore, not with RMextract.