Merge pull request #11 from radionets-project/sampling_masks

Sampling masks

.github/workflows/ci.yml

@@ -25,6 +25,6 @@ jobs:
           pip install .
 
       - name: tests
-        run: python -m pytest --cov vipy --cov-report=xml
+        run: python -m pytest --cov pyvisgen --cov-report=xml
 
       - uses: codecov/codecov-action@v1

.gitignore

@@ -11,6 +11,7 @@
 # make
 *_done
 */build
+*/*/build
 
 *.ipynb_checkpoints
 *.pyc

README.md

@@ -1,4 +1,4 @@
-# vipy [![Actions Status](https://github.com/radionets-project/vipy/workflows/CI/badge.svg)](https://github.com/radionets-project/vipy/actions) [![codecov](https://codecov.io/gh/radionets-project/vipy/branch/main/graph/badge.svg?token=84ATCQQAMN)](https://codecov.io/gh/radionets-project/vipy)
+# pyvisgen [![Actions Status](https://github.com/radionets-project/pyvisgen/workflows/CI/badge.svg)](https://github.com/radionets-project/pyvisgen/actions) [![codecov](https://codecov.io/gh/radionets-project/pyvisgen/branch/main/graph/badge.svg?token=84ATCQQAMN)](https://codecov.io/gh/radionets-project/pyvisgen)
 
 Possible names:
 * pyvisgen

config/data_set.toml

@@ -0,0 +1,24 @@
+[sampling_options]
+mode = "basic"
+layout = "vlba"
+img_size = 64
+fov_center_ra = [00:00:00.0, 23:59:59.59]
+fov_center_dec = [00:00:00.0, 11:59:59.59]
+fov_size = 0.0064 # max res 0.1
+corr_int_time = 10.0
+scan_start = ["01-01-2020 00:00:01", "31-12-2021 23:59:59"]
+scan_duration = [50, 300]
+scans = [30, 72]
+interval_length = 1200
+base_freq = 15.21e9
+frequsel = [0e8, 0.8e8, 1.44e8, 2.08e8]
+bandwidths = [6.4e7, 6.4e7, 6.4e7, 6.4e7]
+
+[bundle_options]
+num_bundles = 5
+size_bundles = 10
+in_path = "../../../test_radiosim/build/test_data"
+out_path = "./build"
+
+[cluster_options]
+num_jobs = 4

config/default.toml

@@ -3,7 +3,7 @@ fov_center_ra = "12:30:49.423382"
 fov_center_dec = "12:23:28.04366"
 fov_size = 0.00018382
 corr_int_time = 10.0
-scan_start = "2016:95:00:00:00"
+scan_start = "22-05-2021 04:00:01"
 scan_duration = 300
 scans = 72
 interval_length = 1200

examples/baselines.ipynb

@@ -6,7 +6,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "from vipy.simulation.utils import read_config"
+    "from pyvisgen.simulation.utils import read_config"
    ]
   },
   {
@@ -61,7 +61,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "from vipy.layouts.layouts import get_array_layout"
+    "from pyvisgen.layouts.layouts import get_array_layout"
    ]
   },
   {
@@ -101,8 +101,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "from vipy.simulation.scan import get_baselines\n",
-    "from vipy.simulation.utils import calc_time_steps\n",
+    "from pyvisgen.simulation.scan import get_baselines\n",
+    "from pyvisgen.simulation.utils import calc_time_steps\n",
     "import matplotlib.pyplot as plt"
    ]
   },
@@ -215,7 +215,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "from vipy.simulation.scan import create_bgrid"
+    "from pyvisgen.simulation.scan import create_bgrid"
    ]
   },
   {

examples/fits_tests.ipynb

@@ -188,7 +188,7 @@
    "outputs": [],
    "source": [
     "from astropy.time import Time\n",
-    "from vipy.simulation.utils import read_config"
+    "from pyvisgen.simulation.utils import read_config"
    ]
   },
   {
@@ -547,7 +547,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "def create_hdu_list(conf, path=\"../vipy/layouts/eht.txt\"):\n",
+    "def create_hdu_list(conf, path=\"../pyvisgen/layouts/eht.txt\"):\n",
     "    vis_hdu = create_vis_hdu(None, conf)\n",
     "    time_hdu = create_time_hdu(None)\n",
     "    freq_hdu = create_frequency_hdu(conf)\n",

examples/visibility.ipynb

@@ -23,8 +23,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "import vipy.simulation.utils as ut\n",
-    "import vipy.layouts.layouts as layouts\n",
+    "import pyvisgen.simulation.utils as ut\n",
+    "import pyvisgen.layouts.layouts as layouts\n",
     "import astropy.constants as const\n",
     "from astropy import units as un\n",
     "import time as t\n",
@@ -41,7 +41,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "# rc = ut.read_config('/net/nfshome/home/sfroese/vipy/config/default.toml')\n",
+    "# rc = ut.read_config('/net/nfshome/home/sfroese/pyvisgen/config/default.toml')\n",
     "rc = ut.read_config(\"../config/default.toml\")\n",
     "array_layout = layouts.get_array_layout('eht')\n",
     "src_crd = rc['src_coord']\n",
@@ -65,7 +65,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "import vipy.simulation.scan as scan\n",
+    "import pyvisgen.simulation.scan as scan\n",
     "import matplotlib.pyplot as plt"
    ]
   },
@@ -457,7 +457,7 @@
     }
    ],
    "source": [
-    "from vipy.simulation.scan import get_valid_baselines\n",
+    "from pyvisgen.simulation.scan import get_valid_baselines\n",
     "import astropy.units as un\n",
     "from astropy.time import Time\n",
     "visibilities = Visibilities([], [], [], [], [], [], [], [], [], [], [], [])\n",
@@ -545,8 +545,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "import vipy.simulation.utils as ut\n",
-    "import vipy.fits.writer as writer"
+    "import pyvisgen.simulation.utils as ut\n",
+    "import pyvisgen.fits.writer as writer"
    ]
   },
   {

pyvisgen/fits/data.py

@@ -0,0 +1,56 @@
+import numpy as np
+from astropy.io import fits
+from pathlib import Path
+
+
+class fits_data:
+    def __init__(self, data_path):
+        """Handle h5 files simulated with radiosim package.
+
+        Parameters
+        ----------
+        data_path: str
+            path to fits data directory
+
+        """
+        self.path = data_path
+        self.files = self.get_files(Path(data_path))
+        self.num_files = len(self.files)
+
+    def __call__(self):
+        return print("This is the pyvisgen fits files data set class.")
+
+    def __len__(self):
+        """
+        Returns the total number of fits files in this dataset.
+        """
+        return self.num_files
+
+    def __getitem__(self, i):
+        return self.open_file(i)
+
+    def get_files(self, path):
+        return np.sort(np.array([x for x in path.iterdir()]))
+
+    def get_uv_data(self, i):
+        with fits.open(self.files[i]) as hdul:
+            uv_data = hdul[0].data
+            hdul.close()
+            return uv_data
+
+    def get_freq_data(self, i):
+        with fits.open(self.files[i]) as hdul:
+            base_freq = hdul[0].header["CRVAL4"]
+            freq_data = hdul[2].data
+            hdul.close()
+            return freq_data, base_freq
+
+    def open_file(self, i):
+        if isinstance(i, np.ndarray):
+            print("Only one file can be open at the same time!")
+            return
+
+        with fits.open(self.files[i]) as hdul:
+            fits_file = hdul
+            hdul.close()
+            return fits_file.info()

vipy/fits/writer.py → pyvisgen/fits/writer.py

@@ -3,11 +3,12 @@
 from astropy import wcs
 import astropy.units as un
 from astropy.time import Time
-import pandas as pd
 import astropy.constants as const
+import pyvisgen.layouts.layouts as layouts
+import warnings
 
 
-def create_vis_hdu(data, conf, layout="EHT", source_name="sim-source-0"):
+def create_vis_hdu(data, conf, layout="vlba", source_name="sim-source-0"):
     u = data.u
 
     v = data.v
@@ -16,40 +17,45 @@ def create_vis_hdu(data, conf, layout="EHT", source_name="sim-source-0"):
 
     DATE = data.date - int(
         data.date.min()
-    )  # placeholder, julian date of vis, central time in the integration period
+    ) 
 
-    # I think this is not really needed, but dunno, documentation is again insane
-    _DATE = (
-        data._date
-    )  # relative julian date for the observation day??, central time in the integration period
+    _DATE = data._date  # central time in the integration period
 
     BASELINE = data.base_num
 
     INTTIM = np.repeat(np.array(conf["corr_int_time"], dtype=">f4"), len(u))
 
     # visibility data
-    values = data.get_values()
+    values = np.swapaxes(data.get_values(), 0, 1)
+
+    num_ifs = values.shape[1]
+
     vis = np.swapaxes(
         np.swapaxes(
-            np.stack([values.real, values.imag, np.ones(values.shape)], axis=1), 1, 2
-        ),0,1,).reshape(-1, 1, 1, 1, 1, 4, 3)
-    DATA = vis  
+            np.stack([values.real, values.imag, np.ones(values.shape)], axis=2),
+            1,
+            3,
+        ),
+        2,
+        3,
+    ).reshape(-1, 1, 1, num_ifs, 1, 4, 3)
+    DATA = vis
     # in dim 4 = IFs , dim = 1, dim 4 = number of jones, 3 = real, imag, weight
 
     # wcs
-    ra = conf["src_coord"].ra.value
-    dec = conf["src_coord"].dec.value
-    freq, freq_d = freq, freq_d = (
-        (np.array(conf["channel"].split(":")).astype("int") * un.MHz).to(un.Hz).value
-    )
+    ra = conf["fov_center_ra"]
+    dec = conf["fov_center_dec"]
+    freq = (conf["base_freq"] * un.Hz).value
+    freq_d = (conf["bandwidths"][0] * un.Hz).value
+
     ws = wcs.WCS(naxis=7)
     ws.wcs.crpix = [1, 1, 1, 1, 1, 1, 1]
     ws.wcs.cdelt = np.array([1, 1, -1, freq_d, 1, 1, 1])
     ws.wcs.crval = [1, 1, -5, freq, 1, ra, dec]
     ws.wcs.ctype = ["", "COMPLEX", "STOKES", "FREQ", "IF", "RA", "DEC"]
     h = ws.to_header()
 
-    scale = 1  # / freq
+    scale = 1
     u_scale = u / const.c
     v_scale = v / const.c
     w_scale = w / const.c
@@ -88,9 +94,8 @@ def create_vis_hdu(data, conf, layout="EHT", source_name="sim-source-0"):
     hdu_vis.header.comments["PTYPE6"] = "Relative Julian date ?"
     hdu_vis.header.comments["PTYPE7"] = "Integration time"
 
-    date_obs = Time(conf["scan_start"], format="yday").to_value(
-        format="iso", subfmt="date"
-    )
+    date_obs = conf["scan_start"].date().strftime("%Y-%m-%d")
+
     date_map = Time.now().to_value(format="iso", subfmt="date")
 
     # add additional keys
@@ -176,25 +181,22 @@ def create_time_hdu(data):
 
 
 def create_frequency_hdu(conf):
-    freq, freq_d = freq, freq_d = (
-        (np.array(conf["channel"].split(":")).astype("int") * un.MHz).to(un.Hz).value
-    )
-    num_ifs = 1  # at the moment only 1 possible
-
     FRQSEL = np.array([1], dtype=">i4")
     col1 = fits.Column(name="FRQSEL", format="1J", unit=" ", array=FRQSEL)
 
-    IF_FREQ = np.array([[0.00e00]], dtype=">f8")  # start with 0, add ch_with per IF
+    IF_FREQ = np.array(
+        [conf["frequsel"]], dtype=">f8"
+    )  # start with 0, add ch_with per IF
     col2 = fits.Column(
         name="IF FREQ", format=str(IF_FREQ.shape[-1]) + "D", unit="Hz", array=IF_FREQ
     )
 
-    CH_WIDTH = np.repeat(np.array([[freq_d]], dtype=">f4"), 1, axis=1)
+    CH_WIDTH = np.repeat(np.array([conf["bandwidths"]], dtype=">f4"), 1, axis=1)
     col3 = fits.Column(
         name="CH WIDTH", format=str(CH_WIDTH.shape[-1]) + "E", unit="Hz", array=CH_WIDTH
     )
 
-    TOTAL_BANDWIDTH = np.repeat(np.array([[freq_d]], dtype=">f4"), 1, axis=1)
+    TOTAL_BANDWIDTH = np.repeat(np.array([conf["bandwidths"]], dtype=">f4"), 1, axis=1)
     col4 = fits.Column(
         name="TOTAL BANDWIDTH",
         format=str(TOTAL_BANDWIDTH.shape[-1]) + "E",
@@ -232,8 +234,8 @@ def create_frequency_hdu(conf):
     return hdu_freq
 
 
-def create_antenna_hdu(layout_txt, conf, layout="EHT"):
-    array = pd.read_csv(layout_txt, sep=" ")
+def create_antenna_hdu(conf):
+    array = layouts.get_array_layout(conf["layout"], writer=True)
 
     ANNAME = np.chararray(len(array), itemsize=8, unicode=True)
     ANNAME[:] = array["station_name"].values
@@ -300,10 +302,8 @@ def create_antenna_hdu(layout_txt, conf, layout="EHT"):
     )
     hdu_ant = fits.BinTableHDU.from_columns(coldefs_ant)
 
-    freq, freq_d = freq, freq_d = (
-        (np.array(conf["channel"].split(":")).astype("int") * un.MHz).to(un.Hz).value
-    )
-    ref_date = Time(conf["scan_start"], format="yday")
+    freq = (conf["base_freq"] * un.Hz).value
+    ref_date = Time(conf["scan_start"].isoformat(), format="isot")
 
     # add additional keywords
     hdu_ant.header["EXTNAME"] = ("AIPS AN", "AIPS table file")
@@ -335,7 +335,7 @@ def create_antenna_hdu(layout_txt, conf, layout="EHT"):
     hdu_ant.header["UT1UTC"] = (0, "UT1 - UTC (sec)")  # missing
     hdu_ant.header["DATUTC"] = (0, "time system - UTC (sec)")  # missing
     hdu_ant.header["TIMSYS"] = ("UTC", "Time system")
-    hdu_ant.header["ARRNAM"] = (layout, "Array name")
+    hdu_ant.header["ARRNAM"] = (conf["layout"], "Array name")
     hdu_ant.header["XYZHAND"] = ("RIGHT", "Handedness of station coordinates")
     hdu_ant.header["FRAME"] = ("????", "Coordinate frame")
     hdu_ant.header["NUMORB"] = (0, "Number orbital parameters in table (n orb)")
@@ -370,10 +370,11 @@ def create_antenna_hdu(layout_txt, conf, layout="EHT"):
     return hdu_ant
 
 
-def create_hdu_list(data, conf, path="../vipy/layouts/eht.txt"):
+def create_hdu_list(data, conf):
+    warnings.filterwarnings("ignore", module="astropy.io.fits")
     vis_hdu = create_vis_hdu(data, conf)
     time_hdu = create_time_hdu(data)
     freq_hdu = create_frequency_hdu(conf)
-    ant_hdu = create_antenna_hdu(path, conf)
+    ant_hdu = create_antenna_hdu(conf)
     hdu_list = fits.HDUList([vis_hdu, time_hdu, freq_hdu, ant_hdu])
     return hdu_list