twillis449/albus_ionosphere

software to determine ionosphere TEC and RM from GPS receiver data

This software determines the ionosphere total electron content (TEC) over any location on the Earth as a function of location and time. It then uses the TEC and a model of the Earth's magnetc field to compute the ionosphere's effect on the Faraday Rotion Meaure (RM) observed for an astronomical radio source. The ionosphere's contribution to the RM can then be renoved. The software may be of interest to both radio astronomers and to ionosphere scientists. Test observations suggest that our analysis gives results consistent with those found by on-site site experiments that used local GPS receivers.

The software derives the TEC of the ionosphere by using publicly available observation data of Global Positioning System (GPS) satellites. It searches through a database of several thousand ground-based GPS receivers and then gets GPS receiver data from those stations located within a specified distance of the position of the telescope being used for the radio astronomy observation.

In principle, the Faraday rotation angle is a function of source direction and antenna position, but Faraday rotation is usually a large-scale effect and it may have approximately the same value across an entire telescope primary beam field of view (perhaps about one degree). For arrays smaller than a few kilometres, the rotation angle will usually also be the same for all stations. These assumptions reduce the number of independent parameters considerably, but they may break down as the observing wavelength gets longer due to the wavelength squared effect and increasing field of view, as well as when telescope arrays have longer baselines.

In reality, from the ground we usually cannot directly measure the distribution of the electrons along the line of sight nor directly measure the magnetic field strength as a function of position and direction. In order to calculate a rotation measure, many routines place all the electrons at some "standard" height and attach a magnetic field value from a model of the terrestrial field. In contrast, the software that we describe goes beyond this simple algorithm by distributing the electrons along the line of sight taking into account modern understanding of ionospheric physics, and employs a model of the terrestrial magnetic field that accounts for change of intensity and direction with height.

The software makes extensive use of python scripts and should work with both python2 and python3.

You need to install pycurl, astropy, pyephem or python casacore, numpy and matplotlib for the system to work. A number of support programs to handle RINEX files are also needed. These programs are specified in the INSTALL file.

More sites (especially Geosciences Australia) are now producing only RINEX3 files. For analysis, we use RINEX2 files. To convert RINEX3 to RINEX2 you need to get and install gfzrnx, available from https://gnss.gfz-potsdam.de/services/gfzrnx and RX3name (see http://acc.igs.org/software.html) Unfortunately these programs seem to be available only in binary format.

Unfortunately, due to security concerns, sites are also changing access methods from simple anonymous ftp to more secure procedures such as sftp etc and we are currently working on modifying our access procedures to reach such sites. secure access procedures

A somewhat more detained description of the software is given in the, as yet, unpublished paper twillis_ALBUS_paper.pdf available in this directory.

albus_ionosphere's Issues

Unknown ephemeris FTP site

Hello,

More of a question then a real issue,

I've installed ALBUS through docker (the easy way I just downloaded the repo and installed it using the pre-included docker file)

I get some errors when successfully running the South_Africa_experiment.py code (namely, BrokenPipeError: [Errno 32] Broken pipe and FileNotFoundError: [Errno 2] No such file or directory: 'RX3name': 'RX3name') but these seem to have no effect, and the code runs well and outputs the TEC file.

However, If I push the start and end times too close to the current date, I get the title error Unknown ephemeris FTP site, and the code crashes. I interpret this as the RINEX file doesn't exist on the trignet servers. Looking at the website http://www.trignet.co.za/, it says that the RINEX files are posted daily, but as of today (Feb 8 2024) the most recent file I can get is Jan 27 2024.

So my question is, is this grace period standard? And is there any way to get around it?

Thanks for the help and the Dockerfile; it made it really easy to install.
Andrew

Fitting using IRI model returns RM = 0

Switching to mode IRI returns zeros in the RM column. Here is an example:
PIM:

seq  rel_time time_width El         Az         STEC           RM (rad/m2)   VTEC factor   STEC_error
0 : 0 -300.0 300.0 26.968437859474882 15.701915887957162 8.372198927124062 -0.280817685993113 0.5245814881133594 2.880392735586948

IRI:

seq  rel_time time_width El         Az         STEC           RM (rad/m2)   VTEC factor   STEC_error
0 : 0 -300.0 300.0 26.968437859474882 15.701915887957162 20.35657702357253 0.0 0.5245814881133594 1.6350473013541083

Detailed comparison of TEC solutions from RMextract vs ALBUS vs AIPS TECOR

We need to systematically work through differences in these packages. I cannot get the CASA scripts to work so I think I'm just going to focus my time on ALBUS as I'm quite convinced it is the better solution compared to the vTEC IONEX based systems currently available.

There is now fully calibrated data for

• VLA L Moon 2021 (AIPS)
• MeerKAT UHF and L Moon 2021 (AIPS & CASA)
• VLA Venus (AIPS)

On top of having multiple possible RINEX sites (seemingly giving the same results for MeerKAT, multiple IONEX databases exist for the VLA, including jplg, which seems to consistently give reasonable results (better than emrg), although not for MeerKAT it seems. All this assumes a vTEC correction which I doubt will work anywhere far from Zenith (which probably explains why it does not work for 3C286 for MeerKAT)

How to curl the IONEX databases from Eric (AIPS release 31DEC22)

; TECOR
;---------------------------------------------------------------
;! Calculate ionospheric delay and Faraday rotation corrections
;# Task Calibration VLA VLBI
;-----------------------------------------------------------------------
;;  Copyright (C) 1998; 2001; 2003-2007, 2015, 2019-2020, 2022
;;  Associated Universities, Inc. Washington DC, USA.
;;
;;  This program is free software; you can redistribute it and/or
;;  modify it under the terms of the GNU General Public License as
;;  published by the Free Software Foundation; either version 2 of
;;  the License, or (at your option) any later version.
;;
;;  This program is distributed in the hope that it will be useful,
;;  but WITHOUT ANY WARRANTY; without even the implied warranty of
;;  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
;;  GNU General Public License for more details.
;;
;;  You should have received a copy of the GNU General Public
;;  License along with this program; if not, write to the Free
;;  Software Foundation, Inc., 675 Massachusetts Ave, Cambridge,
;;  MA 02139, USA.
;;
;;  Correspondence concerning AIPS should be addressed as follows:
;;         Internet email: [email protected].
;;         Postal address: AIPS Project Office
;;                         National Radio Astronomy Observatory
;;                         520 Edgemont Road
;;                         Charlottesville, VA 22903-2475 USA
;-----------------------------------------------------------------------
TECOR     LLLLLLLLLLLLUUUUUUUUUUUU CCCCCCCCCCCCCCCCCCCCCCCCCCCCC
TECOR     Task to calibrate ionospheric delay and Faraday rotn.
INNAME                             Input UV file name (name)
INCLASS                            Input UV file name (class)
INSEQ             0.0     9999.0   Input UV file name (seq. #)
INDISK            0.0        9.0   Input UV file disk unit #
INFILE                             IONEX file name
NFILES            0.0      999.0   # of INFILES
SUBARR                             Subarray to correct
                                     (0 => all)
ANTENNAS                           Antennas to correct
                                     (all 0 => all)
GAINVER                            Input CL table version
GAINUSE                            Output CL table version
APARM                              Switches
                                     (1) if > 0 correct for
                                         dispersive delay
                                     (2) Follow ionosphere
                                         factor
----------------------------------------------------------------
TECOR
Task:   Derives corrections for ionospheric Faraday rotation and
        dispersive delay from maps of total electron content in
        IONEX format.
        There are two procedures that will download the files
        automatically (including calculating the day number) and
        run TECOR.  They are VLBATECR and VLATECR and they are
        part of the procedure packages VLBAUTIL and VLAPROCS
        respectively.  See HELP VLBATECR or VLATECR for more info.

        Moving sources are supported if a PO table is present.
Adverbs
  INNAME.....Input UV file name (name).      Standard defaults.
  INCLASS....Input UV file name (class).     Standard defaults.
  INSEQ......Input UV file name (seq. #).    0 => highest.
  INDISK.....Disk drive # of input UV file.  0 => any.
  INFILE.....The name of the file containing the TEC data. The
             filename should be given in the usual AIPS style as
             a logical directory name followed by a colon and
             the name of the file or the directory path.
             The file should contain TEC maps spanning the time
             range covered by the input CL table and should be in
             IONEX format.  If NFILES > 1 then the name MUST be in
             standard format CCCCdddC.yyC where C can be any
             character, ddd is the day number, and yy is the year.
             Also the name given in INFILE must be the first day.
             See EXPLAIN TECOR for more details.
  NFILES.....Number of IONEX files to use.  Note that if NFILES > 1
             then the INFILE must be in a standard format and the
             file in INFILE must be the first day of those to be
             loaded.  0 => 1
  SUBARR.....Subarray to correct.            0 => all.
  ANTENNAS...A list of antennas to correct. If all entries are
             zero then all antennas in the subarray or subarrays
             selected by SUBARR will be corrected. If all of the
             entries are positive or zero then only those
             antennas with numbers corresponding to the positive
             entries will be corrected. If any entries are
             negative then all of the antennas in the selected
             subarray or subarrays will be corrected except
             those with numbers corresponding the the absolute
             values of the non-zero entries in this array.
  GAINVER....Version number of the input CL table to use.
                                             0 => highest
  GAINUSE....Version number of output CL table. This must not
             be the number of an existing CL table.
                                             0 => highest + 1
  APARM......Miscellaneous settings and switches.
             (1) Enable (> 0.5) or disable (<= 0.5) dispersive
                 delay corrections.
             (2) In principle, it is thought that the ionosphere
                 remains approximately fixed wrt the Sun.  Thus, the
                 task should predict which ionosphere is now in the
                 direction of the source by applying a time correction
                 to the apparent longitude while interpolating between
                 the table values, which are only given every two
                 hours.  This sometimes seems in fact to do odd
                 things.  This parameter allows you to do only a
                 fraction of the time correction, from epsilon to 1.0.
                 0 -> 1; < 0 -> 0.  A correction of zero is equivalent
                 to a model in which the ionosphere rotates with the
                 earth.
----------------------------------------------------------------
TECOR:   Task that corrects ionospheric Faraday rotation and
         dispersive delay using maps of ionospheric electron
         content.
         There are two procedures that will download the files
         automatically (including calulating the day number) and
         run TECOR.  They are VLBATECR and VLATECR and they are
         part of the procedure packages VLBAUTIL and VLAPROCS
         respectively.  See HELP VLBATECR or VLATECR for more info.
DOCUMENTOR: Chris Flatters, Amy Mioduszewski, NRAO
RELATED PROGRAMS: LDGPS, GPSDL, PCAL, LPCAL, VLBATECR, VLATECR
                  most programs that apply calibration

TECOR reads a set of maps of ionospheric electron content that
covers the time range of the observations and uses this data
to calculate the ionospheric Faraday rotation and, if you
request it, the dispersive delay introduced by the ionosphere.
The dispersive delay that is in put in the CL table (column
DISP) is the dispersive delay at a wavelength of 1 meter.  The
input data is expected to be text files in IONEX format.

The IONEX format is a standard format for the interchange of
ionosphere maps and is used by NASA's crustal dynamics data
interchange system (CDDIS), among others. Each IONEX file
contains a series of maps of the zenith total electron content
of the ionosphere as a function of geographical latitude and
longitude taken at different times. TECOR expects the input
file to contain maps covering the whole range of times contained
in the AIPS data file.

Some IONEX files may contain maps that do not cover all of the
antenna in the array. If this is the case then you should
use the ANTENNAS adverb to prevent TECOR from trying to
determine corrections for antennas outside the covered area
otherwise the calibration records for these antennas will be
marked as invalid. Note that you only need to determine the
ionospheric Faraday rotation for your phase reference
antenna to calibrate polarization; differences between the
antennas are removed by self-calibration.

Dispersive delay corrections depend greatly on the accuracy of
the input file and should therefore be checked carefully if they
are turned on.

Ionospheric models are available from the CDDIS data archive.
There are 5 IONEX files produced by different groups for each
day archived at CDDIS.  These groups are the Jet Propulsion
Laboratory (JPL), the Center for Orbit Determination in Europe
(CODE), the Geodetic Survey Division of Natural resources Canada
(EMR), the ESOC Ionosphere Monitoring Facility (ESA), and the
Technical University of Catalonia (UPC).  At this time, all the
files provide maps every 2 hours.  Since these things change so
rapidly, it is difficult to recommend one file over another.
For the best possible results try all the files and use the one
that seems to work the best (this can be judged using VPLOT).
If you only try one, there is a marginal recommendation that you
use the JPL or CODE files.  As mentioned above, there might come
a time where one or more of these files will become significantly
better than the others so downloading all the files might be a
good idea even if you do not plan of trying them all (they are
not large files).


WHICH FILES TO DOWNLOAD
=======================
The IONEX files being read in must bracket the entire
experiment.

BEFORE November 3, 2002 the IONEX files on the CDDIS archive
spanned the time from 1:00 to 23:00.  So if you have an experiment
which starts before 1:00 or ends after 23:00 then you need to get
IONEX files for the day of the experiment and the previous day or
next day respectively.

AFTER November 3, 2002 the IONEX files on the CDDIS archive
span the time from 0:00 to 24:00.  So you can now download
files for only the days on which your experiment has data.

*********************************************************************
*  NOTE: There are two procedures that will download the files      *
*        automatically (including calulating the day number) and    *
*        run TECOR.  They are VLBATECR and VLATECR and are          *
*        part of the procedure packages VLBAUTIL and VLAPROCS       *
*        respectively.  See HELP VLBATECR or VLATECR for more info. *
*********************************************************************

Procedures VLATECR and VLBATECR determine which file(s) are needed and
attempt to download them automatically.  If there is a problem with the
automatic download then you can download the IONEX file(s) manually and
set INFILE and NFILES:

Sometime in October 2020, for new requirements, you can consult

https://cddis.nasa.gov/Data_and_Derived_Products/CDDIS_Archive_Access.html

You may have to create an account with cddis in order to download in
some other fashion.

The command that does work from VLBATECR is (all on one line!)

curl -u anonymous:[email protected] --ftp-ssl ftp://gdc.cddis.eosdis.nasa.gov/
   gps/products/ionex/YYYY/DDD/jplgDDD0.YYi.Z > /tmp/jplgDDD0.YYi.Z

(no blanks between gov/ and gps) where YYYY is the year e.g. 2015, YY is
the last 2 digits of the year e.g. 15, and DDD is the day number within
the year.



TO DOWNLOAD THE FILES:
======================

1) ftp cddis.gsfc.nasa.gov; login as anonymous.
   -- note that there is a mirror site at igs.ensg.ign.fr

2) cd pub/gps/products/ionex/YYYY/DDD (YYYY-year; DDD-day number)
   -- for igs.ensg.ign.fr: cd pub/igs/iono/YYYY/DDD

3) prompt -1

4) mget jpl*.Z  (or get jplgDDD0.YYi.Z)  (not all *.Z !)

5) uncompress above file(s) before using TECOR

If more than one file is needed to cover the time period of
the experiment:

	-- download all the files

        -- make sure the files have the format CCCCdddC.yyC; where
           C is any character, ddd is the three digit day number and
           yy is the 2 digit year.  This is designed to work with
           the standard file names that the CDDIS data archive
           uses.  Note that the C's must be the same for each file
           and they must all be in the same directory.

        -- use the FIRST file as the INFILE

        -- set NFILES to the number of files

Example:

   INFILE 'FITS:CODG1230.99I'
   NFILES 4

   Will expect to find files:

   CODG1230.99I
   CODG1240.99I
   CODG1250.99I
   CODG1260.99I

   in directory $FITS.

Positron detections with Python 3.x

As reported by @o-smirnov there is an issue at present somewhere due to the porting to python3 where numerical issues crept in.

Current codebase produce

Old codebase produce

I will diff the codebases to try track this down

Potential future python software

I suspect much of the ALBUS python code could be replaced by Astropy.

There are python modules that could probably replace a lot of the James Anderson handing of RINEX files. GeoRinex (see https://pypi.org/project/georinex/) looks to be very powerful.

There are a couple of python packages that interface to the IGRF magnetic field model See https://github.com/klaundal/ppigrf and PyIGRF https://pypi.org/project/pyIGRF/

The International Reference Ionosphere has a python interface - see https://github.com/rilma/pyIRI2016

Add lunar ephemeris

Need a flag in iono functions to take in "special body"

Compare IGRFv13 vs v12 vs AIPS

Making a note to link this to an investigation of upgrading the current IGRFv12 inside casacore as mentioned on casacore/casacore#1114

Fails to build with Python 3.6

Hmm think we have a regression - forgot to raise at the time that the change made recently breaks Python 3.6 (Ubuntu 18.04 default - and the provided Dockerfile that bases off this) support - specifically this: 0074a6b

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/optsoft/ALBUS/share/python/MS_Iono_functions.py", line 56, in <module>
    import GPS_stations
  File "/optsoft/ALBUS/share/python/GPS_stations.py", line 520, in <module>
    GPS_stations = fill_standard_stations()
  File "/optsoft/ALBUS/share/python/GPS_stations.py", line 479, in fill_standard_stations
    GPS_dict = fill_GPS_station_dict(GPS_dict, sys_file)
  File "/optsoft/ALBUS/share/python/GPS_stations.py", line 246, in fill_GPS_station_dict
    raise e
  File "/optsoft/ALBUS/share/python/GPS_stations.py", line 221, in fill_GPS_station_dict
    for line in fp:
  File "/usr/lib/python3.6/encodings/ascii.py", line 26, in decode
    return codecs.ascii_decode(input, self.errors)[0]
UnicodeDecodeError: 'ascii' codec can't decode byte 0xc3 in position 5073: ordinal not in range(128)

Continuous integration tests

@twillis449 can you give me a few examples of ALBUS scripts for say MeerKAT and VLA at a given date/time range and coordinate (say 3C286). I will spend some time in the morning to set up continuous integration testing before I attempt any changes to the code base

Expose ionospheric model parameters for testing

Going to expose parameters to switch between PIM and IRI with the other parameters as envionment varaiables so we don't have to rebuild this each time

RI_G03 issue with single station observation

I've noticed that if we only have data from a single GPS station, the RI_G03 solution produces reasonable TEC parameters but the corresponding RMs are garbage. On the other hand, RI_G01 solutions produce reasonable values. This needs to be investigated.

Dockerize & continuous integration

First pass to get things going in a systematic way would be to Dockerize this into a working shippable installation

I see there is some test suites that we can execute although I'm not sure about coverage. I can set up CI testing on that before making any modification to the codebase