Major refactor and add bulk generation

lofar_random_gen.py

+import numpy
+from argparse import ArgumentParser
+import json
+import logging
+from lofardatagenerator.coords import *
+
+logger = logging.getLogger('LOFAR generator')
+logger.setLevel(logging.DEBUG)
+
+__MINUTE_IN_SECONDS = 60
+
+def parse_arguments():
+    parser = ArgumentParser('Generate uvw and compute eccentricity and '
+                            'filling factor of the uvw plane')
+    parser.add_argument('array_coordinates',
+                        type=str,
+                        help='coordinate list of the array elements')
+
+    parser.add_argument('observation_duration',
+                        type=int,
+                        nargs='+',
+                        help='duration of the observation in minutes')
+    parser.add_argument('--n_antennas',
+                        type=int,
+                        default=-1,
+                        help='antenna sample size')
+
+    parser.add_argument('--min_alt',
+                        type=float,
+                        default='10',
+                        help='source_coordinates')
+    parser.add_argument('--n_pointings',
+                        type=float,
+                        default=100)
+
+    parser.add_argument('--random_seed',
+                        default=numpy.random.seed(),
+                        type=int,
+                        help='seed for the random number generator')
+
+    parser.add_argument('--frequency',
+                        default=150.,
+                        type=float,
+                        help='frequency in MHz')
+
+    parser.add_argument('--max_baseline_length',
+                        default=-1,
+                        type=float,
+                        help='max baseline length in m')
+    parser.add_argument('--save_image_to', type=str, default='out')
+    return parser.parse_args()
+
+
+import os
+def main():
+    start_time = datetime.now()
+    arguments = parse_arguments()
+    numpy.random.seed(arguments.random_seed)
+    for i in range(arguments.n_pointings):
+
+        array_x, array_y, array_z = load_station_layout(arguments.array_coordinates,
+                            sample_size = arguments.n_antennas)
+
+        observation_duration_in_seconds = int(numpy.random.choice(arguments.observation_duration) *
+                                              __MINUTE_IN_SECONDS)
+        alt = numpy.random.uniform([arguments.min_alt, 90], 1)[0]
+        az = numpy.random.uniform([0, 360], 1)[0]
+
+        ra, dec = local_coords_to_radians(alt, az)
+
+        uvw = generate_uvw_per_observation(array_x, array_y, array_z, ra, dec,
+                                           start_time,
+                                           observation_duration_in_seconds)
+        uvw = select_baseline_length(uvw, arguments.max_baseline_length)
+        logger.debug("UVW generated")
+        rotated_uvw, pca, r_pca = rotate_ellipse_points(uvw)
+        ecc, filling_factor, sampling, mean_coverage, mean_d = compute_parameters(rotated_uvw)
+        out_path = os.path.join(arguments.save_image_to, f'POINTING_{i}')
+        os.makedirs(out_path, exist_ok=True)
+        make_plots(uvw, rotated_uvw, sampling, arguments.frequency, pca, r_pca, out_path)
+
+        data = dict(
+                    duration=observation_duration_in_seconds,
+                    ra=ra, dec=dec, out_dir=arguments.save_image_to, e=ecc,
+                    filling_factor=filling_factor, average_coverage=mean_coverage,
+                    average_scale=mean_d)
+        with open(os.path.join(arguments.save_image_to, f'parameters_pointing_{i}.json'), 'w') as f_out:
+            json.dump(data, f_out)
+
+
+if __name__== '__main__':
+    main()

lofardatagenerator/coords.py

+from datetime import datetime, timedelta
+import sys
+from astropy.coordinates import SkyCoord
+import pyuvwsim
+import matplotlib.pyplot as plt
+from sklearn.decomposition import PCA
+import scipy
+import numpy
+import logging
+from astropy.coordinates import AltAz, EarthLocation
+import astropy.units as u
+__DEFAULT_LOFAR_ARRAY_CENTER = EarthLocation.from_geocentric(*[3183230.104, 1276816.818, 5359477.705], unit=u.m)
+
+
+def local_coords_to_radians(alt, az, location=__DEFAULT_LOFAR_ARRAY_CENTER):
+    local_frame = AltAz(location=location, obstime=datetime.now())
+    local_coord = SkyCoord(alt=alt * u.deg, az=az * u.deg, frame=local_frame)
+    sky_coord = local_coord.transform_to('icrs')
+    return sky_coord.ra.value, sky_coord.dec.value
+
+
+def coords_into_radians(coords_string):
+    c = SkyCoord(coords_string, unit=(u.hourangle, u.deg))
+    ra, dec = c.ra.radian, c.dec.radian
+    return ra, dec
+
+
+def load_station_layout(file_path: str, sample_size=10, max_baseline_length=-1):
+
+    coords = numpy.loadtxt(file_path, delimiter=',')
+    max_sample = coords.shape[0]
+    if sample_size > 0 and sample_size <= max_sample:
+        print(max_sample, numpy.arange(0, max_sample, 1))
+        sample = numpy.arange(0, max_sample, 1)
+        numpy.random.shuffle(sample)
+        sample = sample[:sample_size]
+        return coords[sample, 0], coords[sample, 1], coords[sample, 2]
+    else:
+        return coords[:, 0], coords[:, 1], coords[:, 2]
+
+
+def get_mjd_from_datetime(time : datetime):
+    mjd = pyuvwsim.datetime_to_mjd(time.year,
+                                    time.month,
+                                    time.day,
+                                    time.hour,
+                                    time.minute,
+                                    time.second)
+    return mjd
+
+
+def generate_uvw_per_observation(array_x, array_y, array_z, ra, dec, start_time,
+                                 observation_duration):
+
+    n_antennas = len(array_x)
+    n_baselines = int(n_antennas * (n_antennas - 1) // 2)
+
+    uvw = numpy.zeros((observation_duration, n_baselines, 3))
+
+    for i in range(observation_duration):
+        mjd = get_mjd_from_datetime(start_time - timedelta(minutes = i))
+        u, v, w = pyuvwsim.evaluate_baseline_uvw(array_x,
+                                                 array_y,
+                                                 array_z,
+                                                 ra, dec, mjd)
+        uvw[i, :, 0] = u
+        uvw[i, :, 1] = v
+        uvw[i, :, 2] = w
+
+    u, v, w = uvw[:, :, 0].flatten(), \
+              uvw[:, :, 1].flatten(), \
+              uvw[:, :, 2].flatten()
+    points_array = numpy.zeros((len(u), 3), dtype=numpy.float32)
+    points_array[:, 0] = u
+    points_array[:, 1] = v
+    points_array[:, 2] = w
+
+    return points_array
+
+
+def rotate_ellipse_points(uvw):
+    __NDIMS = len(uvw.shape)
+
+    pca = PCA(n_components=3)
+
+    mean_uvw = numpy.mean(uvw, axis=0)
+    mean_uvw = uvw - mean_uvw
+    principal_components = pca.fit_transform(mean_uvw.T)
+
+    principal_component_0 = principal_components[0]
+    principal_component_1 = principal_components[1]
+
+    rot_axis = [1., 0, 0]
+    cos_teta = numpy.dot(principal_component_0, rot_axis) / \
+               numpy.linalg.norm(principal_component_0) / \
+               numpy.linalg.norm(rot_axis)
+
+    sin_teta = numpy.sqrt(1 - cos_teta**2)
+
+    rotation_matrix = numpy.zeros((3,3))
+    rotation_matrix[0, 0] = cos_teta
+    rotation_matrix[0, 1] = sin_teta
+    rotation_matrix[1, 0] = -sin_teta
+    rotation_matrix[1, 1] = cos_teta
+    rotation_matrix[2, 2] = 1
+
+    rotated_matrix = numpy.dot(rotation_matrix.T, mean_uvw.T)
+    rotated_components = numpy.dot(rotation_matrix.T, numpy.array(principal_components).T)
+
+    return rotated_matrix, principal_components, rotated_components
+
+
+def compute_parameters(rotated_matrix):
+    u = rotated_matrix[0, :]
+    v = rotated_matrix[1, :]
+    baseline = numpy.sqrt(u ** 2 + v ** 2)
+    grid_step = numpy.max(baseline) / 100
+
+    max_u, min_u = max(rotated_matrix[0, :]), min(rotated_matrix[0, :])
+    max_v, min_v = max(rotated_matrix[1, :]), min(rotated_matrix[1, :])
+    max_w, min_w = max(rotated_matrix[2, :]), min(rotated_matrix[2, :])
+
+    sampling_u = numpy.floor((max_u - min_u) / grid_step)
+    sampling_v = numpy.floor((max_v - min_v) / grid_step)
+    sampling = int(max([sampling_u, sampling_v]))
+
+    dv = max_v - min_v
+    du = max_u - min_u
+    ecc = min([du, dv]) / max([du, dv])
+
+    hist, ax1, ax2 = numpy.histogram2d(rotated_matrix[0, :],
+                                       rotated_matrix[1, :],
+                                       bins=sampling)
+    d = numpy.sqrt(ax1[:-1, numpy.newaxis] ** 2 + ax2[numpy.newaxis, :-1] ** 2)
+
+    mean_coverage = numpy.mean(hist)
+
+    mean_d = numpy.average(d, weights=hist)
+    filling_factor = len(hist[hist > 0]) / (hist.shape[0] * hist.shape[1])
+    return ecc, filling_factor, sampling, mean_coverage, mean_d
+
+
+def plot_lm_coverage(uvw, sampling, frequency_in_mhz):
+    #TODO: THIS might be incorrect. Check
+    __MHZ_IN_HZ = 1.e6
+    __DEG_IN_ARCSEC = 206264.80624709636
+    hist, ax1, ax2 = numpy.histogram2d(uvw[:, 0], uvw[:, 1], bins=sampling)
+    ubl = ax1[1] - ax1[0]
+    vbl = ax2[1] - ax2[0]
+
+    freq = frequency_in_mhz * __MHZ_IN_HZ
+
+    xcoeff = 1./(numpy.pi * ubl / scipy.constants.c * freq)
+    ycoeff = 1./(numpy.pi * vbl / scipy.constants.c * freq)
+
+    values = numpy.zeros_like(hist)
+    values[hist > 0] = 1
+
+    plt.imshow(abs(numpy.fft.fftshift(numpy.fft.fft2(values))),
+               extent=[-xcoeff, xcoeff,
+                       -ycoeff, ycoeff])
+    plt.xlabel('m')
+    plt.ylabel('l')
+
+
+def make_uv_hist_plot(uvw, sampling, pca):
+    hist, ax1, ax2 = numpy.histogram2d(uvw[:, 0], uvw[:, 1], bins=sampling)
+    plt.minorticks_on()
+    plt.imshow(numpy.log10(hist), extent=[
+        ax2[0] * 1.e-6, ax2[-1] * 1.e-6,
+        ax1[0] * 1.e-6, ax1[-1] * 1.e-6
+    ], aspect='equal')
+
+    plt.quiver(0, 0, pca[0][0], pca[0][1], scale=5)
+    plt.quiver(0, 0, pca[1][0], pca[1][1], scale=5)
+    plt.xlabel('x [Mm]')
+    plt.ylabel('y [Mm]')
+
+    plt.colorbar()
+
+def make_baseline_hist_plot(uvw, sampling):
+    u = uvw[:, 0]
+    v = uvw[:, 1]
+    w = uvw[:, 2]
+    baselines = numpy.sqrt(u ** 2 + v ** 2 + w ** 2) / 1.e6
+    plt.hist(baselines, bins=sampling)
+    plt.semilogy()
+    plt.ylabel('N samples')
+    plt.xlabel('$\sqrt{u^2 + v^2 + w^z}$ [Mm]')
+
+
+def make_plots(original_uvw, rotated_uvw, sampling, frequency_in_mhz, pca, r_pca, path):
+    plt.figure('original_uvw')
+    make_uv_hist_plot(original_uvw, sampling, pca)
+    plt.savefig(path + '/original_uvw.png')
+    plt.close()
+
+    plt.figure('lm coverage')
+    plot_lm_coverage(original_uvw, sampling, frequency_in_mhz)
+    plt.savefig(path + '/lm_coverage.png')
+    plt.close()
+
+    plt.figure('baselines histogram')
+    make_baseline_hist_plot(original_uvw, sampling)
+    plt.savefig(path + '/baseline_histogram.png')
+    plt.close()
+
+    plt.figure('rotated_uvw')
+    make_uv_hist_plot(rotated_uvw.T, sampling, r_pca)
+    plt.savefig(path + '/rotated_uvw.png')
+    plt.close()
+
+
+def select_baseline_length(uvw, max_baseline_length):
+    if max_baseline_length > 0:
+        x = uvw[:, 0]
+        y = uvw[:, 1]
+        z = uvw[:, 2]
+        dist = numpy.sqrt(x**2. + y**2. + z**2.)
+        print(dist)
+        uvw_new = uvw[dist <= max_baseline_length, :]
+        if len(uvw_new) == 0:
+            logging.error('max length %s selects an empty sample of antennas',
+                         max_baseline_length)
+            sys.exit(1)
+
+        return uvw_new
+    else:
+        return uvw

lofargen.py

-import pyuvwsim
 import numpy
-from datetime import datetime, timedelta
-import sys
-from astropy.coordinates import SkyCoord
-from astropy import units as u
-import matplotlib.pyplot as plt
-from sklearn.decomposition import PCA
-import scipy
 from argparse import ArgumentParser
-
+import json
 import logging
+from lofardatagenerator.coords import *
 
 logger = logging.getLogger('LOFAR generator')
-logging.basicConfig(level=logging.DEBUG)
+logger.setLevel(logging.DEBUG)
 
+__MINUTE_IN_SECONDS = 60
 
 def parse_arguments():
     parser = ArgumentParser('Generate uvw and compute eccentricity and '
@@ -49,195 +43,10 @@ def parse_arguments():
                         default=-1,
                         type=float,
                         help='max baseline length in m')
-
+    parser.add_argument('--save_image_to', type=str, default='.')
     return parser.parse_args()
 
 
-def coords_into_radians(coords_string):
-    c = SkyCoord(coords_string, unit=(u.hourangle, u.deg))
-    ra, dec = c.ra.radian, c.dec.radian
-    return ra, dec
-
-
-def load_station_layout(file_path: str, sample_size=10, max_baseline_length=-1):
-
-    coords = numpy.loadtxt(file_path, delimiter=',')
-    max_sample = coords.shape[0]
-    if sample_size > 0 and sample_size <= max_sample:
-        print(max_sample, numpy.arange(0, max_sample, 1))
-        sample = numpy.arange(0, max_sample, 1)
-        numpy.random.shuffle(sample)
-        sample = sample[:sample_size]
-        return coords[sample, 0], coords[sample, 1], coords[sample, 2]
-    else:
-        return coords[:, 0], coords[:, 1], coords[:, 2]
-
-
-def get_mjd_from_datetime(time : datetime):
-    mjd = pyuvwsim.datetime_to_mjd(time.year,
-                                    time.month,
-                                    time.day,
-                                    time.hour,
-                                    time.minute,
-                                    time.second)
-    return mjd
-
-
-def generate_uvw_per_observation(array_x, array_y, array_z, ra, dec, start_time,
-                                 observation_duration):
-
-    n_antennas = len(array_x)
-    n_baselines = int(n_antennas * (n_antennas - 1) // 2)
-
-    uvw = numpy.zeros((observation_duration, n_baselines, 3))
-
-    for i in range(observation_duration):
-        mjd = get_mjd_from_datetime(start_time - timedelta(minutes = i))
-        u, v, w = pyuvwsim.evaluate_baseline_uvw(array_x,
-                                                 array_y,
-                                                 array_z,
-                                                 ra, dec, mjd)
-        uvw[i, :, 0] = u
-        uvw[i, :, 1] = v
-        uvw[i, :, 2] = w
-
-    u, v, w = uvw[:, :, 0].flatten(), \
-              uvw[:, :, 1].flatten(), \
-              uvw[:, :, 2].flatten()
-    points_array = numpy.zeros((len(u), 3), dtype=numpy.float32)
-    points_array[:, 0] = u
-    points_array[:, 1] = v
-    points_array[:, 2] = w
-
-    return points_array
-
-
-def rotate_ellipse_points(uvw):
-    __NDIMS = len(uvw.shape)
-
-    pca = PCA(n_components=3)
-    principal_components = pca.fit_transform(uvw.T)
-    logger.debug("Found PCA")
-    logger.debug("PCAs %s", principal_components)
-    principal_component_0 = principal_components[0]
-    principal_component_1 = principal_components[1]
-
-    cos_teta = numpy.dot(principal_component_0, principal_component_1) / \
-               numpy.linalg.norm(principal_component_0) / \
-               numpy.linalg.norm(principal_component_1)
-
-    sin_teta = numpy.sqrt(1 - cos_teta**2)
-
-    rotation_matrix = numpy.zeros_like(principal_components)
-    rotation_matrix[0, 0] = cos_teta
-    rotation_matrix[0, 1] = sin_teta
-    rotation_matrix[1, 0] = -sin_teta
-    rotation_matrix[1, 1] = cos_teta
-    rotation_matrix[2, 2] = 1
-
-    rotated_matrix = numpy.dot(rotation_matrix, uvw.T)
-    return rotated_matrix
-
-
-def compute_parameters(rotated_matrix):
-    u = rotated_matrix[0, :]
-    v = rotated_matrix[1, :]
-    baseline = numpy.sqrt(u ** 2 + v ** 2)
-    grid_step = numpy.max(baseline) / 100
-
-    max_u, min_u = max(rotated_matrix[0, :]), min(rotated_matrix[0, :])
-    max_v, min_v = max(rotated_matrix[1, :]), min(rotated_matrix[1, :])
-    max_w, min_w = max(rotated_matrix[2, :]), min(rotated_matrix[2, :])
-
-    sampling_u = numpy.floor((max_u - min_u) / grid_step)
-    sampling_v = numpy.floor((max_v - min_v) / grid_step)
-    sampling = int(max([sampling_u, sampling_v]))
-
-    dv = max_v - min_v
-    du = max_u - min_u
-    ecc = min([du, dv]) / max([du, dv])
-
-    hist, ax1, ax2 = numpy.histogram2d(rotated_matrix[0, :],
-                                       rotated_matrix[1, :],
-                                       bins=sampling)
-
-    filling_factor = len(hist[hist > 0]) / (hist.shape[0] * hist.shape[1])
-    return ecc, filling_factor, sampling
-
-
-def plot_lm_coverage(uvw, sampling, frequency_in_mhz):
-    #TODO: THIS might be incorrect. Check
-    __MHZ_IN_HZ = 1.e6
-    __DEG_IN_ARCSEC = 206264.80624709636
-    hist, ax1, ax2 = numpy.histogram2d(uvw[:, 0], uvw[:, 1], bins=sampling)
-    ubl = ax1[1] - ax1[0]
-    vbl = ax2[1] - ax2[0]
-
-    freq = frequency_in_mhz * __MHZ_IN_HZ
-
-    xcoeff = 1./(numpy.pi * ubl / scipy.constants.c * freq)
-    ycoeff = 1./(numpy.pi * vbl / scipy.constants.c * freq)
-
-    values = numpy.zeros_like(hist)
-    values[hist > 0] = 1
-
-    plt.imshow(abs(numpy.fft.fftshift(numpy.fft.fft2(values))),
-               extent=[-xcoeff, xcoeff,
-                       -ycoeff, ycoeff])
-    plt.xlabel('m')
-    plt.ylabel('l')
-
-
-def make_uv_hist_plot(uvw, sampling):
-    hist, ax1, ax2 = numpy.histogram2d(uvw[:, 0], uvw[:, 1], bins=sampling)
-    plt.imshow(numpy.log10(hist), extent=[
-        ax2[0], ax2[-1],
-        ax1[0], ax1[-1]
-    ])
-    plt.xlabel('x [km]')
-    plt.ylabel('y [km]')
-
-    plt.colorbar()
-
-def make_baseline_hist_plot(uvw, sampling):
-    u = uvw[:, 0]
-    v = uvw[:, 1]
-    w = uvw[:, 2]
-    baselines = numpy.sqrt(u ** 2 + v ** 2 + w ** 2)
-    plt.hist(baselines, bins=sampling)
-    plt.xlabel('baseline length [km]')
-
-def make_plots(original_uvw, rotated_uvw, sampling, frequency_in_mhz):
-    plt.figure('original_uvw')
-    make_uv_hist_plot(original_uvw, sampling)
-    plt.figure('lm coverage')
-    plot_lm_coverage(original_uvw, sampling, frequency_in_mhz)
-    plt.figure('baselines histogram')
-    make_baseline_hist_plot(original_uvw, sampling)
-    plt.figure('rotated_uvw')
-    make_uv_hist_plot(rotated_uvw.T, sampling)
-    plt.show()
-
-__MINUTE_IN_SECONDS = 60
-
-
-def select_baseline_length(uvw, max_baseline_length):
-    if max_baseline_length > 0:
-        x = uvw[:, 0]
-        y = uvw[:, 1]
-        z = uvw[:, 2]
-        dist = numpy.sqrt(x**2. + y**2. + z**2.)
-        print(dist)
-        uvw_new = uvw[dist <= max_baseline_length, :]
-        if len(uvw_new) == 0:
-            logger.error('max length %s selects an empty sample of antennas',
-                         max_baseline_length)
-            sys.exit(1)
-
-        return uvw_new
-    else:
-        return uvw
-
 
 def main():
     start_time = datetime.now()
@@ -255,14 +64,18 @@ def main():
                                        observation_duration_in_seconds)
     uvw = select_baseline_length(uvw, arguments.max_baseline_length)
     logger.debug("UVW generated")
-    rotated_uvw = rotate_ellipse_points(uvw)
-    ecc, filling_factor, sampling = compute_parameters(rotated_uvw)
+    rotated_uvw, pca, r_pca = rotate_ellipse_points(uvw)
+    ecc, filling_factor, sampling, mean_coverage, mean_d = compute_parameters(rotated_uvw)
 
     if(arguments.make_uvw_plot):
-        make_plots(uvw, rotated_uvw, sampling, arguments.frequency)
-
-    print("Eccentricity: ", ecc)
-    print("Filling factor: ", filling_factor)
+        make_plots(uvw, rotated_uvw, sampling, arguments.frequency, pca, r_pca, arguments.save_image_to)
+
+    data = dict(frequency=arguments.frequency, duration=observation_duration_in_seconds,
+                ra=ra, dec=dec, out_dir=arguments.save_image_to, e=ecc,
+                filling_factor=filling_factor, average_coverage=mean_coverage,
+                average_scale=mean_d)
+    with open(arguments.save_image_to + '/parameters.json', 'w') as f_out:
+        json.dump(data, f_out)
 
 
 if __name__== '__main__':