From 260e8819b736bcd48b15ddae593520c04a927e8a Mon Sep 17 00:00:00 2001 From: Fasil Kebede <46403603+FasilGibdaw@users.noreply.github.com> Date: Fri, 9 Feb 2024 17:57:53 +0200 Subject: [PATCH] Add files via upload --- scripts/ulf_waves_filter.py | 128 ++++++++++++++++++++++++++++++++++++ 1 file changed, 128 insertions(+) create mode 100644 scripts/ulf_waves_filter.py diff --git a/scripts/ulf_waves_filter.py b/scripts/ulf_waves_filter.py new file mode 100644 index 00000000..3c504158 --- /dev/null +++ b/scripts/ulf_waves_filter.py @@ -0,0 +1,128 @@ +import os # nopep8 +os.environ['PTNOLATEX'] = '1' # nopep8 +import pytools as pt # nopep8 +import glob +from scipy.signal import butter, filtfilt +import numpy as np + + +# fileNames={i : "/wrk-vakka/group/spacephysics/vlasiator/3D/EGI/bulk/dense_cold_hall1e5_afterRestart374/bulk1.{:07d}.vlsv".format(i) for i in range(662,1506)} +# fileNames={i : "/wrk-vakka/group/spacephysics/vlasiator/3D/FHA/bulk1/bulk1.{:07d}.vlsv".format(i) for i in range(501,1612)} +# fileNames={i :'/wrk-vakka/group/spacephysics/vlasiator/3D/EGL/bulk/bulk1.egl.{:07d}.vlsv'.format(i) for i in range(621,1760)} +# print(fileNames) + + +def ulf_filter( + filesDirectory="/wrk-vakka/group/spacephysics/vlasiator/3D/EGL/bulk/", + var_to_filter="vg_b_vol", + window_pad=50, + timestate=1000, + filter_order=5, + target_wave="Pc2", + fs=1, + outputDirectory="temp/", + run="EGL", +): + """Function to compute the moving average of electromagentic fields (or any quantities for that matter), and apply band pass filter + + Args: + filesDirectory (str): vlsv files directory + var_to_filter (str, optional): variables (vg_b_vol, vg_e_vol, fg_e (computationally costy), fg_b (computationally costy)). Defaults to "vg_b_vol". + window_pad (int, optional): moving average window. Defaults to 50 for 100 seconds moving average. + timestate (int, optional): (start_of_simulation_time + window_pad < timestate < end_of_simulation_time - window_pad). Defaults to 1000. + filter_order (int, optional): band pass filter order. Defaults to 5. + target_wave (str, optional): ULF wave target (Pc2, Pc3, Pc4), currently resolved waves in Vlasiator. Defaults to "Pc2". + fs (int, optional): samapling frequency. Defaults to 1. + outputDirectory (str, optional): directory where the sidecars will be saved. Defaults to "/proj/kebedefa/temp/EGL/". + output: + vlsv file saved in the output directory + """ + # var_to_filter = "vg_b_vol" + fileNames = glob.glob(filesDirectory + "*.vlsv") + fileNames.sort() + # start of the run (assumes continous simulation data) + run_start = int(fileNames[0].split('.')[-2]) + f0 = pt.vlsvfile.VlsvReader(fileNames[0]) + cids0 = f0.read_variable("CellID") + ncells = len(cids0) + nelems = len(f0.read_variable(var_to_filter, [1])) + windowlength = 2 * window_pad + 1 + B_all = np.zeros((windowlength, ncells, nelems)) + print(B_all.shape) + # timestate = 1000 # this can take from 713 to 1450 ( for example for EGI) + average_range = (timestate - window_pad, timestate + window_pad) + # collecting data for moving_averaging based on the window_pad/average_range + for i, t in enumerate(range(average_range[0], average_range[1] + 1)): + print("loading in ", (i, t), fileNames[t-run_start]) + reader = pt.vlsvfile.VlsvReader(fileNames[t-run_start]) + cids = reader.read_variable("CellID") + sorti = np.argsort(cids) + var = reader.read_variable(var_to_filter)[sorti, :] + B_all[i, :, :] = var + # filter_order =5 + if target_wave == "Pc2": + lowcut = 0.1 + highcut = 0.45 + elif target_wave == "Pc3": + lowcut = 0.02 + highcut = 0.1 + elif target_wave == "Pc4": + lowcut = 0.006 + highcut = 0.02 + else: + print("no defined band filter") + b, a = butter(filter_order, [lowcut, highcut], fs=fs, btype="band") + # b, a = butter(filter_order, 0.05, btype='high', fs=fs) + filtered_dataPc2 = filtfilt(b, a, B_all, axis=0) + averaged = np.mean(B_all, axis=0) + reader = pt.vlsvfile.VlsvReader(fileNames[timestate]) + + if not os.path.exists(outputDirectory): + # If it doesn't exist, create it + os.makedirs(outputDirectory) + + writer = pt.vlsvfile.VlsvWriter( + reader, + "{:s}{:s}_test_movingaverage_{:07d}.vlsv".format( + outputDirectory, run, timestate + ), + copy_meshes="SpatialGrid", + ) + # Uncomment this if you want to copy all the vlsv varibale, but not needed for now + # writer.copy_variables_list(reader, reader.get_all_variables()) + sorti = np.argsort(reader.read_variable("CellID")) + rev_sorti = np.argsort(sorti) + data_arr = averaged[rev_sorti, :] + delta_var = reader.read_variable(var_to_filter) + delta_var = delta_var - data_arr + # moving average of var + varinfo = pt.calculations.VariableInfo( + data_arr, + "{:s}_move_ave_{:d}".format(var_to_filter, windowlength), + units="T", + latex="$B_{\mathrm{ave}}$", + latexunits="$\mathrm{T}$", + ) + writer.write_variable_info(varinfo, "SpatialGrid", unitConversion=1) + varinfo = pt.calculations.VariableInfo( + delta_var, + "{:s}_move_ave_{:d}_delta".format(var_to_filter, windowlength), + units="T", + latex="$\delta{}B_{\mathrm{ave}}$", + latexunits="$\mathrm{T}$", + ) + writer.write_variable_info(varinfo, "SpatialGrid", unitConversion=1) + # window_pad =50 see above this refers to the target file (centered file) + varinfo = pt.calculations.VariableInfo( + filtered_dataPc2[window_pad, rev_sorti, :], + "vg_b_vol_{:s}_{:d}_xyz".format(target_wave, windowlength), + units="T^2Hz^-1", + latex="$\delta{}B_{\mathrm{ave}}$", + latexunits="$\mathrm{T}^2\,\mathrm{Hz}^{-1}\mathrm{orsomething}$", + ) + writer.write_variable_info(varinfo, "SpatialGrid", unitConversion=1) + del reader + + +if __name__ == "__main__": + ulf_filter()