diff --git a/lofar_random_gen.py b/lofar_random_gen.py
new file mode 100644
index 0000000000000000000000000000000000000000..313509118bad04adef8e10ea8ae496f8cd1f7735
--- /dev/null
+++ b/lofar_random_gen.py
@@ -0,0 +1,92 @@
+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()
diff --git a/lofardatagenerator/__init__.py b/lofardatagenerator/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/lofardatagenerator/coords.py b/lofardatagenerator/coords.py
new file mode 100644
index 0000000000000000000000000000000000000000..30f8e670de08076ea9c65374076d8481156d219e
--- /dev/null
+++ b/lofardatagenerator/coords.py
@@ -0,0 +1,229 @@
+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
diff --git a/lofargen.py b/lofargen.py
index 044fd6f7ab40fb17df072934a1111bfd77fe5116..033156af5d4b36862dc1d16463b54036d6d8c895 100644
--- a/lofargen.py
+++ b/lofargen.py
@@ -1,19 +1,13 @@
-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__':