# -*- coding: utf-8 -*- """ This Python program implements the algorithm used to determine the average polar- ization properties of SKR. In contrast to the get dataframe function in the previous program, which merges only the N2 and N3d datasets for spectrogram construction, the get dataframe function here additionally incorporates the beta-angle file, thereby enabling the application of geometry-dependent selection criteria. On this basis, the program filters the data and computes hourly statistical measures of the linear, circular, and total polarization used in the subsequent analysis. @author: Dorian Jost """ import readafile from concurrent.futures import ProcessPoolExecutor from itertools import repeat import pandas as pd import os import numpy as np def get_dataframe(n2_file,n3d_file, beta_file): try: n2 = readafile.read_n2(n2_file) n3d = readafile.read_n3d(n3d_file) beta_csv = pd.read_csv(beta_file, sep = '\s+', skiprows = 2) print(f"Loaded {n2_file}, {n3d_file} and {beta_file}") except FileNotFoundError: print("Missing file") #turn n2 into pandas dataframe n2_df = pd.DataFrame.from_records(n2) #turn n3d into pandas dataframe snr_data = np.vstack(n3d['sn']) snr_df = pd.DataFrame.from_records(snr_data,columns = ['sn0','sn1']) names = n3d.dtype.names keep = names[:-1] n3d_df_ = pd.DataFrame.from_records({ n: n3d[n] for n in keep })#, axis = 1) # ignore snr data to avoid collision n3d_df = pd.concat([n3d_df_,snr_df],axis = 1) #merge back together with two new columns sn0 and sn1 DF = pd.merge(n2_df, n3d_df, on='num', how='inner') # to merge DF that contains n2 and n3d data with beta file we need to normalize t97 value to 10 digits after comma eg: 3193.2506047453703 -> 3193.2506047454, 3193.250234375 -> 3193.2502343750 t97 = DF['t97'] t97_padded = t97.map(lambda x: f"{x:.10f}") DF['t97'] = t97_padded t97_beta = beta_csv['t97'] t97_padded = t97_beta.map(lambda x: f"{x:.10f}") beta_csv['t97'] = t97_padded DF['ydh_x'] = DF['ydh_x'].astype(str) # convert to string for later DF_complete = pd.merge(DF,beta_csv,on = 't97', how = 'inner')#,how = 'outer', indicator = True) return DF_complete def stat_analysis_hourly(n2_file,n3d_file, beta_file): df = get_dataframe(n2_file,n3d_file, beta_file) # implement filter: freq_min, freq_max = 100, 1200 snr_threshold = 200 freq, sn0, sn1 = df['f'], df['sn0'], df['sn1'] mask_freq1 = (freq >=100) & (freq <= 1200) #Frequency mask: #exclude frequencies within the intervals [-7,7],[97,107], ... bad_freqs = [[i*100 -7,i*100+7] for i in np.arange(13)] check_ = [(freq >interval[0]) & (freq < interval[1]) for interval in bad_freqs] reduced = np.logical_or.reduce(check_) mask_freq2 = ~reduced mask_freq = mask_freq1 & mask_freq2 #SNR mask mask_SNR = (sn0 >= snr_threshold) & (sn1 >= snr_threshold) #Beta mask: mask_beta_1 = (df['ant'] == 3) & (abs(df['beta_xw[\circ]']) >= 30) # if direction finding is off mask_beta_2 = (df['ant'] == 11 ) & (abs(df['beta_uw[\circ]']) >= 30) mask_beta_3 = (df['ant'] == 12 ) & (abs(df['beta_vw[\circ]']) >= 30) mask_beta = mask_beta_1 | mask_beta_2 | mask_beta_3 mask = mask_freq & mask_SNR & mask_beta df = df[mask] #number of datavalues num = len(df) # compute linear/circular/total polarisation for mask df['lp'] = np.sqrt(df['q']**2 + df['u']**2) df['circ_abs'] = np.abs(df['v']) df['tot'] = np.sqrt(df['q']**2 + df['u']**2 + df['v']**2) #vertical interferences mask: mask_vert = (df['lp'] <= 1.2) & (df['circ_abs'] <= 1.2) & (df['tot'] <= 1.2) df = df[mask_vert] lp = df['lp'] #linear pol circ = df['v'] # circular pol tot = df['tot'] # total pol s = df['s'] # intensity num_lp = len(df[lp > 0.2]) num_nolp = len(df[lp <= 0.2]) # ignore all cases where there are less than 10 datavalues if abs(num) < 10: return None cpf, lpf = num_nolp/num, num_lp/num mlin = np.mean(lp) mdlin = np.median(lp) dlin = np.std(lp) mcirc = np.mean(circ) mdcirc = np.median(circ) dcirc = np.std(circ) mtot = np.mean(tot) mdtot = np.median(tot) dtot = np.std(tot) ms = np.mean(s) mds = np.median(s) ds = np.std(s) Year, DOY, HH, MM = int(df['Year'].iat[0]), int(df['DOY'].iat[0]), int(df['HH'].iat[0]), int(df['MM'].iat[0]) new_row = {'Year': Year, 'DOY': DOY, 'HH': HH, 'num': num,'num_lp': num_lp, 'cpf': cpf,'lpf': lpf, 'mlin': mlin, 'mdlin': mdlin, 'dlin': dlin, 'mcirc': mcirc, 'mdcirc': mdcirc, 'dcirc': dcirc, 'mtot': mtot, 'mdtot': mdtot, 'dtot': dtot,'ms': ms, 'mds': mds, 'ds': ds} return new_row if __name__ == "__main__": analysis_csv = "C:/Users/doria/Desktop/Cassini/cassini-main/complete_analysis.csv" # output file check = 0 for year in range(2004,2018): beta_file = f"E:/CassiniData/ephemeris/{year}/n2times_beta_{year}.dat" for year_quarter in ["001_090","091_180","181_270","271_366"]: folder = "E:\HFRdata" folder_path = os.path.join(folder,f"{year}_{year_quarter}") folder_path_n2 = os.path.join(folder_path,"n2") n2_list,n3d_list = [],[] for n2_filename in os.listdir(folder_path_n2): ydh = n2_filename[1:] n3d_filename = f"N3d_dsq{ydh}" n2_file = os.path.join(folder_path_n2,n2_filename) n3d_file = os.path.join(folder_path,"n3d",n3d_filename) try: n3d_size = os.path.getsize(n3d_file) except FileNotFoundError: n3d_size = 0 print("Missing n3d-file") if n3d_size <= 2048: continue # ignore empty n3d files to avoid error n2_list.append(n2_file) n3d_list.append(n3d_file) with ProcessPoolExecutor() as executor: #Use multithreading to reduce computation time significantly rows = list(executor.map(stat_analysis_hourly,n2_list,n3d_list,repeat(beta_file))) rows = [x for x in rows if x is not None] df_quarter = pd.DataFrame(rows) df_quarter.to_csv(f'analysis_new_{year}_{year_quarter}.csv ') # final outputfile for yearly quarter