Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add report on daytime proportion of the observation #3

Merged
merged 2 commits into from
Dec 13, 2023
Merged

Add report on daytime proportion of the observation #3

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
b2b9bee
Parse args before logging start time
jd-au Feb 25, 2022
c6b9d4d
Add solar elevation, daytime fraction, solar separation and target el…
jd-au Dec 13, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
163 changes: 144 additions & 19 deletions gaskap-hi-validation.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -26,11 +26,12 @@

import aplpy
from astropy.constants import k_B
from astropy.coordinates import SkyCoord
from astropy.coordinates import solar_system_ephemeris, EarthLocation, Angle, SkyCoord, AltAz, FK5, get_body
from astropy.io import ascii, fits
from astropy.io.votable import parse, from_table, writeto
from astropy.io.votable.tree import Info
from astropy.table import Table, Column
from astropy.time import Time, TimezoneInfo
import astropy.units as u
from astropy.utils.exceptions import AstropyWarning
from astropy.wcs import WCS
Expand Down Expand Up @@ -172,7 +173,7 @@ def convert_slab_to_jy(slab, header):
elif 'RESTFRQ' in header.keys():
restfreq = header['RESTFRQ']*u.Hz

if slab.unmasked_data[0,0,0].unit != u.Jy:
if slab.unmasked_data[0,0,0].unit != u.Jy and slab.unmasked_data[0,0,0].unit != u.Jy/u.beam:
print ("Converting slab from {} to Jy".format(slab.unmasked_data[0,0,0].unit) )
print (slab)
slab.allow_huge_operations=True
Expand Down Expand Up @@ -300,7 +301,78 @@ def calc_velocity_res(hdr):
return spec_res_km_s


def report_observation(image, reporter, input_duration, sched_info, obs_metadata):
def get_observing_location():
# Coordinates for Inyarrimanha Ilgari Bundara, the CSIRO Murchison Radio-astronomy Observatory
latitude = Angle("-26:41:46.0", unit=u.deg)
longitude = Angle("116:38:13.0", unit=u.deg)
return EarthLocation(lat=latitude, lon=longitude)


def plot_solar_elevation(fig_folder, obs_time_utc, duration, field_pos, sbid):
obs_end_utc = obs_time_utc + duration
tmjd = np.arange(obs_time_utc.mjd, obs_end_utc.mjd, 1./24.0/60.0/6.0)
times = Time(tmjd, format='mjd', scale='utc')

observing_location = get_observing_location()

solar_system_ephemeris.set('builtin')
sun_sc = get_body("sun", times, location=observing_location)
#sun_sc_median = get_body("sun", Time(np.median(tmjd), format='mjd', scale='utc'), location=observing_location)

sun_altaz = sun_sc.transform_to(AltAz(obstime=times, location=observing_location))
sun_altdeg = sun_altaz.alt.deg
field_sun_sep_deg = sun_sc.separation(field_pos).degree

sns.set()
fig = plt.figure()
ax1 = fig.add_subplot(111)
plot, = ax1.plot((tmjd-tmjd[0])*24.0, sun_altdeg, marker='None', color="black")
#print("Solar position = %s" %(sun_sc_median.to_string(style='hmsdms')))
ax1.set_title("Solar elevation for SBID %s" %(sbid))
ax1.set_xlabel("Time since start of observation (h)")
ax1.set_ylabel("Elevation (deg)")
ax1.axhline(0, color='g')
ax1.axhline(-6, color='blue')
ax1.axhline(-12, color='darkblue')
ax1.axhline(-18, color='indigo')
fig_path = fig_folder + 'solar_el.png'
fig.savefig(fig_path)
plt.close()

# Calculate the daylight percent of the observation
sun_up_filter = sun_altdeg >= 0
steps_sun_up = np.sum(sun_up_filter)
pct_daylight = steps_sun_up/len(sun_altdeg)*100


return fig_path, pct_daylight, field_sun_sep_deg


def plot_target_elevation(fig_folder, obs_time_utc, duration, field_pos, sbid):
obs_end_utc = obs_time_utc + duration
tmjd = np.arange(obs_time_utc.mjd, obs_end_utc.mjd, 1./24.0/60.0/6.0)
times = Time(tmjd, format='mjd', scale='utc')

observing_location = get_observing_location()

field_altaz = field_pos.transform_to(AltAz(obstime=times, location=observing_location))
field_altdeg = field_altaz.alt.deg

fig = plt.figure()
ax1 = fig.add_subplot(111)
plot, = ax1.plot((tmjd-tmjd[0])*24.0, field_altdeg, marker='None', color="black")
ax1.set_title("Target elevation for SBID %s" %(sbid))
ax1.set_xlabel("Time since start of observation (h)")
ax1.set_ylabel("Elevation (deg)")
ax1.axhline(15, color='g')
fig_path = fig_folder + 'target_el.png'
fig.savefig(fig_path)
plt.close()

return fig_path


def report_observation(image, dest_folder, reporter, input_duration, sched_info, obs_metadata):
print('\nReporting observation based on ' + image)

hdr = fits.getheader(image)
Expand All @@ -312,18 +384,20 @@ def report_observation(image, reporter, input_duration, sched_info, obs_metadata
if project.startswith('AS'):
proj_link = "https://confluence.csiro.au/display/askapsst/{0}+Data".format(project)

date = hdr['DATE-OBS']
obs_date = hdr['DATE-OBS'] if 'DATE-OBS' in hdr else 'Unknown'
duration = float(hdr['DURATION'])/3600 if 'DURATION' in hdr else input_duration

naxis1 = int(hdr['NAXIS1'])
naxis2 = int(hdr['NAXIS2'])
pixcrd = np.array([[naxis1/2, naxis2/2]])
centre = w.all_pix2world(pixcrd,1)
centre = SkyCoord(ra=centre[0][0], dec=centre[0][1], unit="deg,deg").to_string(style='hmsdms',sep=':')
centre_coord = SkyCoord(ra=centre[0][0], dec=centre[0][1], unit="deg,deg")
centre = centre_coord.to_string(style='hmsdms',sep=':', precision=1)

# spectral axis
spectral_unit = 'None'
spectral_range = ''
spec_title = 'Spectral Range'
for i in range(3,int(hdr['NAXIS'])+1):
ctype = hdr['CTYPE'+str(i)]
if (ctype.startswith('VEL') or ctype.startswith('VRAD') or ctype.startswith('FREQ')):
Expand Down Expand Up @@ -361,17 +435,41 @@ def report_observation(image, reporter, input_duration, sched_info, obs_metadata

footprint = sched_info.footprint
if footprint and sched_info.pitch:
footprint = "{}_{}".format(footprint, sched_info.pitch)
footprint = "{}_{}".format(footprint, sched_info.pitch)

section = ReportSection('Observation')
section.add_item('SBID', value=sbid)
section.add_item('Project', value=project, link=proj_link)
section.add_item('Date', value=date)
section.add_item('Date (UTC)', value=obs_date)
section.add_item('Duration<br/>(hours)', value='{:.2f}'.format(duration))
section.add_item('Field(s)', value=field_names)
section.add_item('Field Centre(s)', value=field_centres)
section.add_item('Correlator<br/>Mode', value=sched_info.corr_mode)
section.add_item('Footprint', value=footprint)
section.add_item('{}<br/>({})'.format(spec_title, spectral_unit), value=spectral_range)

if obs_date != 'Unknown':
murchison_zone = TimezoneInfo(utc_offset=8*u.hour)
obs_time_utc = Time(obs_date, format='isot', scale='utc')
local_start = obs_time_utc.to_datetime(timezone=murchison_zone)
local_end = (obs_time_utc + duration*u.hour).to_datetime(timezone=murchison_zone)
local_time = '{:%Y-%m-%d %H:%M:%S%z}<br/>{:%Y-%m-%d %H:%M:%S%z}'.format(local_start, local_end)

fig_folder= get_figures_folder(dest_folder)
field_pos = SkyCoord(obs_metadata.fields[0].ra, obs_metadata.fields[0].dec, frame=FK5, unit=(u.hourangle, u.deg)) if obs_metadata else centre_coord
solr_el_path, pct_daylight, field_sun_sep_deg = plot_solar_elevation(fig_folder, obs_time_utc, duration*u.hour, field_pos, sbid)

tgt_el_path = plot_target_elevation(fig_folder, obs_time_utc, duration*u.hour, field_pos, sbid)

section.start_new_row()
section.add_item('Local Time', value=local_time)
#section.add_item('Solar Elevation', link=solr_el_path_rel, image=solr_el_path_rel) # Second should be a thumbnail
add_opt_image_section('Solar Elevation', solr_el_path, fig_folder, dest_folder, section)
section.add_item('Daytime Proportion (%)', value="{:.1f}".format(pct_daylight))
section.add_item('Solar Separation (deg)', value="{:.3f} to {:.3f}".format(field_sun_sep_deg[0], field_sun_sep_deg[1]))
add_opt_image_section('Target Elevation', tgt_el_path, fig_folder, dest_folder, section)


reporter.add_section(section)

reporter.project = project
Expand All @@ -387,7 +485,7 @@ def report_cube_stats(cube, reporter):
# Cube information
askapSoftVer = 'N/A'
askapPipelineVer = 'N/A'
history = hdr['history']
history = hdr['history'] if 'history' in hdr else []
askapSoftVerPrefix = 'Produced with ASKAPsoft version '
askapPipelinePrefix = 'Processed with ASKAP pipeline version '
for row in history:
Expand All @@ -401,6 +499,10 @@ def report_cube_stats(cube, reporter):
beam_min = hdr['BMIN'] * 60 * 60
beam = '{:.1f} x {:.1f}'.format(beam_maj, beam_min)

units = 'N/A'
if 'BUNIT' in hdr:
units = hdr['BUNIT']

dims = []
for i in range(1,int(hdr['NAXIS'])+1):
dims.append(str(hdr['NAXIS'+str(i)]))
Expand All @@ -415,6 +517,7 @@ def report_cube_stats(cube, reporter):
section.add_item('Synthesised Beam<br/>(arcsec)', value=beam)
section.add_item('Sky Area<br/>(deg2)', value='')
section.add_item('Dimensions', value=dimensions)
section.add_item('Units', value=units)
reporter.add_section(section)
return

Expand Down Expand Up @@ -569,7 +672,7 @@ def measure_spectral_line_noise(slab, cube, vel_start, vel_end, reporter, dest_f
# Scale the noise to mJy / 5 kHz channel
std_data = np.nanstd(slab.unmasked_data[:], axis=0)
noise_5kHz = std_data / np.sqrt(1 / spec_res_km_s)
noise_5kHz = noise_5kHz.to(u.mJy) # Jy => mJy
noise_5kHz = noise_5kHz*1000 # Jy => mJy

# Extract the spectral line noise map
mom0_prefix = build_fname(cube, '_mom0_off')
Expand All @@ -581,7 +684,7 @@ def measure_spectral_line_noise(slab, cube, vel_start, vel_end, reporter, dest_f

# Produce the noise plots
cube_name = os.path.basename(cube)
plot_map(folder+prefix, "Spectral axis noise map for " + cube_name, cmap='mako_r', stretch='arcsinh',
plot_map(folder+prefix, "Spectral axis noise map for " + cube_name, cmap='mako', stretch='arcsinh',
colorbar_label=r'Noise level per 5 kHz channel (mJy beam$^{-1}$)')
plot_histogram(folder+prefix, r'Noise level per 5 kHz channel (mJy beam$^{-1}$)', 'Spectral axis noise for ' + cube_name)
median_noise_5kHz = np.nanmedian(noise_5kHz.value[noise_5kHz.value!=0.0])
Expand Down Expand Up @@ -705,7 +808,7 @@ def report_image_stats(image, noise_file, reporter, dest_folder, diagnostics_dir
#img_flux = np.sum(img_data[~np.isnan(img_data)]) / (1.133*((beam_maj * beam_min) / (img.raPS * img.decPS))) #divide by beam area

# Copy pipleine plots
field_src_plot = copy_existing_image(diagnostics_dir+'/image.i.SB*.cont.restored_sources.png', fig_folder)
field_src_plot = copy_existing_image(diagnostics_dir+'/image.i.SB*.cont.restored_sources.png', fig_folder) if diagnostics_dir else None

image_name = os.path.basename(image)
section = ReportSection('Image', image_name)
Expand All @@ -728,7 +831,7 @@ def set_velocity_range(emvelstr, nonemvelstr):
raise ValueError('Velocity {} is not one of the supported GASS velocity steps e.g. 165, 200.'.format(emvel))
nonemvel = int(nonemvelstr)
if not nonemvel in vel_steps:
raise ValueError('Velocity {} is not one of the supported GASS velocity steps e.g. 165, 200.'.format(emvel))
raise ValueError('Velocity {} is not one of the supported GASS velocity steps e.g. 165, 200.'.format(nonemvel))

idx = vel_steps.index(emvel)
if idx +1 >= len(vel_steps):
Expand Down Expand Up @@ -961,9 +1064,12 @@ def extract_spectra(cube, source_cat, dest_folder, reporter, num_spectra, beam_l
unit = slab.unit

for j, pos in enumerate(pix_pos_bright):
data = slab[:,pos[1], pos[0]]
#data = convert_data_to_jy(data, header)
spectra_bright[j][i:max_idx] = data.value
if 0 <= pos[1] < slab.shape[1] and 0 <= pos[0] < slab.shape[2]:
data = slab[:,pos[1], pos[0]]
#data = convert_data_to_jy(data, header)
spectra_bright[j][i:max_idx] = data.value
else:
print ("Skipping src {} as y:{} x:{} is outside cube".format(j, pos[1], pos[0]))
for j, pos in enumerate(pix_pos_beam):
data = slab[:,pos[1], pos[0]]
spectra_beam[j][i:max_idx] = data.value
Expand Down Expand Up @@ -1056,7 +1162,7 @@ def add_opt_image_section(title, image_path, fig_folder, dest_folder, section, t
section.add_item(title, value='N/A')
return

img_thumb, img_thumb_rel = Diagnostics.make_thumbnail(image_path, fig_folder, dest_folder, size_x=thumb_size_y,
img_thumb, img_thumb_rel = Diagnostics.make_thumbnail(image_path, fig_folder, dest_folder, size_x=thumb_size_x,
size_y=thumb_size_y)
image_path_rel = os.path.relpath(image_path, dest_folder)
section.add_item(title, link=image_path_rel, image=img_thumb_rel)
Expand Down Expand Up @@ -1243,22 +1349,41 @@ def report_self_cal(cube, image, obs_metadata, dest_folder, reporter):
reporter.add_metric(metric)


def log_config(args, dest_folder, figures_folder):
print ('Configuration:')
print (' {: >20} : {}'.format('cube', args.cube))
print (' {: >20} : {}'.format('image', args.image))
print (' {: >20} : {}'.format('source_cat', args.source_cat))
print (' {: >20} : {}'.format('beam_list', args.beam_list))
print (' {: >20} : {}'.format('duration', args.duration))
print (' {: >20} : {}'.format('emvel', args.emvel))
print (' {: >20} : {}'.format('nonemvel', args.nonemvel))
print (' {: >20} : {}'.format('noise', args.noise))
print (' {: >20} : {}'.format('redo', args.redo))
print (' {: >20} : {}'.format('num_spectra', args.num_spectra))
print (' {: >20} : {}'.format('dest_folder', dest_folder))
print (' {: >20} : {}'.format('figures_folder', figures_folder))
print ('')


def main():
# Parse command line options
args = parseargs()

start = time.time()
print("#### Started validation at {} ####".format(
(time.strftime('%Y-%m-%d %H:%M:%S', time.localtime(start)))))

#ignore astropy warnings
warnings.simplefilter('ignore', AstropyWarning)

# Parse command line options
args = parseargs()
dest_folder = args.output
if not os.path.exists(dest_folder):
os.makedirs(dest_folder)
figures_folder = dest_folder + '/figures'
if not os.path.exists(figures_folder):
os.makedirs(figures_folder)
log_config(args, dest_folder, figures_folder)

if args.cube and (not os.path.exists(args.cube) or not os.path.isfile(args.cube)):
raise ValueError('Cube {} could not be found or is not a file.'.format(args.cube))
Expand Down Expand Up @@ -1287,7 +1412,7 @@ def main():
sched_info = Diagnostics.get_sched_info(obs_img)
diagnostics_dir = Diagnostics.find_diagnostics_dir(args.cube, args.image)
obs_metadata = Diagnostics.get_metadata(diagnostics_dir) if diagnostics_dir else None
sbid = report_observation(obs_img, reporter, args.duration, sched_info, obs_metadata)
sbid = report_observation(obs_img, dest_folder, reporter, args.duration, sched_info, obs_metadata)

if args.cube:
report_cube_stats(args.cube, reporter)
Expand Down
Loading