From 8da39eb8e5880427ab78352d75f711dc5c118c62 Mon Sep 17 00:00:00 2001 From: "Alec Thomson (S&A, Kensington WA)" Date: Mon, 11 Dec 2023 16:50:49 +0800 Subject: [PATCH 1/2] Add config --- .pre-commit-config.yaml | 28 ++++++++++++++++++++++++++++ 1 file changed, 28 insertions(+) create mode 100644 .pre-commit-config.yaml diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml new file mode 100644 index 0000000..0f743fd --- /dev/null +++ b/.pre-commit-config.yaml @@ -0,0 +1,28 @@ +repos: +- repo: https://github.com/pre-commit/pre-commit-hooks + rev: v4.5.0 + hooks: + - id: trailing-whitespace + - id: end-of-file-fixer + - id: check-yaml + - id: check-added-large-files +- repo: https://github.com/psf/black + rev: 23.11.0 + hooks: + - id: black +- repo: https://github.com/PyCQA/isort + rev: 5.12.0 + hooks: + - id: isort + args: ["--profile=black"] +ci: + autofix_commit_msg: | + [pre-commit.ci] auto fixes from pre-commit.com hooks + + for more information, see https://pre-commit.ci + autofix_prs: true + autoupdate_branch: '' + autoupdate_commit_msg: '[pre-commit.ci] pre-commit autoupdate' + autoupdate_schedule: weekly + skip: [] + submodules: false From 5a43c5f9aff7b163877a722b8f93c693fc229d3d Mon Sep 17 00:00:00 2001 From: "Alec Thomson (S&A, Kensington WA)" Date: Mon, 11 Dec 2023 16:51:19 +0800 Subject: [PATCH 2/2] Fix formatting --- FRion/correct.py | 456 +++++++++-------- FRion/download_IONEX_CDDIS.py | 114 +++-- FRion/predict.py | 908 ++++++++++++++++++++-------------- README.md | 1 - docs/source/conf.py | 23 +- docs/source/index.rst | 62 +-- setup.py | 56 +-- 7 files changed, 909 insertions(+), 711 deletions(-) diff --git a/FRion/correct.py b/FRion/correct.py index 13859b3..1e8e9d6 100644 --- a/FRion/correct.py +++ b/FRion/correct.py @@ -6,29 +6,32 @@ to data cubes (either in memory, or in FITS files.) The complex polarization (Q + iU) is divided by the predicted ionospheric -modulation to produce corrected values that should have the effect of the +modulation to produce corrected values that should have the effect of the ionosphere removed. These can then be saved to new Stokes Q and U FITS files. The main functions below do not do anything specific to handle -very large FITS files gracefully. It may not perform efficiently when file +very large FITS files gracefully. It may not perform efficiently when file sizes are comparable to the amount of available RAM. An alternative mode, -`apply_correction_large_cube()`, has been developed to reduce the memory +`apply_correction_large_cube()`, has been developed to reduce the memory footprint required. """ -import numpy as np -from astropy.io import fits as pf import os import sys from math import floor -def apply_correction_to_files(Qfile,Ufile,predictionfile,Qoutfile,Uoutfile, - overwrite=False): - """ This function combines all the individual steps needed to apply a +import numpy as np +from astropy.io import fits as pf + + +def apply_correction_to_files( + Qfile, Ufile, predictionfile, Qoutfile, Uoutfile, overwrite=False +): + """This function combines all the individual steps needed to apply a correction to a set of Q and U FITS cubes and save the results. The user should supply the paths to all the files as specified. - + Args: Qfile (str): filename of uncorrected Stokes Q FITS cube Ufile (str): filename of uncorrected Stokes U FITS cube @@ -36,45 +39,46 @@ def apply_correction_to_files(Qfile,Ufile,predictionfile,Qoutfile,Uoutfile, Qoutfile (str): filename for corrected Stokes Q FITS cube. Uoutfile (str): filename for corrected Stokes U FITS cube. overwrite (bool): overwrite Stokes Q/U files if they already exist? [False] - + """ - - #Get all data: - frequencies,theta=read_prediction(predictionfile) - Qdata,Udata,Qheader,Uheader=readData(Qfile,Ufile) - - #Checks for data consistency. - if (Qdata.shape != Udata.shape): + + # Get all data: + frequencies, theta = read_prediction(predictionfile) + Qdata, Udata, Qheader, Uheader = readData(Qfile, Ufile) + + # Checks for data consistency. + if Qdata.shape != Udata.shape: raise Exception("Q and U files don't have same dimensions.") if Qdata.shape[0] != theta.size: - raise Exception("Prediction file does not have same number of channels as FITS cube.") - #Currently this doesn't actually check that the frequencies are the same, - #just that the number of channels is the same. Should this be a more - #strict check? + raise Exception( + "Prediction file does not have same number of channels as FITS cube." + ) + # Currently this doesn't actually check that the frequencies are the same, + # just that the number of channels is the same. Should this be a more + # strict check? + # Apply correction + Qcorr, Ucorr = correct_cubes(Qdata, Udata, theta) - #Apply correction - Qcorr,Ucorr=correct_cubes(Qdata,Udata,theta) - - #Save results - write_corrected_cubes(Qoutfile,Uoutfile,Qcorr,Ucorr,Qheader,Uheader, - overwrite=overwrite) - + # Save results + write_corrected_cubes( + Qoutfile, Uoutfile, Qcorr, Ucorr, Qheader, Uheader, overwrite=overwrite + ) def read_prediction(filename): """Read in frequencies and ionospheric predictions from text file. - + Returns: tuple containing - - -frequencies (array): frequencies of each channel (Hz); + + -frequencies (array): frequencies of each channel (Hz); -theta (array): ionospheric modulation for each channel - + """ - (frequencies,real,imag)=np.genfromtxt(filename,unpack=True) - theta=real+1.j*imag + (frequencies, real, imag) = np.genfromtxt(filename, unpack=True) + theta = real + 1.0j * imag return frequencies, theta @@ -83,58 +87,59 @@ def find_freq_axis(header): Input: header: a Pyfits header object. Returns the axis number (as recorded in the FITS file, **NOT** in numpy ordering.) Returns 0 if the frequency axis cannot be found. - + """ - freq_axis=0 #Default for 'frequency axis not identified' - #Check for frequency axes. Because I don't know what different formatting - #I might get ('FREQ' vs 'OBSFREQ' vs 'Freq' vs 'Frequency'), convert to - #all caps and check for 'FREQ' anywhere in the axis name. - for i in range(1,header['NAXIS']+1): #Check each axis in turn. + freq_axis = 0 # Default for 'frequency axis not identified' + # Check for frequency axes. Because I don't know what different formatting + # I might get ('FREQ' vs 'OBSFREQ' vs 'Freq' vs 'Frequency'), convert to + # all caps and check for 'FREQ' anywhere in the axis name. + for i in range(1, header["NAXIS"] + 1): # Check each axis in turn. try: - if 'FREQ' in header['CTYPE'+str(i)].upper(): - freq_axis=i + if "FREQ" in header["CTYPE" + str(i)].upper(): + freq_axis = i except: - pass #The try statement is needed for if the FITS header does not - # have CTYPE keywords. + pass # The try statement is needed for if the FITS header does not + # have CTYPE keywords. return freq_axis -def readData(Qfilename,Ufilename): - """Open the Stokes Q and U input cubes (from the supplied - file names) and return data-access variables and the header. +def readData(Qfilename, Ufilename): + """Open the Stokes Q and U input cubes (from the supplied + file names) and return data-access variables and the header. Axes are re-ordered so that frequency is first, beyond that the number and ordering of axes doesn't matter. Uses the memmap functionality so that data isn't read into data; variables are just handles to access the data on disk. Returns the header from the Q file, the U file's header is ignored. - - """ - - hdulistQ=pf.open(Qfilename,memmap=True) - Qheader=hdulistQ[0].header - Qdata=hdulistQ[0].data - hdulistU=pf.open(Ufilename,memmap=True) - Udata=hdulistU[0].data - Uheader=hdulistU[0].header - - N_dim=Qheader['NAXIS'] #Get number of axes - - freq_axis=find_freq_axis(Qheader) - #If the frequency axis isn't the last one, rotate the array until it is. - #Recall that pyfits reverses the axis ordering, so we want frequency on - #axis 0 of the numpy array. + + """ + + hdulistQ = pf.open(Qfilename, memmap=True) + Qheader = hdulistQ[0].header + Qdata = hdulistQ[0].data + hdulistU = pf.open(Ufilename, memmap=True) + Udata = hdulistU[0].data + Uheader = hdulistU[0].header + + N_dim = Qheader["NAXIS"] # Get number of axes + + freq_axis = find_freq_axis(Qheader) + # If the frequency axis isn't the last one, rotate the array until it is. + # Recall that pyfits reverses the axis ordering, so we want frequency on + # axis 0 of the numpy array. if freq_axis != 0 and freq_axis != N_dim: - Qdata=np.moveaxis(Qdata,N_dim-freq_axis,0) - Udata=np.moveaxis(Udata,N_dim-freq_axis,0) + Qdata = np.moveaxis(Qdata, N_dim - freq_axis, 0) + Udata = np.moveaxis(Udata, N_dim - freq_axis, 0) - return Qdata, Udata, Qheader, Uheader -def write_corrected_cubes(Qoutputname,Uoutputname,Qcorr,Ucorr,Qheader,Uheader,overwrite=False): - """ Write the corrected Q and U data to FITS files. Copies the supplied +def write_corrected_cubes( + Qoutputname, Uoutputname, Qcorr, Ucorr, Qheader, Uheader, overwrite=False +): + """Write the corrected Q and U data to FITS files. Copies the supplied header, adding a note to the history saying that the correction was applied. - + Inputs: Qoutputname (str): filename to write corrected Stoke Q data to. Uoutputname (str): filename to write corrected Stoke U data to. @@ -142,54 +147,57 @@ def write_corrected_cubes(Qoutputname,Uoutputname,Qcorr,Ucorr,Qheader,Uheader,ov Ucorr (array): corrected Stokes U data [Q/U]header: Astropy FITS header objects that describes the data overwrite (bool): overwrite Stokes Q/U files if they already exist? [False] - - """ - Qoutput_header=Qheader.copy() - Qoutput_header.add_history('Corrected for ionospheric Faraday rotation using FRion.') - Uoutput_header=Uheader.copy() - Uoutput_header.add_history('Corrected for ionospheric Faraday rotation using FRion.') - - #Get data back to original axis order, if necessary. - N_dim=Qoutput_header['NAXIS'] #Get number of axes - freq_axis=find_freq_axis(Qoutput_header) + """ + Qoutput_header = Qheader.copy() + Qoutput_header.add_history( + "Corrected for ionospheric Faraday rotation using FRion." + ) + Uoutput_header = Uheader.copy() + Uoutput_header.add_history( + "Corrected for ionospheric Faraday rotation using FRion." + ) + + # Get data back to original axis order, if necessary. + N_dim = Qoutput_header["NAXIS"] # Get number of axes + freq_axis = find_freq_axis(Qoutput_header) if freq_axis != 0: - Qcorr=np.moveaxis(Qcorr,0,N_dim-freq_axis) - Ucorr=np.moveaxis(Ucorr,0,N_dim-freq_axis) - + Qcorr = np.moveaxis(Qcorr, 0, N_dim - freq_axis) + Ucorr = np.moveaxis(Ucorr, 0, N_dim - freq_axis) - pf.writeto(Qoutputname,Qcorr.astype('float32'),Qoutput_header,overwrite=overwrite) - pf.writeto(Uoutputname,Ucorr.astype('float32'),Uoutput_header,overwrite=overwrite) + pf.writeto( + Qoutputname, Qcorr.astype("float32"), Qoutput_header, overwrite=overwrite + ) + pf.writeto( + Uoutputname, Ucorr.astype("float32"), Uoutput_header, overwrite=overwrite + ) - - -def correct_cubes(Qdata,Udata,theta): +def correct_cubes(Qdata, Udata, theta): """Applies the ionospheric Faraday rotation correction to the Stokes Q/U data, derotating the polarization angle and renormalizing to remove depolarization. Note that this will amplify the noise present in the data, particularly if the depolarization is large (\|theta\| is small). - + Inputs: Qdata (array): uncorrected Stokes Q data, frequency axis first Udata (array): uncorrected Stokes U data, frequency axis first theta (1D array): ionospheric modulation, per frequency - + Returns: Qcorr (array): corrected Stokes Q data, same axis ordering Ucorr (array): corrected Stokes U data, same axis ordering """ - - Pdata=Qdata+1.j*Udata #Input complex polarization - arrshape=np.array(Pdata.shape) #the correction needs the same number of - arrshape[:]=1 #axes as the input data - arrshape[0]=theta.size #(but they can all be degenerate) - Pcorr=np.true_divide(Pdata,np.reshape(theta,arrshape)) - Qcorr=Pcorr.real - Ucorr=Pcorr.imag - - return Qcorr,Ucorr + Pdata = Qdata + 1.0j * Udata # Input complex polarization + arrshape = np.array(Pdata.shape) # the correction needs the same number of + arrshape[:] = 1 # axes as the input data + arrshape[0] = theta.size # (but they can all be degenerate) + Pcorr = np.true_divide(Pdata, np.reshape(theta, arrshape)) + Qcorr = Pcorr.real + Ucorr = Pcorr.imag + + return Qcorr, Ucorr def progress(width, percent): @@ -200,27 +208,27 @@ def progress(width, percent): marks = floor(width * (percent / 100.0)) spaces = floor(width - marks) - loader = ' [' + ('=' * int(marks)) + (' ' * int(spaces)) + ']' + loader = " [" + ("=" * int(marks)) + (" " * int(spaces)) + "]" sys.stdout.write("%s %d%%\r" % (loader, percent)) if percent >= 100: sys.stdout.write("\n") sys.stdout.flush() - -def apply_correction_large_cube(Qfile,Ufile,predictionfile,Qoutfile,Uoutfile, - overwrite=False): +def apply_correction_large_cube( + Qfile, Ufile, predictionfile, Qoutfile, Uoutfile, overwrite=False +): """Functions as apply_correction_to_files, but for files too large to hold in memory. Combines the correct_cubes() and write_corrected_cubes() steps into a single function so that it can operate on smaller pieces of data at one time. - - This function combines all the individual steps needed to apply a + + This function combines all the individual steps needed to apply a correction to a set of Q and U FITS cubes and save the results. The user should supply the paths to all the files as specified. This function is less flexible in terms of axis ordering: it assumes two spatial axes, and a frequency axis in position 3 or 4. - + Args: Qfile (str): filename of uncorrected Stokes Q FITS cube Ufile (str): filename of uncorrected Stokes U FITS cube @@ -228,84 +236,97 @@ def apply_correction_large_cube(Qfile,Ufile,predictionfile,Qoutfile,Uoutfile, Qoutfile (str): filename for corrected Stokes Q FITS cube. Uoutfile (str): filename for corrected Stokes U FITS cube. overwrite (bool): overwrite Stokes Q/U files if they already exist? [False] - + """ + # Get all data: + frequencies, theta = read_prediction(predictionfile) + + hdulistQ = pf.open(Qfile, memmap=True) + Qheader = hdulistQ[0].header + Qdata = hdulistQ[0].data + hdulistU = pf.open(Ufile, memmap=True) + Udata = hdulistU[0].data + Uheader = hdulistU[0].header - #Get all data: - frequencies,theta=read_prediction(predictionfile) - - hdulistQ=pf.open(Qfile,memmap=True) - Qheader=hdulistQ[0].header - Qdata=hdulistQ[0].data - hdulistU=pf.open(Ufile,memmap=True) - Udata=hdulistU[0].data - Uheader=hdulistU[0].header - - N_dim=Qheader['NAXIS'] #Get number of axes - freq_axis=find_freq_axis(Qheader) - #Checks for data consistency. - if (Qdata.shape != Udata.shape): + N_dim = Qheader["NAXIS"] # Get number of axes + freq_axis = find_freq_axis(Qheader) + # Checks for data consistency. + if Qdata.shape != Udata.shape: raise Exception("Q and U files don't have same dimensions.") - if Qdata.shape[N_dim-freq_axis] != theta.size: - raise Exception("Prediction file does not have same number of channels as FITS cube.") - #Currently this doesn't actually check that the frequencies are the same, - #just that the number of channels is the same. Should this be a more - #strict check? + if Qdata.shape[N_dim - freq_axis] != theta.size: + raise Exception( + "Prediction file does not have same number of channels as FITS cube." + ) + # Currently this doesn't actually check that the frequencies are the same, + # just that the number of channels is the same. Should this be a more + # strict check? if Qdata.ndim != 3 and Qdata.ndim != 4: - raise Exception("Cube does not have 3 or 4 axes; only these are supported for large files.") - - - #Add correction to header history - Qoutput_header=Qheader.copy() - Qoutput_header.add_history('Corrected for ionospheric Faraday rotation using FRion.') - Uoutput_header=Uheader.copy() - Uoutput_header.add_history('Corrected for ionospheric Faraday rotation using FRion.') - - - #Deal with any existing output files: + raise Exception( + "Cube does not have 3 or 4 axes; only these are supported for large files." + ) + + # Add correction to header history + Qoutput_header = Qheader.copy() + Qoutput_header.add_history( + "Corrected for ionospheric Faraday rotation using FRion." + ) + Uoutput_header = Uheader.copy() + Uoutput_header.add_history( + "Corrected for ionospheric Faraday rotation using FRion." + ) + + # Deal with any existing output files: if (os.path.isfile(Qoutfile) or os.path.isfile(Uoutfile)) and not overwrite: raise Exception("Output file(s) aready exist.") if os.path.isfile(Qoutfile) and overwrite: os.remove(Qoutfile) if os.path.isfile(Uoutfile) and overwrite: os.remove(Uoutfile) - - #Create large blank files. This seems to produce file size complaints - #sometimes, but those seem harmless so far. - shape = tuple(Qoutput_header['NAXIS{0}'.format(ii)] for ii in range(1, Qoutput_header['NAXIS']+1)) - - Qoutput_header.tofile(Qoutfile) - with open(Qoutfile, 'rb+') as fobj: - fobj.seek(len(Qoutput_header.tostring()) + (np.product(shape) * np.abs(Qoutput_header['BITPIX']//8)) - 1) - fobj.write(b'\0') - - Uoutput_header.tofile(Uoutfile) - with open(Uoutfile, 'rb+') as fobj: - fobj.seek(len(Uoutput_header.tostring()) + (np.product(shape) * np.abs(Uoutput_header['BITPIX']//8)) - 1) - fobj.write(b'\0') - - - Qout_hdu=pf.open(Qoutfile,mode='update',memmap=True) - Uout_hdu=pf.open(Uoutfile,mode='update',memmap=True) + # Create large blank files. This seems to produce file size complaints + # sometimes, but those seem harmless so far. + shape = tuple( + Qoutput_header["NAXIS{0}".format(ii)] + for ii in range(1, Qoutput_header["NAXIS"] + 1) + ) + Qoutput_header.tofile(Qoutfile) + with open(Qoutfile, "rb+") as fobj: + fobj.seek( + len(Qoutput_header.tostring()) + + (np.product(shape) * np.abs(Qoutput_header["BITPIX"] // 8)) + - 1 + ) + fobj.write(b"\0") - for i in range(theta.size): #Iterate over channels. - if (N_dim==4) & (freq_axis == 3): - Pdata=Qdata[:,i]+1.j*Udata[:,i] #Input complex polarization + Uoutput_header.tofile(Uoutfile) + with open(Uoutfile, "rb+") as fobj: + fobj.seek( + len(Uoutput_header.tostring()) + + (np.product(shape) * np.abs(Uoutput_header["BITPIX"] // 8)) + - 1 + ) + fobj.write(b"\0") + + Qout_hdu = pf.open(Qoutfile, mode="update", memmap=True) + Uout_hdu = pf.open(Uoutfile, mode="update", memmap=True) + + for i in range(theta.size): # Iterate over channels. + if (N_dim == 4) & (freq_axis == 3): + Pdata = Qdata[:, i] + 1.0j * Udata[:, i] # Input complex polarization else: - Pdata=Qdata[i]+1.j*Udata[i] #Input complex polarization + Pdata = Qdata[i] + 1.0j * Udata[i] # Input complex polarization - Pcorr=np.true_divide(Pdata,theta[i]) - if (N_dim==4) & (freq_axis == 3): - Qout_hdu[0].data[:,i]=Pcorr.real - Uout_hdu[0].data[:,i]=Pcorr.imag + Pcorr = np.true_divide(Pdata, theta[i]) + if (N_dim == 4) & (freq_axis == 3): + Qout_hdu[0].data[:, i] = Pcorr.real + Uout_hdu[0].data[:, i] = Pcorr.imag else: - Qout_hdu[0].data[i]=Pcorr.real - Uout_hdu[0].data[i]=Pcorr.imag + Qout_hdu[0].data[i] = Pcorr.real + Uout_hdu[0].data[i] = Pcorr.imag - progress(40, i/theta.size*100) + progress(40, i / theta.size * 100) Qout_hdu.flush() Uout_hdu.flush() @@ -313,68 +334,93 @@ def apply_correction_large_cube(Qfile,Ufile,predictionfile,Qoutfile,Uoutfile, Uout_hdu.close() - - - - def command_line(): """When invoked from the command line, parse the input options to get the filenames and other parameters, then invoke apply_correction_to_files to run all the steps and save the output cubes. - + """ - + import argparse import os + descStr = """ Apply correction for ionospheric Faraday rotation to Stokes Q and U FITS cubes. Requires the file names for the input cubes, output cubes, and the - prediction file (which contains the ionospheric modulation per channel + prediction file (which contains the ionospheric modulation per channel in the cubes). """ - parser = argparse.ArgumentParser(description=descStr, - formatter_class=argparse.RawTextHelpFormatter) - parser.add_argument("fitsQ",metavar="fitsQ", - help="FITS cube containing (uncorrected) Stokes Q data.") - parser.add_argument("fitsU",metavar="fitsU", - help="FITS cube containing (uncorrected) Stokes U data.") - parser.add_argument("predictionfile",metavar="predictionfile", - help="File containing ionospheric prediction to be applied.") - parser.add_argument("outQ",metavar="Qcorrected", - help="Output filename for corrected Stokes Q cube.") - parser.add_argument("outU",metavar="Ucorrected", - help="Output filename for corrected Stokes U cube.") - parser.add_argument("-o",dest="overwrite",action="store_true", - help="Overwrite exising output files? [False]") - parser.add_argument("-L",dest="large",action="store_true", - help="Use large-file mode? (Reduced memory footprint) [False]") + parser = argparse.ArgumentParser( + description=descStr, formatter_class=argparse.RawTextHelpFormatter + ) + parser.add_argument( + "fitsQ", + metavar="fitsQ", + help="FITS cube containing (uncorrected) Stokes Q data.", + ) + parser.add_argument( + "fitsU", + metavar="fitsU", + help="FITS cube containing (uncorrected) Stokes U data.", + ) + parser.add_argument( + "predictionfile", + metavar="predictionfile", + help="File containing ionospheric prediction to be applied.", + ) + parser.add_argument( + "outQ", + metavar="Qcorrected", + help="Output filename for corrected Stokes Q cube.", + ) + parser.add_argument( + "outU", + metavar="Ucorrected", + help="Output filename for corrected Stokes U cube.", + ) + parser.add_argument( + "-o", + dest="overwrite", + action="store_true", + help="Overwrite exising output files? [False]", + ) + parser.add_argument( + "-L", + dest="large", + action="store_true", + help="Use large-file mode? (Reduced memory footprint) [False]", + ) args = parser.parse_args() - #Check for file existence. + # Check for file existence. if not os.path.isfile(args.fitsQ): raise Exception("Stokes Q file not found.") if not os.path.isfile(args.fitsU): raise Exception("Stokes U file not found.") - - #Pass file names into do-everything function (either basic or large-file, as - #set by user. + + # Pass file names into do-everything function (either basic or large-file, as + # set by user. if args.large: - apply_correction_large_cube(args.fitsQ,args.fitsU,args.predictionfile, - args.outQ,args.outU,overwrite=args.overwrite) + apply_correction_large_cube( + args.fitsQ, + args.fitsU, + args.predictionfile, + args.outQ, + args.outU, + overwrite=args.overwrite, + ) else: - apply_correction_to_files(args.fitsQ,args.fitsU,args.predictionfile, - args.outQ,args.outU,overwrite=args.overwrite) - - - - - - + apply_correction_to_files( + args.fitsQ, + args.fitsU, + args.predictionfile, + args.outQ, + args.outU, + overwrite=args.overwrite, + ) if __name__ == "__main__": command_line() - - diff --git a/FRion/download_IONEX_CDDIS.py b/FRion/download_IONEX_CDDIS.py index 1a882e0..158ca80 100644 --- a/FRion/download_IONEX_CDDIS.py +++ b/FRion/download_IONEX_CDDIS.py @@ -2,10 +2,10 @@ # -*- coding: utf-8 -*- """ This file contains the code necessary to download IONEX files from CDDIS. -It is intended to replace the getIONEX functionality of RMextract, since +It is intended to replace the getIONEX functionality of RMextract, since that code does not currently support CDDIS (which requires authentication). -The majority of the code is taken directly from RMextract. This has been +The majority of the code is taken directly from RMextract. This has been modified by Art Davydov and Anna Ordog to be compatible with CDDIS, and further modified by me to streamline the process. @@ -22,21 +22,21 @@ import datetime import logging import os + # AO, AD added these: -#import requests +# import requests logging.basicConfig(level=logging.ERROR) -def get_CDDIS_IONEXfile(time="2023/03/23/02:20:10.01", - prefix="jplg", - outpath='./', - overwrite=False): +def get_CDDIS_IONEXfile( + time="2023/03/23/02:20:10.01", prefix="jplg", outpath="./", overwrite=False +): """Get IONEX file with prefix from server for a given day Downloads files with given prefix from CDDIS, unzips and stores the data. - + CDDIS requires authentication to download files. Thus must be done as a .netrc file with the credentials in them. (https://cddis.nasa.gov/Data_and_Derived_Products/CreateNetrcFile.html) @@ -47,20 +47,18 @@ def get_CDDIS_IONEXfile(time="2023/03/23/02:20:10.01", outpath (string) : path where the data is stored overwrite (bool) : Do (not) overwrite existing data """ - - server="https://cddis.nasa.gov" - - + + server = "https://cddis.nasa.gov" + # AO, AD changed this to lower from upper to comply with cddis: - prefix=prefix.lower() + prefix = prefix.lower() if outpath[-1] != "/": outpath += "/" if not os.path.isdir(outpath): try: os.makedirs(outpath) except: - print("cannot create output dir for IONEXdata: %s", - outpath) + print("cannot create output dir for IONEXdata: %s", outpath) try: yy = int(time[2:4]) @@ -74,63 +72,73 @@ def get_CDDIS_IONEXfile(time="2023/03/23/02:20:10.01", day = time[2] mydate = datetime.date(year, month, day) dayofyear = mydate.timetuple().tm_yday - if not overwrite and os.path.isfile("%s%s%03d0.%02dI"%(outpath,prefix.upper(),dayofyear,yy)): - logging.info("FILE exists: %s%s%03d0.%02dI",outpath,prefix,dayofyear,yy) - return "%s%s%03d0.%02dI"%(outpath,prefix,dayofyear,yy) - - #If proxy url is given, enable proxy using pysocks + if not overwrite and os.path.isfile( + "%s%s%03d0.%02dI" % (outpath, prefix.upper(), dayofyear, yy) + ): + logging.info("FILE exists: %s%s%03d0.%02dI", outpath, prefix, dayofyear, yy) + return "%s%s%03d0.%02dI" % (outpath, prefix, dayofyear, yy) + + # If proxy url is given, enable proxy using pysocks try: from urllib import request except ImportError: import urllib2 as request - - #Naming conventions changed, at different times for different data sets. - #Depending on the data source/prefix, the filename has to be completely - #different for different date ranges. This needs to be manually coded, - #which means supporting a limited range of data sources. I'll try to support - #the main ones stored at CDDIS. - - if prefix == 'jplg' and mydate > datetime.date(2023,8,7): - filename=f"{prefix[0:3].upper()}0OPSFIN_{year}{dayofyear}0000_01D_02H_GIM.INX.gz" - elif prefix == 'codg' and mydate > datetime.date(2022,11,26): - filename=f"{prefix[0:3].upper()}0OPSFIN_{year}{dayofyear}0000_01D_01H_GIM.INX.gz" - elif prefix == 'igsg' and mydate > datetime.date(2022,11,26): - filename=f"{prefix[0:3].upper()}0OPSFIN_{year}{dayofyear}0000_01D_02H_GIM.INX.gz" - elif prefix == 'casg' and mydate > datetime.date(2022,12,31): - filename=f"{prefix[0:3].upper()}0OPSFIN_{year}{dayofyear}0000_01D_30M_GIM.INX.gz" - elif prefix == 'esag': - if mydate <= datetime.date(2023,2,4): - raise Exception("ESA did not publish GIM.INX.gz files prior to 2023 Feb 04.") + # Naming conventions changed, at different times for different data sets. + # Depending on the data source/prefix, the filename has to be completely + # different for different date ranges. This needs to be manually coded, + # which means supporting a limited range of data sources. I'll try to support + # the main ones stored at CDDIS. + + if prefix == "jplg" and mydate > datetime.date(2023, 8, 7): + filename = ( + f"{prefix[0:3].upper()}0OPSFIN_{year}{dayofyear}0000_01D_02H_GIM.INX.gz" + ) + elif prefix == "codg" and mydate > datetime.date(2022, 11, 26): + filename = ( + f"{prefix[0:3].upper()}0OPSFIN_{year}{dayofyear}0000_01D_01H_GIM.INX.gz" + ) + elif prefix == "igsg" and mydate > datetime.date(2022, 11, 26): + filename = ( + f"{prefix[0:3].upper()}0OPSFIN_{year}{dayofyear}0000_01D_02H_GIM.INX.gz" + ) + elif prefix == "casg" and mydate > datetime.date(2022, 12, 31): + filename = ( + f"{prefix[0:3].upper()}0OPSFIN_{year}{dayofyear}0000_01D_30M_GIM.INX.gz" + ) + elif prefix == "esag": + if mydate <= datetime.date(2023, 2, 4): + raise Exception( + "ESA did not publish GIM.INX.gz files prior to 2023 Feb 04." + ) else: - filename=f"{prefix[0:3].upper()}0OPSFIN_{year}{dayofyear}0000_01D_02H_GIM.INX.gz" + filename = ( + f"{prefix[0:3].upper()}0OPSFIN_{year}{dayofyear}0000_01D_02H_GIM.INX.gz" + ) else: - filename=f"{prefix}{dayofyear:03d}0.{yy:02d}i.Z" + filename = f"{prefix}{dayofyear:03d}0.{yy:02d}i.Z" - - url = "https://cddis.nasa.gov/archive/gnss/products/ionex/%4d/%03d/%s"%(year,dayofyear,filename) + url = "https://cddis.nasa.gov/archive/gnss/products/ionex/%4d/%03d/%s" % ( + year, + dayofyear, + filename, + ) # Download IONEX file, make sure output format is old-format name and uppercase. - fname = outpath+'/'+(f"{prefix}{dayofyear:03d}0.{yy:02d}i.Z").upper() - print("Downloading ",url) - + fname = outpath + "/" + (f"{prefix}{dayofyear:03d}0.{yy:02d}i.Z").upper() + print("Downloading ", url) os.system(f'wget --auth-no-challenge -O {fname} "{url}"') - ###### gunzip files - if fname[-2:].upper()==".Z": + if fname[-2:].upper() == ".Z": command = "gunzip -dc %s > %s" % (fname, fname[:-2]) retcode = os.system(command) if retcode: raise RuntimeError("Could not run '%s'" % command) else: os.remove(fname) - fname=fname[:-2] - #returns filename of uncompressed file + fname = fname[:-2] + # returns filename of uncompressed file return fname - - - - diff --git a/FRion/predict.py b/FRion/predict.py index ece03e2..5a9ac73 100644 --- a/FRion/predict.py +++ b/FRion/predict.py @@ -12,7 +12,7 @@ Other ionosphere RM codes are available (ionFR, ALBUS) are available, but RMextract was selected for its comparative ease of install and use. -RMextract relies on external maps of Total Electron Content (TEC). As of +RMextract relies on external maps of Total Electron Content (TEC). As of version 1.1, the default is to get the TEC data from CDDIS (https://cddis.nasa.gov/Data_and_Derived_Products/GNSS/atmospheric_products.html). This requires an account and a .netrc file with the credentials @@ -36,20 +36,21 @@ import RMextract.getRM as RME from RMextract import getIONEX as ionex except: - #This is the easiest solution to the documentation problem: - #This code needs to be importable when RMextract isn't installed, for - #ReadTheDocs to work. Getting RMextract to install properly in RTD is too - #much work, so this is my workaround. - print('Cannot import RMextract. Continuing import, but will fail if called.') -from astropy.time import Time,TimeDelta -import numpy as np -from astropy.coordinates import EarthLocation,SkyCoord,Angle, UnknownSiteException + # This is the easiest solution to the documentation problem: + # This code needs to be importable when RMextract isn't installed, for + # ReadTheDocs to work. Getting RMextract to install properly in RTD is too + # much work, so this is my workaround. + print("Cannot import RMextract. Continuing import, but will fail if called.") +import astropy import astropy.units as u +import numpy as np +from astropy.coordinates import Angle, EarthLocation, SkyCoord, UnknownSiteException +from astropy.time import Time, TimeDelta + from FRion.correct import find_freq_axis -import astropy from FRion.download_IONEX_CDDIS import get_CDDIS_IONEXfile -C = 2.99792458e8 # Speed of light [m/s] +C = 2.99792458e8 # Speed of light [m/s] def predict(): @@ -63,6 +64,7 @@ def predict(): Call with the -h flag to see command line options. """ import argparse + descStr = """ Calculate ionospheric Faraday rotation and predict time-integrated effect as a function of frequency. @@ -70,61 +72,111 @@ def predict(): parameters from a supplied FITS cube or PSRFITS file, if it has the correct keywords, otherwise from those parameters must be supplied on the command line. - + By default, gets ionosphere TEC data from CDDIS (https://cddis.nasa.gov/Data_and_Derived_Products/GNSS/atmospheric_products.html) This requires an account and a .netrc file with the credentials to run. (https://cddis.nasa.gov/Data_and_Derived_Products/CreateNetrcFile.html) """ - parser = argparse.ArgumentParser(description=descStr, - formatter_class=argparse.RawTextHelpFormatter) - parser.add_argument("-F",dest="fits",default=None,metavar='FILENAME', - help="FITS cube relevant information in header.") - parser.add_argument("-d", dest="times", nargs=2,type=str,default=None, - metavar=('START','END'), - help="start and end time strings.") - parser.add_argument("-c", dest=("freq_parms"),nargs=3,default=None,type=float, - metavar=('MINFREQ','MAXFREQ','CHANNELWIDTH'), - help=("Generate channel frequencies (in Hz): \n" - " minfreq, maxfreq, channel_width")) - parser.add_argument("-t",dest='telescope_name',type=str,default=None, - help="Telescope name") - parser.add_argument("-T",dest='telescope_coords',nargs=3,type=float,default=None, - metavar=('LONG','LAT','ALT'), - help="Telescope coordinates:\n lat[deg],long[deg], altitude[m].") - parser.add_argument('-p',dest='pointing',nargs=2,type=float,default=None, - metavar=('RA','DEC'), - help="Pointing center: RA[deg], Dec[deg]") - parser.add_argument("-s", dest='savefile',type=str,default=None,metavar='POLFILE', - help="Filename to save ionosphere data to.") - parser.add_argument("-S", dest='savefig',type=str,default=None,metavar='FIGFILE', - help="Filename to save the plots to. Entering 'screen' plots to the screen.") - parser.add_argument("--timestep",dest='timestep',default=600.,type=float, - help="Timestep for ionospheric prediction, in seconds. Default = 600") + parser = argparse.ArgumentParser( + description=descStr, formatter_class=argparse.RawTextHelpFormatter + ) + parser.add_argument( + "-F", + dest="fits", + default=None, + metavar="FILENAME", + help="FITS cube relevant information in header.", + ) + parser.add_argument( + "-d", + dest="times", + nargs=2, + type=str, + default=None, + metavar=("START", "END"), + help="start and end time strings.", + ) + parser.add_argument( + "-c", + dest=("freq_parms"), + nargs=3, + default=None, + type=float, + metavar=("MINFREQ", "MAXFREQ", "CHANNELWIDTH"), + help=( + "Generate channel frequencies (in Hz): \n" + " minfreq, maxfreq, channel_width" + ), + ) + parser.add_argument( + "-t", dest="telescope_name", type=str, default=None, help="Telescope name" + ) + parser.add_argument( + "-T", + dest="telescope_coords", + nargs=3, + type=float, + default=None, + metavar=("LONG", "LAT", "ALT"), + help="Telescope coordinates:\n lat[deg],long[deg], altitude[m].", + ) + parser.add_argument( + "-p", + dest="pointing", + nargs=2, + type=float, + default=None, + metavar=("RA", "DEC"), + help="Pointing center: RA[deg], Dec[deg]", + ) + parser.add_argument( + "-s", + dest="savefile", + type=str, + default=None, + metavar="POLFILE", + help="Filename to save ionosphere data to.", + ) + parser.add_argument( + "-S", + dest="savefig", + type=str, + default=None, + metavar="FIGFILE", + help="Filename to save the plots to. Entering 'screen' plots to the screen.", + ) + parser.add_argument( + "--timestep", + dest="timestep", + default=600.0, + type=float, + help="Timestep for ionospheric prediction, in seconds. Default = 600", + ) args = parser.parse_args() + start_time = None + end_time = None + freq_arr = None + telescope = None + ra = None + dec = None - start_time=None - end_time=None - freq_arr=None - telescope=None - ra=None - dec=None - - - #If a FITS file is present, try to fill in any missing keywords. - #But since FITS headers can be very different, it may be that not all - #keywords can be found. + # If a FITS file is present, try to fill in any missing keywords. + # But since FITS headers can be very different, it may be that not all + # keywords can be found. if args.fits is not None: - start_time,end_time,freq_arr,telescope,ra,dec=get_parms_from_FITS(args.fits) + start_time, end_time, freq_arr, telescope, ra, dec = get_parms_from_FITS( + args.fits + ) - #Any parametrs taken from FITS header can be overridden by manual inputs: + # Any parametrs taken from FITS header can be overridden by manual inputs: if args.times is not None: - start_time=Time(args.times[0]) - end_time=Time(args.times[1]) + start_time = Time(args.times[0]) + end_time = Time(args.times[1]) if args.freq_parms is not None: - freq_arr=np.arange(args.freq_parms[0],args.freq_parms[1],args.freq_parms[2]) + freq_arr = np.arange(args.freq_parms[0], args.freq_parms[1], args.freq_parms[2]) if args.telescope_name is not None: telescope = get_telescope_coordinates(args.telescope_name) @@ -132,33 +184,57 @@ def predict(): telescope = get_telescope_coordinates(args.telescope_coords) if args.pointing is not None: - ra=args.pointing[0] - dec=args.pointing[1] - - #Check that all parameters are set: - missing_parms=[] - if (start_time is None): missing_parms.append('Start time') - if (end_time is None): missing_parms.append('End time') - if (freq_arr is None): missing_parms.append('Frequency array') - if (telescope is None): missing_parms.append('Telescope') - if (ra is None): missing_parms.append('Pointing center') + ra = args.pointing[0] + dec = args.pointing[1] + + # Check that all parameters are set: + missing_parms = [] + if start_time is None: + missing_parms.append("Start time") + if end_time is None: + missing_parms.append("End time") + if freq_arr is None: + missing_parms.append("Frequency array") + if telescope is None: + missing_parms.append("Telescope") + if ra is None: + missing_parms.append("Pointing center") if len(missing_parms) > 0: - print("\n\nMissing parameters:",missing_parms,'\n') + print("\n\nMissing parameters:", missing_parms, "\n") raise Exception("Cannot continue without parameters listed above.") - times,RMs,theta=calculate_modulation(start_time, end_time, freq_arr, telescope, - ra,dec, timestep=args.timestep,ionexPath='./IONEXdata/') + times, RMs, theta = calculate_modulation( + start_time, + end_time, + freq_arr, + telescope, + ra, + dec, + timestep=args.timestep, + ionexPath="./IONEXdata/", + ) if args.savefile is not None: - write_modulation(freq_arr,theta,args.savefile) + write_modulation(freq_arr, theta, args.savefile) if args.savefig is not None: - generate_plots(times,RMs,theta,freq_arr,position=[ra,dec],savename=args.savefig) + generate_plots( + times, RMs, theta, freq_arr, position=[ra, dec], savename=args.savefig + ) -def calculate_modulation(start_time, end_time, freq_array, telescope_location, - ra,dec, timestep=600.,ionexPath='./IONEXdata/', **kwargs): +def calculate_modulation( + start_time, + end_time, + freq_array, + telescope_location, + ra, + dec, + timestep=600.0, + ionexPath="./IONEXdata/", + **kwargs +): """Calculate the ionospheric FR modulation (time-averaged effect), as a function of frequency, for a given observation (time, location, target direction). @@ -184,7 +260,7 @@ def calculate_modulation(start_time, end_time, freq_array, telescope_location, ionexPath (str, default='./IONEXdata/'): path to download IONEX files to for ionosphere calculations. **kwargs: additional keyword arguments to pass to RMextract.getRM() - + Returns: tuple containing @@ -197,38 +273,44 @@ def calculate_modulation(start_time, end_time, freq_array, telescope_location, """ - #Convert frequencies to have units if needed. + # Convert frequencies to have units if needed. if type(freq_array) == astropy.units.quantity.Quantity: - frequencies=freq_array + frequencies = freq_array else: - frequencies=freq_array*u.Hz - - #Calculation of the time-dependent RMs. - times,RMs=get_RM(start_time, end_time, telescope_location, - ra,dec, timestep=timestep,ionexPath=ionexPath, **kwargs) + frequencies = freq_array * u.Hz + # Calculation of the time-dependent RMs. + times, RMs = get_RM( + start_time, + end_time, + telescope_location, + ra, + dec, + timestep=timestep, + ionexPath=ionexPath, + **kwargs, + ) - #Compute the time-integrated change in polarization. - theta=numeric_integration(times.mjd*86400.,RMs,frequencies.to(u.Hz).value) + # Compute the time-integrated change in polarization. + theta = numeric_integration(times.mjd * 86400.0, RMs, frequencies.to(u.Hz).value) - #Verify that we are not in a regime where numerical instabilities might occur. - check_numeric_problems(RMs, frequencies.to(u.Hz).value,theta) + # Verify that we are not in a regime where numerical instabilities might occur. + check_numeric_problems(RMs, frequencies.to(u.Hz).value, theta) - - return times,RMs, theta + return times, RMs, theta def get_RM( - start_time, - end_time, + start_time, + end_time, telescope_location, ra, - dec, - timestep=600., - ionexPath='./IONEXdata/', + dec, + timestep=600.0, + ionexPath="./IONEXdata/", pre_download=True, - prefix='jplg', - server='http://cddis.gsfc.nasa.gov', + prefix="jplg", + server="http://cddis.gsfc.nasa.gov", **kwargs ): """ @@ -258,7 +340,7 @@ def get_RM( server (str, default='http://cddis.gsfc.nasa.gov'): server to download IONEX files from. **kwargs: additional keyword arguments to pass to RMextract.getRM() - + Returns: tuple containing @@ -267,108 +349,113 @@ def get_RM( -RMs (array): vector of RM values computed for each time step """ - #If necessary, convert telescope name into telescope location object: + # If necessary, convert telescope name into telescope location object: if type(telescope_location) != EarthLocation: - telescope_location=get_telescope_coordinates(telescope_location) + telescope_location = get_telescope_coordinates(telescope_location) - #RMextract wants time ranges in MJD seconds: - timerange=[Time(start_time).mjd*86400.0, - Time(end_time).mjd*86400.0] + # RMextract wants time ranges in MJD seconds: + timerange = [Time(start_time).mjd * 86400.0, Time(end_time).mjd * 86400.0] - #Extract telescope coordinates into expected format (geodetic x,y,z): - telescope_coordinates=[telescope_location.x.value, - telescope_location.y.value, - telescope_location.z.value] + # Extract telescope coordinates into expected format (geodetic x,y,z): + telescope_coordinates = [ + telescope_location.x.value, + telescope_location.y.value, + telescope_location.z.value, + ] - #Handle all forms of angle input: + # Handle all forms of angle input: if type(ra) == float or type(ra) == int: - ra_angle=Angle(ra,'deg') + ra_angle = Angle(ra, "deg") elif type(ra) == astropy.units.quantity.Quantity: - ra_angle=Angle(ra) + ra_angle = Angle(ra) elif type(ra) == astropy.coordinates.angles.Angle: - ra_angle=ra + ra_angle = ra else: - raise Exception("""RA input object type not recognized. - Only astropy.coordinates Angle, astropy.units Quantity, or float (in deg) allowed.""") + raise Exception( + """RA input object type not recognized. + Only astropy.coordinates Angle, astropy.units Quantity, or float (in deg) allowed.""" + ) if type(dec) == float or type(dec) == int: - dec_angle=Angle(dec,'deg') + dec_angle = Angle(dec, "deg") elif type(dec) == astropy.units.quantity.Quantity: - dec_angle=Angle(dec) + dec_angle = Angle(dec) elif type(dec) == astropy.coordinates.angles.Angle: - dec_angle=dec + dec_angle = dec else: - raise Exception("""Dec input object type not recognized. - Only astropy.coordinates Angle, astropy.units Quantity, or float (in deg) allowed.""") + raise Exception( + """Dec input object type not recognized. + Only astropy.coordinates Angle, astropy.units Quantity, or float (in deg) allowed.""" + ) - #Handle all forms of time input: + # Handle all forms of time input: if type(timestep) == float or type(timestep) == int: - timestep_Delta=TimeDelta(timestep*u.second) + timestep_Delta = TimeDelta(timestep * u.second) elif type(timestep) == astropy.units.quantity.Quantity: - timestep_Delta=TimeDelta(timestep) + timestep_Delta = TimeDelta(timestep) elif type(timestep) == astropy.time.core.TimeDelta: - timestep_Delta=timestep + timestep_Delta = timestep else: - raise Exception("""Timestep input object type not recognized. - Only astropy.time TimeDelta, astropy.units Quantity, or float (in seconds) allowed.""") - + raise Exception( + """Timestep input object type not recognized. + Only astropy.time TimeDelta, astropy.units Quantity, or float (in seconds) allowed.""" + ) - #Pre-download the IONEX data from CDDIS, to work around the RMextract lack - #of support for CDDIS downloads. + # Pre-download the IONEX data from CDDIS, to work around the RMextract lack + # of support for CDDIS downloads. if pre_download: - _predownload_CDDIS(start_time, end_time,prefix,outpath=ionexPath) - - #Get RMExtract to generate its RM predictions - predictions=RME.getRM(ionexPath=ionexPath, - radec=[ra_angle.rad,dec_angle.rad], - timestep=timestep_Delta.sec, - timerange = timerange, - stat_positions=[telescope_coordinates,], - prefix=prefix, - server=server, - **kwargs, - ) - - #predictions dictionary contains STEC, Bpar, BField, AirMass, elev, azimuth + _predownload_CDDIS(start_time, end_time, prefix, outpath=ionexPath) + + # Get RMExtract to generate its RM predictions + predictions = RME.getRM( + ionexPath=ionexPath, + radec=[ra_angle.rad, dec_angle.rad], + timestep=timestep_Delta.sec, + timerange=timerange, + stat_positions=[ + telescope_coordinates, + ], + prefix=prefix, + server=server, + **kwargs, + ) + + # predictions dictionary contains STEC, Bpar, BField, AirMass, elev, azimuth # RM, times, timestep, station_names, stat_pos, flags, reference_time - times=Time(predictions['times']/86400.,format='mjd') - RMs=np.squeeze(predictions['RM']['st1']) - + times = Time(predictions["times"] / 86400.0, format="mjd") + RMs = np.squeeze(predictions["RM"]["st1"]) return times, RMs -def _predownload_CDDIS(start_time,end_time,prefix='jplg',outpath='./IONEXdata/'): +def _predownload_CDDIS(start_time, end_time, prefix="jplg", outpath="./IONEXdata/"): """Downloads the IONEX maps for the required dates from CDDIS. This must occur before RMextract is invoked, otherwise RMextract will try to use its own download tool which is broken. - + Will try to download all the days between the start and end dates. Is slightly greedy (downloading more days than) may be necessary) as a safety against edge cases. - + """ from math import ceil, floor - - - #Work out how many days need to be downloaded: - start_date=Time(start_time) - end_date=Time(end_time) - - #Download each day one by one: - for day_mjd in range(floor(start_date.mjd)-1,ceil(end_date.mjd)+1): - day = Time(day_mjd,format='mjd') - fname=get_CDDIS_IONEXfile(time=day.to_value('isot'), - prefix=prefix, - outpath=outpath, - overwrite=False) - if fname == -1: - raise Exception('Something has gone wrong in downloading IONEX data.') + # Work out how many days need to be downloaded: + start_date = Time(start_time) + end_date = Time(end_time) + # Download each day one by one: + for day_mjd in range(floor(start_date.mjd) - 1, ceil(end_date.mjd) + 1): + day = Time(day_mjd, format="mjd") + fname = get_CDDIS_IONEXfile( + time=day.to_value("isot"), prefix=prefix, outpath=outpath, overwrite=False + ) + if fname == -1: + raise Exception("Something has gone wrong in downloading IONEX data.") -def numeric_integration(times,RMs,freq_array): + +def numeric_integration(times, RMs, freq_array): """Numerical integration of the time-varying ionospheric polariation modulation. Testing has shown that numerical integration is accurate to better than 1% accuracy except where depolarization is extreme (>99%). @@ -382,19 +469,19 @@ def numeric_integration(times,RMs,freq_array): """ from scipy.integrate import simps - l2_arr=(C/freq_array)**2 - z=np.exp(2.j*np.outer(l2_arr,RMs)) - #Scipy's numerical integrators can't handle complex numbers, so the + l2_arr = (C / freq_array) ** 2 + z = np.exp(2.0j * np.outer(l2_arr, RMs)) + + # Scipy's numerical integrators can't handle complex numbers, so the # integral needs to be broken into real and complex components. - real=simps(z.real,times,axis=1) - imag=simps(z.imag,times,axis=1) - theta=(real+1.j*imag)/(times[-1]-times[0]) + real = simps(z.real, times, axis=1) + imag = simps(z.imag, times, axis=1) + theta = (real + 1.0j * imag) / (times[-1] - times[0]) return theta - -def check_numeric_problems(RMs, freq_array,theta): +def check_numeric_problems(RMs, freq_array, theta): """Checks for conditions that might cause numeric instability in the time-integration, and warns the user if there might be concerns. @@ -409,41 +496,49 @@ def check_numeric_problems(RMs, freq_array,theta): """ import warnings - #Check for large jumps in RM/polarization angle between steps. - #These can cause the numeric integrator to not catch angle wraps. - longest_l2=(C/np.min(freq_array))**2 - max_deltaRM=np.max(np.diff(RMs)) - max_delta_polangle=longest_l2*max_deltaRM #in radians + # Check for large jumps in RM/polarization angle between steps. + # These can cause the numeric integrator to not catch angle wraps. + longest_l2 = (C / np.min(freq_array)) ** 2 + max_deltaRM = np.max(np.diff(RMs)) + max_delta_polangle = longest_l2 * max_deltaRM # in radians if max_delta_polangle > 0.5: - warnings.warn(("\nLarge variations in RM between points, which may " - "introduce numerical errors.\n" - "Consider trying a smaller timestep.")) - - #Warn about very low values of theta (very strong depolarization) + warnings.warn( + ( + "\nLarge variations in RM between points, which may " + "introduce numerical errors.\n" + "Consider trying a smaller timestep." + ) + ) + + # Warn about very low values of theta (very strong depolarization) # as these can probably not be corrected reliably. if np.min(np.abs(theta)) < 0.02: - warnings.warn(("\nExtreme depolarization predicted (>98%). " - "Corrected polarization will almost certainly not " - "be trustworthy in affected channels.")) + warnings.warn( + ( + "\nExtreme depolarization predicted (>98%). " + "Corrected polarization will almost certainly not " + "be trustworthy in affected channels." + ) + ) elif np.min(np.abs(theta)) < 0.1: - warnings.warn(("\nSignificant depolarization predicted (>90%). " - "Errors in corrected polarization are likely to be " - "very large in some channels.")) - + warnings.warn( + ( + "\nSignificant depolarization predicted (>90%). " + "Errors in corrected polarization are likely to be " + "very large in some channels." + ) + ) -def write_modulation(freq_array,theta,filename): +def write_modulation(freq_array, theta, filename): """Saves predicted ionospheric modulation to a text file. File has two columns, whitespace-delimited. Args: freq_array (array): channel frequencies (in Hz) theta (array): ionospheric (complex) modulation at each frequency - filename (str): file path to save data to. -""" - np.savetxt(filename, list(zip(freq_array,theta.real,theta.imag))) - - + filename (str): file path to save data to.""" + np.savetxt(filename, list(zip(freq_array, theta.real, theta.imag))) def get_telescope_coordinates(telescope): @@ -460,26 +555,27 @@ def get_telescope_coordinates(telescope): manually coded in as a temporary measure. """ - if type(telescope) == EarthLocation: #Pass EarthLocations through without processing + if ( + type(telescope) == EarthLocation + ): # Pass EarthLocations through without processing return telescope - elif type(telescope) == str: #Hardcoded coordinates for some telescopes. - if telescope == 'ASKAP': - lat = -1*26+42/60+15/3600 #degree - long = +1*116+39/60+32/3600 # degree - height = 381.0 # + elif type(telescope) == str: # Hardcoded coordinates for some telescopes. + if telescope == "ASKAP": + lat = -1 * 26 + 42 / 60 + 15 / 3600 # degree + long = +1 * 116 + 39 / 60 + 32 / 3600 # degree + height = 381.0 # else: return EarthLocation.of_site(telescope) - return EarthLocation(lat=lat*u.deg, lon=long*u.deg, height=height*u.m) + return EarthLocation(lat=lat * u.deg, lon=long * u.deg, height=height * u.m) elif (type(telescope) == tuple) or (type(telescope) == list): - return EarthLocation(lat=telescope[0]*u.deg, lon=telescope[1]*u.deg, - height=telescope[2]*u.m) - - - + return EarthLocation( + lat=telescope[0] * u.deg, + lon=telescope[1] * u.deg, + height=telescope[2] * u.m, + ) - -def generate_plots(times,RMs,theta,freq_array,position=None,savename=None): +def generate_plots(times, RMs, theta, freq_array, position=None, savename=None): """Makes a figure with two plots: the RM variation over time, and the (modulus of the) modulation as a function of frequency. If savename contains a string it will save the plots to that filename, @@ -495,33 +591,30 @@ def generate_plots(times,RMs,theta,freq_array,position=None,savename=None): savename (str): File path to save plot to; if 'screen' will send to display. """ - from matplotlib import pyplot as plt from matplotlib import dates as mdates - plot_times=times.plot_date + from matplotlib import pyplot as plt + plot_times = times.plot_date - fig,(ax1,ax2)=plt.subplots(2,1,figsize=(8,8)) - ax1.plot_date(plot_times,RMs,fmt='k.') - locator=mdates.AutoDateLocator(minticks=3,maxticks=7) + fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(8, 8)) + ax1.plot_date(plot_times, RMs, fmt="k.") + locator = mdates.AutoDateLocator(minticks=3, maxticks=7) formatter = mdates.ConciseDateFormatter(locator) ax1.xaxis.set_major_locator(locator) ax1.xaxis.set_major_formatter(formatter) - ax1.set_ylabel(r'$\phi_\mathrm{ion}$ [rad m$^{-2}$]') + ax1.set_ylabel(r"$\phi_\mathrm{ion}$ [rad m$^{-2}$]") - ax2.plot(freq_array,np.abs(theta),'k.') - ax2.set_xlabel('Frequency [Hz]') - ax2.set_ylabel(r'|$\Theta(\lambda^2)$|') + ax2.plot(freq_array, np.abs(theta), "k.") + ax2.set_xlabel("Frequency [Hz]") + ax2.set_ylabel(r"|$\Theta(\lambda^2)$|") if position is not None: - ax1.set_title("RA: {:.2f}°, Dec: {:.2f}°".format(position[0],position[1])) + ax1.set_title("RA: {:.2f}°, Dec: {:.2f}°".format(position[0], position[1])) if savename is not None: if savename == "screen": plt.show() else: - plt.savefig(savename,bbox_inches='tight') - - - + plt.savefig(savename, bbox_inches="tight") def get_parms_from_FITS(filename): @@ -542,92 +635,105 @@ def get_parms_from_FITS(filename): dec """ import astropy.io.fits as pf - hdulist=pf.open(filename) - header=hdulist[0].header - - - start_time=None - end_time=None - freq_arr=None - telescope=None - ra=None - dec=None - - #PSRFITS files have a different format, so split prpcesing depending on - #PSRFITS vs FITS image: - if 'FITSTYPE' in header.keys() and header['FITSTYPE']=='PSRFITS': - #PRSFITS code adapted from example by David Kaplan. - #I'm trying to avoid using a PSRFITS reader module, to minimize dependencies, - #at risk of not handling all variations of PSRFITS properly. - - if 'STT_IMJD' in header.keys(): - start_time=Time(float(header['STT_IMJD']) + - (float(header['STT_SMJD'])+ - float(header['STT_OFFS']))/float(86400), - format='mjd') + + hdulist = pf.open(filename) + header = hdulist[0].header + + start_time = None + end_time = None + freq_arr = None + telescope = None + ra = None + dec = None + + # PSRFITS files have a different format, so split prpcesing depending on + # PSRFITS vs FITS image: + if "FITSTYPE" in header.keys() and header["FITSTYPE"] == "PSRFITS": + # PRSFITS code adapted from example by David Kaplan. + # I'm trying to avoid using a PSRFITS reader module, to minimize dependencies, + # at risk of not handling all variations of PSRFITS properly. + + if "STT_IMJD" in header.keys(): + start_time = Time( + float(header["STT_IMJD"]) + + (float(header["STT_SMJD"]) + float(header["STT_OFFS"])) + / float(86400), + format="mjd", + ) try: - end_time = start_time + np.sum(hdulist['SUBINT'].data['TSUBINT']) * u.s + end_time = start_time + np.sum(hdulist["SUBINT"].data["TSUBINT"]) * u.s except: pass - if 'TELESCOP' in header.keys(): + if "TELESCOP" in header.keys(): try: - telescope = EarthLocation.of_site(header['TELESCOP']) + telescope = EarthLocation.of_site(header["TELESCOP"]) except UnknownSiteException: - telescope = EarthLocation.from_geocentric(header['ANT_X'], - header['ANT_Y'], - header['ANT_Z'], - unit=u.m) - - if 'RA' in header.keys(): - ra=Angle(hdulist[0].header['RA'],unit='hour') - if 'DEC' in header.keys(): - dec=Angle(hdulist[0].header['DEC'],unit='deg') - - if 'OBSFREQ' in header.keys() and 'OBSNCHAN' in header.keys() and 'OBSBW' in header.keys(): - chan_bw=header['OBSBW']/header['OBSNCHAN'] - - startchan= header['OBSFREQ'] - (header['OBSNCHAN']/2 - 1)*chan_bw - freq_arr=np.arange(header['OBSNCHAN'])*chan_bw+startchan - - - #FITS images/cubes: + telescope = EarthLocation.from_geocentric( + header["ANT_X"], header["ANT_Y"], header["ANT_Z"], unit=u.m + ) + + if "RA" in header.keys(): + ra = Angle(hdulist[0].header["RA"], unit="hour") + if "DEC" in header.keys(): + dec = Angle(hdulist[0].header["DEC"], unit="deg") + + if ( + "OBSFREQ" in header.keys() + and "OBSNCHAN" in header.keys() + and "OBSBW" in header.keys() + ): + chan_bw = header["OBSBW"] / header["OBSNCHAN"] + + startchan = header["OBSFREQ"] - (header["OBSNCHAN"] / 2 - 1) * chan_bw + freq_arr = np.arange(header["OBSNCHAN"]) * chan_bw + startchan + + # FITS images/cubes: else: - if 'DATE-OBS' in header.keys(): - start_time=Time(header['DATE-OBS']) - if 'DATE-OBS' in header.keys() and 'DURATION' in header.keys(): - end_time=Time(header['DATE-OBS'])+TimeDelta(header['DURATION'],format="sec") - - #For coordinates, the code will always use the middle pixel to derive - #the RA and Dec. Assumes position coordinates are in first 2 axes. - if 'RA' in header['CTYPE1']: - ra=header['CRVAL1']+header['CDELT1']*(header['NAXIS1']/2-header['CRPIX1']) - if 'DEC' in header['CTYPE2']: - dec=header['CRVAL2']+header['CDELT2']*(header['NAXIS2']/2-header['CRPIX2']) - if 'GLON' in header['CTYPE1'] and 'GLAT' in header['CTYPE2']: - #Support galactic coordinates, just in case: - gl=header['CRVAL1']+header['CDELT1']*(header['NAXIS1']/2-header['CRPIX1']) - gb=header['CRVAL2']+header['CDELT2']*(header['NAXIS2']/2-header['CRPIX2']) - position=SkyCoord(gl,gb,frame='galactic',unit='deg') - ra=position.fk5.ra.deg - dec=position.fk5.dec.deg - - if 'TELESCOP' in header.keys(): - telescope=header['TELESCOP'] - - freq_axis=str(find_freq_axis(header)) - if freq_axis != '0': - chan0=header['CRVAL'+freq_axis]-header['CDELT'+freq_axis]*(header['CRPIX'+freq_axis]-1) - chan_final=chan0+(header['NAXIS'+freq_axis]-1)*header['CDELT'+freq_axis] - freq_arr=np.linspace(chan0,chan_final,header['NAXIS'+freq_axis]) - - return start_time,end_time,freq_arr,telescope,ra,dec - - - - - + if "DATE-OBS" in header.keys(): + start_time = Time(header["DATE-OBS"]) + if "DATE-OBS" in header.keys() and "DURATION" in header.keys(): + end_time = Time(header["DATE-OBS"]) + TimeDelta( + header["DURATION"], format="sec" + ) + + # For coordinates, the code will always use the middle pixel to derive + # the RA and Dec. Assumes position coordinates are in first 2 axes. + if "RA" in header["CTYPE1"]: + ra = header["CRVAL1"] + header["CDELT1"] * ( + header["NAXIS1"] / 2 - header["CRPIX1"] + ) + if "DEC" in header["CTYPE2"]: + dec = header["CRVAL2"] + header["CDELT2"] * ( + header["NAXIS2"] / 2 - header["CRPIX2"] + ) + if "GLON" in header["CTYPE1"] and "GLAT" in header["CTYPE2"]: + # Support galactic coordinates, just in case: + gl = header["CRVAL1"] + header["CDELT1"] * ( + header["NAXIS1"] / 2 - header["CRPIX1"] + ) + gb = header["CRVAL2"] + header["CDELT2"] * ( + header["NAXIS2"] / 2 - header["CRPIX2"] + ) + position = SkyCoord(gl, gb, frame="galactic", unit="deg") + ra = position.fk5.ra.deg + dec = position.fk5.dec.deg + + if "TELESCOP" in header.keys(): + telescope = header["TELESCOP"] + + freq_axis = str(find_freq_axis(header)) + if freq_axis != "0": + chan0 = header["CRVAL" + freq_axis] - header["CDELT" + freq_axis] * ( + header["CRPIX" + freq_axis] - 1 + ) + chan_final = ( + chan0 + (header["NAXIS" + freq_axis] - 1) * header["CDELT" + freq_axis] + ) + freq_arr = np.linspace(chan0, chan_final, header["NAXIS" + freq_axis]) + + return start_time, end_time, freq_arr, telescope, ra, dec def timeseries(): @@ -642,6 +748,7 @@ def timeseries(): Call with the -h flag to see command line options. """ import argparse + descStr = """ Calculate ionospheric Faraday rotation as a function of time. Can determine the observation time, direction, and location parameters @@ -649,52 +756,98 @@ def timeseries(): otherwise from those parameters must be supplied on the command line. """ - parser = argparse.ArgumentParser(description=descStr, - formatter_class=argparse.RawTextHelpFormatter) - parser.add_argument("-F",dest="fits",default=None,metavar='FILENAME', - help="FITS cube relevant information in header.") - parser.add_argument("-d", dest="times", nargs=2,type=str,default=None, - metavar=('START','END'), - help="start and end time strings.") - parser.add_argument("-t",dest='telescope_name',type=str,default=None, - help="Telescope name") - parser.add_argument("-T",dest='telescope_coords',nargs=3,type=float,default=None, - metavar=('LONG','LAT','ALT'), - help="Telescope coordinates:\n lat[deg],long[deg], altitude[m].") - parser.add_argument('-p',dest='pointing',nargs=2,type=float,default=None, - metavar=('RA','DEC'), - help="Pointing center: RA[deg], Dec[deg]") - parser.add_argument("-s", dest='savefile',type=str,default=None,metavar='POLFILE', - help="Filename to save ionosphere data to.") - parser.add_argument("-f", dest='timeformat',type=str,default='mjd', - metavar='TIMEFORMAT', - help="Format for times, must be one from list at https://docs.astropy.org/en/stable/time/index.html#time-format \n Default is mjd, for human-readable try fits.") - parser.add_argument("-S", dest='savefig',type=str,default=None,metavar='FIGFILE', - help="Filename to save the plots to. Entering 'screen' plots to the screen.") - parser.add_argument("--timestep",dest='timestep',default=600.,type=float, - help="Timestep for ionospheric prediction, in seconds. Default = 600") + parser = argparse.ArgumentParser( + description=descStr, formatter_class=argparse.RawTextHelpFormatter + ) + parser.add_argument( + "-F", + dest="fits", + default=None, + metavar="FILENAME", + help="FITS cube relevant information in header.", + ) + parser.add_argument( + "-d", + dest="times", + nargs=2, + type=str, + default=None, + metavar=("START", "END"), + help="start and end time strings.", + ) + parser.add_argument( + "-t", dest="telescope_name", type=str, default=None, help="Telescope name" + ) + parser.add_argument( + "-T", + dest="telescope_coords", + nargs=3, + type=float, + default=None, + metavar=("LONG", "LAT", "ALT"), + help="Telescope coordinates:\n lat[deg],long[deg], altitude[m].", + ) + parser.add_argument( + "-p", + dest="pointing", + nargs=2, + type=float, + default=None, + metavar=("RA", "DEC"), + help="Pointing center: RA[deg], Dec[deg]", + ) + parser.add_argument( + "-s", + dest="savefile", + type=str, + default=None, + metavar="POLFILE", + help="Filename to save ionosphere data to.", + ) + parser.add_argument( + "-f", + dest="timeformat", + type=str, + default="mjd", + metavar="TIMEFORMAT", + help="Format for times, must be one from list at https://docs.astropy.org/en/stable/time/index.html#time-format \n Default is mjd, for human-readable try fits.", + ) + parser.add_argument( + "-S", + dest="savefig", + type=str, + default=None, + metavar="FIGFILE", + help="Filename to save the plots to. Entering 'screen' plots to the screen.", + ) + parser.add_argument( + "--timestep", + dest="timestep", + default=600.0, + type=float, + help="Timestep for ionospheric prediction, in seconds. Default = 600", + ) args = parser.parse_args() + start_time = None + end_time = None + freq_arr = None + telescope = None + ra = None + dec = None - start_time=None - end_time=None - freq_arr=None - telescope=None - ra=None - dec=None - - - #If a FITS file is present, try to fill in any missing keywords. - #But since FITS headers can be very different, it may be that not all - #keywords can be found. + # If a FITS file is present, try to fill in any missing keywords. + # But since FITS headers can be very different, it may be that not all + # keywords can be found. if args.fits is not None: - start_time,end_time,freq_arr,telescope,ra,dec=get_parms_from_FITS(args.fits) + start_time, end_time, freq_arr, telescope, ra, dec = get_parms_from_FITS( + args.fits + ) - #Any parametrs taken from FITS header can be overridden by manual inputs: + # Any parametrs taken from FITS header can be overridden by manual inputs: if args.times is not None: - start_time=Time(args.times[0]) - end_time=Time(args.times[1]) - + start_time = Time(args.times[0]) + end_time = Time(args.times[1]) if args.telescope_name is not None: telescope = get_telescope_coordinates(args.telescope_name) @@ -702,40 +855,46 @@ def timeseries(): telescope = get_telescope_coordinates(args.telescope_coords) if args.pointing is not None: - ra=args.pointing[0] - dec=args.pointing[1] - - #Check that all parameters are set: - missing_parms=[] - if (start_time is None): missing_parms.append('Start time') - if (end_time is None): missing_parms.append('End time') - if (telescope is None): missing_parms.append('Telescope') - if (ra is None): missing_parms.append('Pointing center') + ra = args.pointing[0] + dec = args.pointing[1] + + # Check that all parameters are set: + missing_parms = [] + if start_time is None: + missing_parms.append("Start time") + if end_time is None: + missing_parms.append("End time") + if telescope is None: + missing_parms.append("Telescope") + if ra is None: + missing_parms.append("Pointing center") if len(missing_parms) > 0: - print("\n\nMissing parameters:",missing_parms,'\n') + print("\n\nMissing parameters:", missing_parms, "\n") raise Exception("Cannot continue without parameters listed above.") - times,RMs=get_RM(start_time, end_time, telescope, - ra,dec, timestep=args.timestep,ionexPath='./IONEXdata/') + times, RMs = get_RM( + start_time, + end_time, + telescope, + ra, + dec, + timestep=args.timestep, + ionexPath="./IONEXdata/", + ) if args.savefile is not None: - write_timeseries(times,RMs,args.savefile,timeformat=args.timeformat) + write_timeseries(times, RMs, args.savefile, timeformat=args.timeformat) else: - print('Times: RM:') - for tm, rm in zip(times.to_value(args.timeformat),RMs): - print(tm,rm) - + print("Times: RM:") + for tm, rm in zip(times.to_value(args.timeformat), RMs): + print(tm, rm) if args.savefig is not None: - generate_plots(times,RMs,position=[ra,dec],savename=args.savefig) - - - - + generate_plots(times, RMs, position=[ra, dec], savename=args.savefig) -def write_timeseries(times,RMs,filename,timeformat='mjd'): +def write_timeseries(times, RMs, filename, timeformat="mjd"): """Saves the predicted ionospheric RMs as a function of time to a text file. File has two columns, whitespace-delimited (Be aware that some time formats will add whitespace to the time column.) @@ -751,14 +910,11 @@ def write_timeseries(times,RMs,filename,timeformat='mjd'): """ fout = open(filename, "w") - for tm, rm in zip(times.to_value(timeformat),RMs): - print(tm,rm, file=fout) + for tm, rm in zip(times.to_value(timeformat), RMs): + print(tm, rm, file=fout) - - - -def plot_timeseries(times,RMs,position=None,savename=None): +def plot_timeseries(times, RMs, position=None, savename=None): """Makes a figure plotting the RM variation over time, If savename contains a string it will save the plots to that filename (unless the name is 'screen', in which case it will plot on screen), @@ -772,38 +928,28 @@ def plot_timeseries(times,RMs,position=None,savename=None): savename (str): File path to save plot to; if 'screen' will send to display. """ - from matplotlib import pyplot as plt from matplotlib import dates as mdates - plot_times=times.plot_date + from matplotlib import pyplot as plt + plot_times = times.plot_date - fig,ax1=plt.subplots(1,1,figsize=(8,8)) - ax1.plot_date(plot_times,RMs,fmt='k.') - locator=mdates.AutoDateLocator(minticks=3,maxticks=7) + fig, ax1 = plt.subplots(1, 1, figsize=(8, 8)) + ax1.plot_date(plot_times, RMs, fmt="k.") + locator = mdates.AutoDateLocator(minticks=3, maxticks=7) formatter = mdates.ConciseDateFormatter(locator) ax1.xaxis.set_major_locator(locator) ax1.xaxis.set_major_formatter(formatter) - ax1.set_ylabel(r'$\phi_\mathrm{ion}$ [rad m$^{-2}$]') + ax1.set_ylabel(r"$\phi_\mathrm{ion}$ [rad m$^{-2}$]") if position is not None: - ax1.set_title("RA: {:.2f}°, Dec: {:.2f}°".format(position[0],position[1])) + ax1.set_title("RA: {:.2f}°, Dec: {:.2f}°".format(position[0], position[1])) if savename is not None: if savename == "screen": plt.show() else: - plt.savefig(savename,bbox_inches='tight') - + plt.savefig(savename, bbox_inches="tight") if __name__ == "__main__": predict() - - - - - - - - - diff --git a/README.md b/README.md index 70e9ad5..1c43268 100644 --- a/README.md +++ b/README.md @@ -13,4 +13,3 @@ This packages original author, and maintainer as of Mar 2022, is Cameron Van Eck Please submit bug reports and feature requests fo the GitHub issues page, and feel free to email me with questions or comments. More information on the Canadian Initiative for Radio Astronomy Data Analysis (CIRADA) can be found at cirada.ca. - diff --git a/docs/source/conf.py b/docs/source/conf.py index 2ddef53..05964ab 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -12,17 +12,18 @@ # import os import sys -sys.path.insert(0, os.path.abspath('../')) + +sys.path.insert(0, os.path.abspath("../")) # -- Project information ----------------------------------------------------- -project = 'FRion' -copyright = '2021, Cameron Van Eck' -author = 'Cameron Van Eck' +project = "FRion" +copyright = "2021, Cameron Van Eck" +author = "Cameron Van Eck" # The full version, including alpha/beta/rc tags -release = '1.0' +release = "1.0" # -- General configuration --------------------------------------------------- @@ -30,11 +31,10 @@ # Add any Sphinx extension module names here, as strings. They can be # extensions coming with Sphinx (named 'sphinx.ext.*') or your custom # ones. -extensions = ['sphinx.ext.autodoc','sphinx.ext.napoleon' -] +extensions = ["sphinx.ext.autodoc", "sphinx.ext.napoleon"] # Add any paths that contain templates here, relative to this directory. -templates_path = ['_templates'] +templates_path = ["_templates"] # List of patterns, relative to source directory, that match files and # directories to ignore when looking for source files. @@ -47,12 +47,11 @@ # The theme to use for HTML and HTML Help pages. See the documentation for # a list of builtin themes. # -html_theme = 'sphinx_rtd_theme' +html_theme = "sphinx_rtd_theme" # Add any paths that contain custom static files (such as style sheets) here, # relative to this directory. They are copied after the builtin static files, # so a file named "default.css" will overwrite the builtin "default.css". -html_static_path = ['_static'] - -autodoc_member_order = 'bysource' +html_static_path = ["_static"] +autodoc_member_order = "bysource" diff --git a/docs/source/index.rst b/docs/source/index.rst index a5cb02f..359e362 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -22,14 +22,14 @@ FRion focuses on time-averaged effects of the ionosphere, for cases where data gets time-averaged before an ionospheric correction can be applied, but also has functions to produce time-series of the ionospheric Faraday rotation. -This package uses `RMExtract `_ +This package uses `RMExtract `_ for the underlying ionospheric calculations. Users interested in alternative ionospheric Faraday rotation packages can look at `ionFR `_, or `ALBUS `_. -A mathematical derivation of how the time-independent ionospheric Faraday -rotation correction is defined, along with some remarks on its use, +A mathematical derivation of how the time-independent ionospheric Faraday +rotation correction is defined, along with some remarks on its use, can be found :download:`here <./Ionospheric_Correction.pdf>`. The package consists of two parts: @@ -40,13 +40,13 @@ The package consists of two parts: Each part can be imported into Python scripts, or the basic functionality can be accessed through the following terminal commands: -- ``frion_predict`` - Runs the time-averaged prediction script, given user-supplied observation - time, location, direction, and frequencies. Can read this information from +- ``frion_predict`` + Runs the time-averaged prediction script, given user-supplied observation + time, location, direction, and frequencies. Can read this information from a FITS header. -- ``frion_correct`` - Runs the correction script, applying a correction (generated by the predict - functions) to a pair of Stokes Q and U cubes to removed the predicted +- ``frion_correct`` + Runs the correction script, applying a correction (generated by the predict + functions) to a pair of Stokes Q and U cubes to removed the predicted ionospheric modulation and depolarization. Use the ``-h`` flag to get detailed usage instructions for each. @@ -62,14 +62,14 @@ Installation FRion has been released on PyPi, so it can be installed with pip as ``pip install FRion``. -Alternatively, it can be installed by downloading the code from -`this link `_, -unzipping, moving the code directory somewhere convenient, +Alternatively, it can be installed by downloading the code from +`this link `_, +unzipping, moving the code directory somewhere convenient, going into the code directory, then running ``pip install -e .``. This will install the package to the Python packages directory. -RMExtract must be installed seperately. It is now available through pip, -using ``pip install RMextract``, but this will try to install casacore as a dependency. +RMExtract must be installed seperately. It is now available through pip, +using ``pip install RMextract``, but this will try to install casacore as a dependency. casacore can be difficult to install on some systems, so my reccomenadation is to install it without casacore by using ``pip install --no-deps RMextract``. casacore is not required for RMExtract: if casacore is not installed, RMExtract will use the pyephem package instead, which installs @@ -82,51 +82,51 @@ be importable using the statements ``import FRion.predict as predict`` and ``frion_predict``, ``frion_timeseries``, and ``frion_correct``. Note that to use the default mode (including the command line tools) requires -an account with CDDIS and a corresponding .netrc file in order to download TEC +an account with CDDIS and a corresponding .netrc file in order to download TEC data; see below for details. -In some cases users may encounter an error +In some cases users may encounter an error ``RuntimeError: Cannot convert due to missing frame information``. This occurs when RMextract finds casacore and tries to use it, but is missing the casadata module. This can be solved by installing casadata (``pip install --index-url https://casa-pip.nrao.edu/repository/pypi-casa-release/simple casadata``), and updating the ``.casarc`` file to point to the updated install. -Removing casacore +Removing casacore (``pip uninstall python-casacore``) solves the issue by forcing RMextract to rely on the ephem module instead. Usage ------------ -To generate ionospheric Faraday rotation predictions as a function of time, +To generate ionospheric Faraday rotation predictions as a function of time, the user can use the predict timeseries functions. These can be used from the command line using ``frion_timeseries``. This requires the user to supply the -start and end times of the observation, the location of the telescope, and the +start and end times of the observation, the location of the telescope, and the sky coordinates. These values can be supplied from a FITS or PSRFITS file if the correct keywords are in the file's header. The user can choose to save the values to a file and/or generate a plot. -These features can be accessed in a Python script using the functions +These features can be accessed in a Python script using the functions available in the :ref:`predict` module, starting with :py:func:`FRion.predict.get_RM()`. Generating time-averaged predictions can be done similarly from the command line using the ``frion_predict`` command or the functions in the :ref:`predict` module. -These predictions require the same information, plus the frequencies of each +These predictions require the same information, plus the frequencies of each channel. -Warnings about reduced accuracy when using PyEphem can be safely ignored +Warnings about reduced accuracy when using PyEphem can be safely ignored (accuracy is approximately 1 arcsecond, and ionospheric data is so coarse that this has no effect on results). Stokes Q and U FITS cubes can be corrected for the time-averaged Faraday rotation, -using the :ref:`correct` module. Tools for time-dependent corrections to -different data types (pulsar observation, visibilities, etc) are outside the +using the :ref:`correct` module. Tools for time-dependent corrections to +different data types (pulsar observation, visibilities, etc) are outside the scope of this package. The correction functions can be used from the command line using ``frion_correct``, or within a script using :py:func:`FRion.correct.apply_correction_to_files()`. The correction relies on the predictions generated by the predict module, so -the user must run the predict tools first and save the ionospheric modulation +the user must run the predict tools first and save the ionospheric modulation to a text file; this text file is used as an input by the correct tools. Note that the correct tools will create a new pair of Stokes Q and U FITS cubes, @@ -134,13 +134,13 @@ with the same size as the input cubes. The user must ensure that sufficient disk space is available, otherwise this step will fail. The default correct functions require holding the full Q and U cubes in RAM -while processing, which may overflow some systems (requiring the use of much +while processing, which may overflow some systems (requiring the use of much slower virtual memory). A large-file version has also been developed to reduce this memory footprint, and can be enabled in the command-line tool by -setting the ``-L`` flag or in scripts by using the +setting the ``-L`` flag or in scripts by using the :py:func:`FRion.correct.apply_correction_large_cube()` function. - - + + .. _ver1.1: @@ -153,7 +153,7 @@ maps from CODE. Unfortunately, CODE stopped producing new TEC maps in January 2023, so it can no longer be used for observations made after that time. The best alternative source of TEC data was the `NASA CDDIS archive`_, -but this requires a (free) account to access the data. At time of writing +but this requires a (free) account to access the data. At time of writing (Aug 2023), RMextract does not support this. A workaround was developed by Anna Ordog and Art Davydov at DRAO that allows RMextract to download from CDDIS, but it is not clear when this will be deployed into the official release. @@ -167,7 +167,7 @@ was selected based on the results of `Porayko et al. 2019`_ who found that the JPL maps produced the lowest residual errors in LOFAR observations. Downloading data from CDDIS requires an account. Information about applying -for an account, as well as creating a .netrc file that carries your login +for an account, as well as creating a .netrc file that carries your login credentials, can be `found here`_. .. _NASA CDDIS archive: https://cddis.nasa.gov/Data_and_Derived_Products/GNSS/atmospheric_products.html diff --git a/setup.py b/setup.py index 339c365..1f59a46 100644 --- a/setup.py +++ b/setup.py @@ -6,26 +6,24 @@ import sys from shutil import rmtree -from setuptools import find_packages, setup, Command +from setuptools import Command, find_packages, setup -NAME = 'FRion' -DESCRIPTION = 'Ionospheric Faraday rotation prediction and correction for radio astronomy polarization cubes.' -URL = 'https://github.com/Cameron-Van-Eck/FRion' -REQUIRES_PYTHON = '>=3.5.0' -VERSION = '1.1.2' -DOWNLOAD_URL = 'https://github.com/Cameron-Van-Eck/FRion/archive/refs/heads/main.zip' +NAME = "FRion" +DESCRIPTION = "Ionospheric Faraday rotation prediction and correction for radio astronomy polarization cubes." +URL = "https://github.com/Cameron-Van-Eck/FRion" +REQUIRES_PYTHON = ">=3.5.0" +VERSION = "1.1.2" +DOWNLOAD_URL = "https://github.com/Cameron-Van-Eck/FRion/archive/refs/heads/main.zip" -REQUIRED = [ - 'numpy', 'astropy', 'pyephem', 'requests' - ] +REQUIRED = ["numpy", "astropy", "pyephem", "requests"] -extras_require={} +extras_require = {} here = os.path.abspath(os.path.dirname(__file__)) try: - with io.open(os.path.join(here, 'README.md'), encoding='utf-8') as f: - long_description = '\n' + f.read() + with io.open(os.path.join(here, "README.md"), encoding="utf-8") as f: + long_description = "\n" + f.read() except FileNotFoundError: long_description = DESCRIPTION @@ -34,28 +32,30 @@ version=VERSION, description=DESCRIPTION, long_description=long_description, - long_description_content_type='text/markdown', + long_description_content_type="text/markdown", python_requires=REQUIRES_PYTHON, url=URL, download_url=DOWNLOAD_URL, - packages=['FRion'], + packages=["FRion"], entry_points={ - 'console_scripts': ['frion_predict=FRion.predict:predict', - 'frion_timeseries=FRion.predict:timeseries', - 'frion_correct=FRion.correct:command_line'], + "console_scripts": [ + "frion_predict=FRion.predict:predict", + "frion_timeseries=FRion.predict:timeseries", + "frion_correct=FRion.correct:command_line", + ], }, install_requires=REQUIRED, include_package_data=True, - license='MIT', + license="MIT", classifiers=[ - 'Development Status :: 5 - Production/Stable', - 'License :: OSI Approved :: MIT License', - 'Programming Language :: Python', - 'Programming Language :: Python :: 3', - 'Programming Language :: Python :: 3.5', - 'Intended Audience :: Science/Research', - 'Topic :: Scientific/Engineering :: Astronomy', + "Development Status :: 5 - Production/Stable", + "License :: OSI Approved :: MIT License", + "Programming Language :: Python", + "Programming Language :: Python :: 3", + "Programming Language :: Python :: 3.5", + "Intended Audience :: Science/Research", + "Topic :: Scientific/Engineering :: Astronomy", ], - maintainer='Cameron Van Eck', - maintainer_email='cameron.vaneck@anu.edu.au', + maintainer="Cameron Van Eck", + maintainer_email="cameron.vaneck@anu.edu.au", )