Coverage for PyFHD/data_setup/obs.py: 72%
268 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-07-01 10:58 +0800
1import numpy as np
2import logging
3from math import pi, log10
4from PyFHD.pyfhd_tools.pyfhd_utils import (
5 idl_argunique,
6 histogram,
7 angle_difference,
8 parallactic_angle,
9)
10from PyFHD.pyfhd_tools.unit_conv import altaz_to_radec, radec_to_pixel, radec_to_altaz
11from pathlib import Path
12from astropy.io import fits
13from astropy.table import Table
14from astropy.time import Time
17def create_obs(
18 pyfhd_header: dict,
19 params: dict,
20 layout: dict,
21 pyfhd_config: dict,
22 logger: logging.Logger,
23) -> dict:
24 """
25 create_obs takes all the data that has been read in and creates the obs data structure which holds data
26 and metadata of the observation we're doing a PyFHD run on. Inside this function the metafits file will
27 be read as well.
29 Parameters
30 ----------
31 pyfhd_header : dict
32 The data from the UVFITS header
33 params : dict
34 The data from the UVFITS file
35 layout : dict
36 The data dictionary containing data and metadata about the antennas
37 pyfhd_config : dict
38 PyFHD's configuration dictionary
39 logger : logging.Logger
40 The PyFHD logger
42 Returns
43 -------
44 obs : dict
45 The observatiopn data structure for PyFHD containing data from the config and metadata from the observation UVFITS and METAFITS files.
46 """
48 obs = {}
49 baseline_info = {}
51 # Save the data from the header
52 obs["n_pol"] = (
53 pyfhd_config["n_pol"] if pyfhd_config["n_pol"] else pyfhd_header["n_pol"]
54 )
55 obs["n_tile"] = pyfhd_header["n_tile"]
56 obs["n_freq"] = pyfhd_header["n_freq"]
57 obs["n_freq_flag"] = 0
58 obs["instrument"] = pyfhd_config["instrument"]
59 obs["obsname"] = pyfhd_config["obs_id"]
60 # Deal with the times
61 time = params["time"]
62 b0i = idl_argunique(time)
63 obs["n_time"] = b0i.size
64 bin_width = np.empty(obs["n_time"])
65 if obs["n_time"] > 1: 65 ↛ 68line 65 didn't jump to line 68 because the condition on line 65 was always true
66 bin_width[0 : obs["n_time"]] = b0i + 1
67 else:
68 bin_width = time.size
69 b0i_range = np.arange(1, obs["n_time"])
70 bin_width[b0i_range] = b0i[b0i_range] - b0i[b0i_range - 1]
71 baseline_info["bin_offset"] = np.zeros(obs["n_time"], dtype=np.int64)
72 if obs["n_time"] > 1: 72 ↛ 75line 72 didn't jump to line 75 because the condition on line 72 was always true
73 baseline_info["bin_offset"][1:] = np.cumsum(bin_width[: obs["n_time"] - 1])
74 # Deal with the number of visibilities
75 obs["n_baselines"] = int(bin_width[0])
76 obs["n_vis"] = time.size * obs["n_freq"]
77 obs["n_vis_raw"] = obs["n_vis_in"] = obs["n_vis"]
78 obs["nf_vis"] = np.zeros(obs["n_freq"], dtype=np.int64)
80 obs["freq_res"] = pyfhd_header["freq_res"]
81 baseline_info["freq"] = pyfhd_header["frequency_array"]
82 if pyfhd_config["beam_nfreq_avg"] is not None: 82 ↛ 85line 82 didn't jump to line 85 because the condition on line 82 was always true
83 obs["beam_nfreq_avg"] = pyfhd_config["beam_nfreq_avg"]
84 else:
85 obs["beam_nfreq_avg"] = 1
86 freq_bin = obs["beam_nfreq_avg"] * obs["freq_res"]
87 # Let's get the freq_center and the bins we want to use
88 freq_hist, _, freq_ri = histogram(baseline_info["freq"], bin_size=freq_bin)
89 freq_bin_i = np.zeros(obs["n_freq"])
90 for bin in range(freq_hist.size):
91 if freq_ri[bin] < freq_ri[bin + 1]: 91 ↛ 90line 91 didn't jump to line 90 because the condition on line 91 was always true
92 freq_bin_i[freq_ri[freq_ri[bin] : freq_ri[bin + 1]]] = bin
93 baseline_info["fbin_i"] = freq_bin_i.astype(np.int64)
94 obs["freq_center"] = np.median(baseline_info["freq"])
96 antenna_flag = True
97 if np.max(params["antenna1"]) > 0 and (np.max(params["antenna2"]) > 0): 97 ↛ 101line 97 didn't jump to line 101 because the condition on line 97 was always true
98 baseline_info["tile_a"] = params["antenna1"]
99 baseline_info["tile_b"] = params["antenna2"]
100 antenna_flag = False
101 if antenna_flag: 101 ↛ 104line 101 didn't jump to line 104 because the condition on line 101 was never true
102 # 256 tile upper limit is hard-coded in CASA format
103 # these tile numbers have been verified to be correct
104 baseline_min = np.min(params["baseline_arr"])
105 exponent = np.log(np.min(baseline_min)) / np.log(2)
106 antenna_mod_index = 2 ** np.floor(exponent)
107 tile_B_test = baseline_min % antenna_mod_index
108 # Check if a bad fit and if autocorrelations or the first tile are missing
109 if (tile_B_test > 1) and (baseline_min % 2 == 1):
110 antenna_mod_index /= 2 ** np.floor(np.log(np.min(tile_B_test)) / np.log(2))
111 baseline_info["tile_a"] = np.floor(params["baseline_arr"] / antenna_mod_index)
112 baseline_info["tile_b"] = np.fix(params["baseline_arr"] / antenna_mod_index)
113 if (
114 max(np.max(baseline_info["tile_a"]), np.max(baseline_info["tile_b"]))
115 != obs["n_tile"]
116 ):
117 logger.warning(
118 f"Mis-matched n_tiles Header: {obs['n_tile']}, Data: {max(np.max(baseline_info['tile_a']), np.max(baseline_info['tile_b']))}, adjusting n_tiles to be same as data"
119 )
120 obs["n_tile"] = max(
121 np.max(baseline_info["tile_a"]), np.max(baseline_info["tile_b"])
122 )
123 params["antenna1"] = baseline_info["tile_a"]
124 params["antenna2"] = baseline_info["tile_b"]
126 # check that all elements in the antenna1 and antenna2 array exist in the antenna numbers
127 # from the uvfits antenna table
128 all_ants = np.hstack([params["antenna1"], params["antenna2"]])
129 all_ants = np.unique(all_ants)
130 if not (np.all(np.isin(all_ants, layout["antenna_numbers"]))): 130 ↛ 131line 130 didn't jump to line 131 because the condition on line 130 was never true
131 logger.warning("Antenna arrays contain number(s) not found in antenna table")
133 # fhd expects antenna1 and antenna2 arrays containing indices that are one-indexed.
134 # Some uvfits files contain actual antenna numbers in these fields, while others
135 # (particularly, those written by cotter or birli) contain indices.
136 # To account for this, all antenna numbers from the uvfits header are mapped to indices
137 # using the antenna numbers from the uvfits antenna table.
138 # If the antenna numbers were written into the file as indices, they will be mapped to themselves.
139 for tile_i in range(obs["n_tile"]):
140 tile_a_antennas = np.where(
141 layout["antenna_numbers"][tile_i] == params["antenna1"]
142 )
143 if np.size(tile_a_antennas) > 0:
144 baseline_info["tile_a"][tile_a_antennas] = tile_i + 1
145 tile_b_antennas = np.where(
146 layout["antenna_numbers"][tile_i] == params["antenna2"]
147 )
148 if np.size(tile_b_antennas) > 0:
149 baseline_info["tile_b"][tile_b_antennas] = tile_i + 1
150 # Change the type to int to avoid issues with numba
151 baseline_info["tile_a"] = baseline_info["tile_a"].astype(np.int64)
152 baseline_info["tile_b"] = baseline_info["tile_b"].astype(np.int64)
153 params["antenna1"] = baseline_info["tile_a"]
154 params["antenna2"] = baseline_info["tile_b"]
156 baseline_info["freq_use"] = np.ones(obs["n_freq"], dtype=np.int64)
158 # Calculate kx and ky for each baseline at high precision to get most accurate observation information
159 kx_arr = np.outer(baseline_info["freq"], params["uu"])
160 ky_arr = np.outer(baseline_info["freq"], params["vv"])
161 kr_arr = np.sqrt(kx_arr**2 + ky_arr**2)
162 max_baseline = max(np.max(np.abs(kx_arr)), np.max(np.abs(ky_arr)))
164 # Determine the imaging parameters to use
165 if pyfhd_config["FoV"] is not None: 165 ↛ 166line 165 didn't jump to line 166 because the condition on line 165 was never true
166 obs["kpix"] = (180 / pi) / pyfhd_config["FoV"]
167 if pyfhd_config["kbinsize"] is None: 167 ↛ 168line 167 didn't jump to line 168 because the condition on line 167 was never true
168 obs["kpix"] = 0.5
169 else:
170 obs["kpix"] = pyfhd_config["kbinsize"]
172 # Determine observation resolution/extent parameters given number of pixels in x direction (dimension)
173 if pyfhd_config["dimension"] is None and pyfhd_config["elements"] is None: 173 ↛ 174line 173 didn't jump to line 174 because the condition on line 173 was never true
174 obs["dimension"] = 2 ** int(
175 (log10((2 * max_baseline) / pyfhd_config["kpix"]) / log10(2))
176 )
177 obs["elements"] = obs["dimension"]
178 elif pyfhd_config["dimension"] is not None and pyfhd_config["elements"] is None: 178 ↛ 179line 178 didn't jump to line 179 because the condition on line 178 was never true
179 obs["dimension"] = pyfhd_config["dimension"]
180 obs["elements"] = pyfhd_config["dimension"]
181 elif pyfhd_config["dimension"] is None and pyfhd_config["elements"] is not None: 181 ↛ 182line 181 didn't jump to line 182 because the condition on line 181 was never true
182 obs["dimension"] = pyfhd_config["elements"]
183 obs["elements"] = pyfhd_config["elements"]
184 else:
185 obs["dimension"] = pyfhd_config["dimension"]
186 obs["elements"] = pyfhd_config["elements"]
187 # Ensure both dimension and elements are ints to prevent issues down the pipeline
188 obs["dimension"] = int(obs["dimension"])
189 obs["elements"] = int(obs["elements"])
190 obs["degpix"] = (180 / pi) / (obs["kpix"] * pyfhd_config["dimension"])
192 # Set the max and min baseline
193 max_baseline_inds = np.where(
194 (np.abs(kx_arr) / obs["kpix"] < obs["dimension"] / 2)
195 & (np.abs(ky_arr) / obs["kpix"] < obs["elements"] / 2)
196 )
197 obs["max_baseline"] = np.max(np.abs(kr_arr[max_baseline_inds]))
198 if pyfhd_config["min_baseline"] is None: 198 ↛ 199line 198 didn't jump to line 199 because the condition on line 198 was never true
199 obs["min_baseline"] = np.min(kr_arr[np.nonzero(kr_arr)])
200 else:
201 obs["min_baseline"] = max(
202 pyfhd_config["min_baseline"], np.min(kr_arr[np.nonzero(kr_arr)])
203 )
205 meta = read_metafits(obs, pyfhd_header, params, pyfhd_config, logger)
207 baseline_info["time_use"] = np.ones(obs["n_time"], dtype=np.int8)
208 # time cut is specified in seconds to cut (rounded up to next time integration point).
209 # Specify negative time_cut to cut time off the end. Specify a vector to cut at both the start and end
210 if pyfhd_config["time_cut"] is not None: 210 ↛ 211line 210 didn't jump to line 211 because the condition on line 210 was never true
211 time_cut = pyfhd_config["time_cut"]
212 for ti in time_cut:
213 if time_cut[ti] < 0:
214 ti_start = min(
215 obs["n_time"]
216 - max(np.ceil(np.abs(time_cut[ti]) / meta["time_res"]), 0),
217 obs["n_time"] - 1,
218 )
219 ti_end = obs["n_time"] - 1
220 else:
221 ti_start = 0
222 ti_end = min(
223 np.ceil(np.abs(time_cut[ti])) / meta["time_res"] - 1,
224 obs["n_time"] - 1,
225 )
226 if ti_end >= ti_start:
227 baseline_info["time_use"][ti_start : ti_end + 1] = 0
228 obs["n_time_flag"] = obs["n_time"] - np.sum(baseline_info["time_use"])
230 # Flag tiles based on meta data
231 baseline_info["tile_use"] = 1 - meta["tile_flag"]
232 obs["n_tile_flag"] = np.count_nonzero(baseline_info["tile_use"] == 0)
234 # Set the last of obs values
235 if pyfhd_config["dft_threshold"]: 235 ↛ 236line 235 didn't jump to line 236 because the condition on line 235 was never true
236 obs["dft_threshold"] = 1 / (2 * pi) ** 2 * obs["dimension"]
237 else:
238 obs["dft_threshold"] = 0
239 obs["degrid_spectral_terms"] = 0
240 obs["grid_spectral_terms"] = 0
241 obs["alpha"] = -0.8
242 obs["pol_names"] = ["XX", "YY", "XY", "YX", "I", "Q", "U", "V"]
243 obs["residual"] = False
245 # Setup healpix structure for obs
246 healpix = {}
247 healpix["nside"] = 0
248 healpix["n_pix"] = 0
249 # May be none!
250 healpix["ind_list"] = pyfhd_config["healpix_inds"]
251 healpix["n_zero"] = -1
252 obs["healpix"] = healpix
254 # Save the baseline_info into obs
255 baseline_info["jdate"] = meta["jdate"]
256 baseline_info["tile_names"] = meta["tile_names"]
257 baseline_info["tile_height"] = meta["tile_height"]
258 baseline_info["tile_flag"] = meta["tile_flag"]
259 obs["baseline_info"] = baseline_info
261 # Save the last of the metadata into obs
262 for key in meta.keys():
263 if key not in baseline_info.keys():
264 obs[key] = meta[key]
266 return obs
269def read_metafits(
270 obs: dict,
271 pyfhd_header: dict,
272 params: dict,
273 pyfhd_config: dict,
274 logger: logging.Logger,
275) -> dict:
276 """
277 Reads the metafits file provided inside the same input directory as the UVFITS file.
278 It will process the data found in the METAFITS file and then returns a meta dictionary.
279 Which will eventually end up inside the obs dictionary.
281 Parameters
282 ----------
283 obs : dict
284 The current obs structure without the metadata
285 pyfhd_header : dict
286 The data from the UVFITS header
287 params : dict
288 The data from the UVFITS file
289 pyfhd_config : dict
290 PyFHD's configuration dictionary
291 logger : logging.Logger
292 PyFHD's logger
294 Returns
295 -------
296 meta : dict
297 The dictionary holding the metadata from the UVFITS and METAFITS files
298 """
300 meta = {}
301 time = params["time"]
302 b0i = idl_argunique(time)
303 meta["jdate"] = time[
304 b0i
305 ] # Time is already in julian. No need to add the bzero (or pzero) value
306 meta["obsx"] = obs["dimension"] // 2
307 meta["obsy"] = obs["elements"] // 2
308 meta["jd0"] = np.min(meta["jdate"])
309 meta["epoch"] = Time(meta["jd0"], format="jd").to_value("decimalyear")
310 meta_path = Path(pyfhd_config["input_path"], pyfhd_config["obs_id"] + ".metafits")
311 if meta_path.is_file(): 311 ↛ 340line 311 didn't jump to line 340 because the condition on line 311 was always true
312 metadata = fits.open(meta_path)
313 hdr = metadata[0].header
314 data = metadata[1].data
315 # Sort the data by antenna using a stable sort, astropy Table is required to access Antenna column for sorting
316 # Standard Astropy does not do stable sorting, hence use of argsort to do stable sorting
317 data = data[np.array(Table(data).argsort("Antenna", kind="stable"))]
318 single_i = np.where(data["pol"] == data["pol"][0])
319 meta["tile_names"] = data["tile"][single_i]
320 meta["tile_height"] = data["height"][single_i] - pyfhd_header["alt"]
321 meta["tile_flag"] = data["flag"][single_i]
322 if np.sum(meta["tile_flag"]) == meta["tile_flag"].size - 1: 322 ↛ 323line 322 didn't jump to line 323 because the condition on line 322 was never true
323 if pyfhd_config["run_simulation"]:
324 logger.warning("All tiles flagged in metadata")
325 else:
326 logger.error("All tiles flagged in metadata")
327 exit()
328 meta["obsra"] = hdr["RA"]
329 meta["obsdec"] = hdr["DEC"]
330 meta["phasera"] = hdr["RAPHASE"]
331 meta["phasedec"] = hdr["DECPHASE"]
332 meta["time_res"] = hdr["INTTIME"]
333 delays = hdr["DELAYS"].split(",")
334 meta["delays"] = (
335 np.asarray(delays, np.int64)
336 .repeat(obs["n_pol"])
337 .reshape((obs["n_pol"], len(delays)))
338 )
339 else:
340 logger.warning(
341 "METAFITS file has not been found, Calculating obs meta settings from the uvfits header instead"
342 )
343 hdr = None
344 # Simulate the flagging of tiles by taking where tiles don't exist
345 tile_A1 = params["antenna1"]
346 tile_B1 = params["antenna2"]
347 hist_A1, _, _ = histogram(tile_A1, min=1, max=obs["n_tile"])
348 hist_B1, _, _ = histogram(tile_B1, min=1, max=obs["n_tile"])
349 hist_AB = hist_A1 + hist_B1
350 meta["tile_names"] = np.arange(1, obs["n_tile"] + 1)
351 meta["tile_height"] = np.zeros(obs["n_tile"])
352 tile_use = np.where(hist_AB == 0)[0]
353 meta["tile_flag"] = np.zeros(obs["n_tile"], dtype=np.int8)
354 if tile_use.size > 0:
355 meta["tile_flag"][tile_use] = 1
356 if b0i.size > 1:
357 meta["time_res"] = (time[b0i[1]] - time[b0i[0]]) * 24.0 * 3600.0
358 else:
359 meta["time_res"] = 1
360 meta["obsra"] = pyfhd_header["obsra"]
361 meta["obsdec"] = pyfhd_header["obsdec"]
362 meta["phasera"] = pyfhd_header["obsra"]
363 meta["phasedec"] = pyfhd_header["obsdec"]
364 meta["delays"] = None
366 # Store an origin/target phase ra/dec
367 meta["orig_phasera"] = (
368 pyfhd_config["override_target_phasera"]
369 if pyfhd_config["override_target_phasera"] is not None
370 else meta["phasera"]
371 )
372 meta["orig_phasedec"] = (
373 pyfhd_config["override_target_phasedec"]
374 if pyfhd_config["override_target_phasedec"] is not None
375 else meta["phasedec"]
376 )
378 # Get the Zenith RA and DEC from the location and time
379 zenra, zendec = altaz_to_radec(
380 90,
381 0,
382 pyfhd_header["lat"],
383 pyfhd_header["lon"],
384 pyfhd_header["alt"],
385 meta["jd0"],
386 )
387 meta["zenra"] = zenra
388 meta["zendec"] = zendec
390 # Project Slant Orthographic
391 # astr is basically a WCS data structure.
392 # zenx and zeny are the Zenith pixel coordinates
393 meta["astr"], meta["zenx"], meta["zeny"] = project_slant_orthographic(meta, obs)
395 # Get the alt and azimuth of the observation
396 meta["obsalt"], meta["obsaz"] = radec_to_altaz(
397 meta["obsra"],
398 meta["obsdec"],
399 pyfhd_header["lat"],
400 pyfhd_header["lon"],
401 pyfhd_header["alt"],
402 meta["jd0"],
403 )
405 if hdr is not None: 405 ↛ 419line 405 didn't jump to line 419 because the condition on line 405 was always true
406 # Save the raw header and data into the meta dictionary
407 # Save the header as a Python dictionary
408 meta["meta_hdr"] = {}
409 for key in hdr.keys():
410 # Check if they is HISTORY or COMMENT which will be changed to a list for ease of use with hdf5 files
411 if key in ["HISTORY", "COMMENT"]:
412 meta["meta_hdr"][key] = list(hdr[key])
413 else:
414 meta["meta_hdr"][key] = hdr[key]
415 # The astropy FITS_rec class is based off a numpy record array so saving as is should be fine
416 # If so desired a tolist() function to turn the data into list of lists, but you lose the column names
417 meta["meta_data"] = data
419 return meta
422def project_slant_orthographic(meta: dict, obs: dict, epoch=2000) -> dict:
423 """
424 Create an astrometry data structure holding key astrometry information.
425 It's essentially a WCS data structure, done as a Python dictionary allowing greater compatibility
426 with other packages.
428 Parameters
429 ----------
430 meta : dict
431 The current metadata dictionary
432 obs : dict
433 The current obs dictionary
434 epoch : float
435 The equinox used for the dictionary structure
437 Returns
438 -------
439 astr : dict
440 An astrometry structure built from meta and obs
441 """
443 if abs(meta["phasera"] - meta["zenra"]) > 90: 443 ↛ 450line 443 didn't jump to line 450 because the condition on line 443 was always true
444 lon_offset = (
445 meta["phasera"]
446 - (360 if meta["phasera"] > meta["zenra"] else -360)
447 - meta["zenra"]
448 )
449 else:
450 lon_offset = meta["phasera"] - meta["zenra"]
451 zenith_ang = angle_difference(
452 meta["phasera"], meta["phasedec"], meta["zenra"], meta["zendec"], degree=True
453 )
454 parallactic_ang = parallactic_angle(meta["zendec"], lon_offset, meta["phasedec"])
456 xi = -1 * np.tan(np.radians(zenith_ang)) * np.sin(np.radians(parallactic_ang))
457 eta = np.tan(np.radians(zenith_ang)) * np.cos(np.radians(parallactic_ang))
459 # Replicate MAKE_ASTR return dictionary structure from astrolib
460 # We don't have to do it perfectly as it's only used for this function with the above as inputs
461 # This is essentially a WCS in a dictionary for use with other libraries other than Astropy
462 astr = {}
463 astr["naxis"] = np.array([obs["dimension"], obs["elements"]])
464 astr["cd"] = np.identity(2)
465 astr["cdelt"] = np.full(2, obs["degpix"])
466 astr["crpix"] = np.array([meta["obsx"], meta["obsy"]]) + 1
467 astr["crval"] = np.array([meta["phasera"], meta["phasedec"]])
468 projection_name = "SIN"
469 astr["ctype"] = ["RA---" + projection_name, "DEC--" + projection_name]
470 astr["longpole"] = 180
471 astr["latpole"] = 0
472 astr["pv2"] = np.array([xi, eta])
473 # The PV1 array in Astrolib ASTR contains 5 projection parameters associated with longitude axis
474 # [xyoff, phi0, theta0, longpole, latpole]
475 # xyoff and phi0 are 0 as default
476 # The third number [i = 2] is determined by the fact we are using SIN zenithal projections
477 # The last are the longpole and latpole we set earlier
478 astr["pv1"] = np.array([0, 0, 90, 180, 0], dtype=np.float64)
479 astr["axes"] = np.array([1, 2])
480 astr["reverse"] = (
481 0 # Since Axes are always valid Celestial, we don't need to reverse them
482 )
483 astr["coord_sys"] = "C" # Celestial Coordinate System in MAKE_ASTR
484 astr["projection"] = projection_name
485 astr["known"] = np.array([1]) # The projection name is guaranteed to be known
486 astr["radecsys"] = "ICRS" # Using ICRS instead of FK5
487 astr["equinox"] = epoch
488 astr["dateobs"] = Time(meta["jd0"], format="jd").to_value("fits")
489 astr["mjdobs"] = meta["jd0"] - 2400000.5
490 astr["x0y0"] = np.zeros(2, dtype=np.float64)
491 # Get the pixel coordinates of zenra and zendec
492 zenx, zeny = radec_to_pixel(meta["zenra"], meta["zendec"], astr)
494 return astr, zenx, zeny
497def update_obs(
498 obs: dict,
499 dimension: int,
500 kbinsize: float | int,
501 beam_nfreq_avg: float | None = None,
502 fov: float | None = None,
503) -> dict:
504 """
505 Inside the quickview function for exporting files we need to update the obs
506 dictionary based on the new dimension and kbinsize given. This differs slightly from
507 FHD as we only adjust the exact things required for this as we only use this function
508 once in quickview.
510 Parameters
511 ----------
512 obs : dict
513 The original observation dictionary
514 dimension : int
515 The new dimension for the size of each axes
516 kbinsize : float | int
517 The new kbin
518 beam_nfreq_avg: float | None
519 Set the new factor to average up the frequency resolution,by default None
520 fov: float | None
521 Set a new field of view, by default None
523 Returns
524 -------
525 obs: dict
526 The new updated obs dictionary
527 """
528 if beam_nfreq_avg is None: 528 ↛ 532line 528 didn't jump to line 532 because the condition on line 528 was always true
529 beam_nfreq_avg = np.round(
530 obs["n_freq"] / (np.max(obs["baseline_info"]["fbin_i"]) + 1)
531 )
532 freq_bin = beam_nfreq_avg * obs["freq_res"]
533 freq_hist, _, freq_ri = histogram(obs["baseline_info"]["freq"], bin_size=freq_bin)
534 freq_bin_i = np.zeros(obs["n_freq"], dtype=np.int64)
535 for bin in range(freq_hist.size - 1):
536 if freq_ri[bin] < freq_ri[bin + 1]: 536 ↛ 535line 536 didn't jump to line 535 because the condition on line 536 was always true
537 freq_bin_i[freq_ri[freq_ri[bin] : freq_ri[bin + 1]]] = bin
538 # Adjust the obs dictionary based on the new dimension and kbinsize
539 obs["dimension"] = dimension
540 obs["elements"] = dimension
541 obs["obsx"] = dimension / 2
542 obs["obsy"] = dimension / 2
543 obs["kpix"] = kbinsize if fov is None else (180 / np.pi) / fov
544 obs["degpix"] = (180 / np.pi) / (obs["kpix"] * dimension)
545 obs["max_baseline"] = min(
546 obs["max_baseline"], (dimension * obs["kpix"]) / np.sqrt(2)
547 )
548 obs["astr"]["naxis"] = np.array([dimension, dimension])
549 obs["astr"]["cdelt"] = np.full(2, obs["degpix"])
550 obs["astr"]["crpix"] = np.array([obs["obsx"] + 1, obs["obsy"] + 1])
551 obs["baseline_info"]["fbin_i"] = freq_bin_i
553 return obs