{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.14","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":70367,"databundleVersionId":9188054,"sourceType":"competition"}],"dockerImageVersionId":30761,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"from pathlib import Path\n\nimport numpy as np\nimport holoviews as hv\nfrom scipy.signal import savgol_filter, find_peaks\n\nhv.extension('bokeh')","metadata":{"execution":{"iopub.status.busy":"2024-08-22T13:31:02.195013Z","iopub.execute_input":"2024-08-22T13:31:02.195605Z","iopub.status.idle":"2024-08-22T13:31:02.609748Z","shell.execute_reply.started":"2024-08-22T13:31:02.195556Z","shell.execute_reply":"2024-08-22T13:31:02.60789Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"airs_dir=Path('/kaggle/working/')\nairs_paths = [i for i in airs_dir.rglob('airs_signal_*.npy')]","metadata":{"execution":{"iopub.status.busy":"2024-08-22T13:31:02.897321Z","iopub.execute_input":"2024-08-22T13:31:02.897776Z","iopub.status.idle":"2024-08-22T13:31:02.905564Z","shell.execute_reply.started":"2024-08-22T13:31:02.897735Z","shell.execute_reply":"2024-08-22T13:31:02.904Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"I used the second derivative of a signal to find possible locations corresponding to the inside of the transition zone. By signal, I mean the sum of the **AIRS-CH0** calibrated data of shape (N, 32, 282) over spatial and wavelegth dimensions, N is the number of time points.\nThe two maximum values of the second derivative mark the beginning and end of the transition zone inside it.\nIts two minimum mark the beginnign and end of the transition zone just outside of it.\n\nThree Parameters to tune:\nThese peaks are found and then are slightly shiftted for a better estimation using **inner_shift** and **outter_shift** parameters. \nThe other parameter to tune is the **height** parameter in **find_peaks**.","metadata":{}},{"cell_type":"code","source":"#Overlay plot of signal with its two derivatives\nairs_path=airs_paths[0]\ndata=np.load(airs_path)\nplanet_id=airs_path.stem.split('_')[-1]\nmask=np.load(airs_dir/f'airs_mask_{planet_id}.npy')\n# The mask belogns to hot and dead pixels. Invert it\n# and multiply by the signal to set the value of those pixels to zero.\ntime_series_data=(~mask*data).sum(axis=(1,2))\ntime_series_data/=time_series_data.max()\ngradient=np.gradient(time_series_data)\ngradient/=gradient.max()\ngradient_2=np.gradient(gradient)\ngradient_2/=gradient_2.max()\n(hv.Curve(time_series_data, vdims='signal', label='signal').opts(width=800, ylabel='normalized intensity', xlabel='time index',autorange='y')*\\\nhv.Curve(gradient, label='1st derivative').opts(color='black', ylabel='normalized derivative amplitude')*hv.Curve(gradient_2, label='2nd derivative')).opts(\n    multi_y=True, legend_position='bottom_left')","metadata":{"execution":{"iopub.status.busy":"2024-08-22T13:31:03.537133Z","iopub.execute_input":"2024-08-22T13:31:03.537578Z","iopub.status.idle":"2024-08-22T13:31:03.825193Z","shell.execute_reply.started":"2024-08-22T13:31:03.537539Z","shell.execute_reply":"2024-08-22T13:31:03.823838Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def get_peaks(airs_mask:np.array, airs_sig:np.array, inner_shift:int, outter_shift:int)->dict[tuple[int,int], tuple[int, int]]:\n    \"\"\"\n    Finds four indices in airs_sig. Two of those mark the transition inside the transiztion zone.\n    The other two mark the transition zone just outside of it.\n    Parameters\n    airs_mask: The mask of hot and dead pixels to be applied on the calibrated signal.\n    airs_sig: Calibrated signal. I did not apply the mask of hot and dead pixels during calibration so here I need to do it.\n    inner_shift: shift the two indecies in airs_sig identified to exist inside the the transition zone inward, that is to the right \n    (for the start of the transition) and to the left (for the end of the transition) to \n    make the identified transition zone tigher if it needs it. inner_shift's unit is integer index.\n    outter_shift: Similar to inner_shift, outter_shift pushes the locations found just outside of the transitions zone outward\n    to compensate for inaccurate estimation of transition zone.\n    \n    \"\"\"\n    # appllying a mask on the calibrated data. I did avoid using masked arrays during calibration.\n    masked_signal=~airs_mask*airs_sig\n    # Smoothing the signal for a better derivative.\n    smoothed_signal=savgol_filter((masked_signal).sum(axis=(1,2)), 10, 3)\n    # first derivative\n    gradient=np.gradient(smoothed_signal)\n    gradient/=gradient.max()\n    # second derivative\n    gradient_2 = np.gradient(gradient)\n    gradient_2/=gradient_2.max()\n    # find the inner peaks inside the transition zone\n    # If you get more than two peaks, change the height value\n    in_transition_peaks,_=find_peaks(gradient_2, height=0.6)\n    # Multiply by -1, to get the minmimum peaks by using the height option\n    # which seems only works for positive values.\n    # find outter peaks just outside of the transition zone\n    out_of_transition_peaks, _ =find_peaks(-1*gradient_2, height=0.6)\n    return {'inner_peaks': (in_transition_peaks[0]+inner_shift, in_transition_peaks[1]-inner_shift), \n            'outter_peaks': (out_of_transition_peaks[0]-outter_shift, out_of_transition_peaks[1]+outter_shift)}\n\ndef get_plots(planet_id:str, smoothed_signal:np.array, masked_signal:np.array, \n              in_transition_peaks:tuple[int, int], \n              out_of_transition_peaks:tuple[int,int])->hv.Overlay:\n    return (hv.Curve(smoothed_signal, label='smoothed signal')*\\\n    hv.Curve((masked_signal).sum(axis=(1,2)), label='signal')*\\\n    hv.VLines([in_transition_peaks[0],in_transition_peaks[1]]).opts(color='green', line_width=2, width=600)*\\\n    hv.VLines([out_of_transition_peaks[0],out_of_transition_peaks[1]]).opts(color='black', line_width=2, width=600)).opts(\n        legend_position='bottom_left', ylabel='intensity', xlabel='time index', title=f'planet_id:{planet_id}',axiswise=True)","metadata":{"execution":{"iopub.status.busy":"2024-08-22T13:31:03.8277Z","iopub.execute_input":"2024-08-22T13:31:03.828127Z","iopub.status.idle":"2024-08-22T13:31:03.844864Z","shell.execute_reply.started":"2024-08-22T13:31:03.828083Z","shell.execute_reply":"2024-08-22T13:31:03.843245Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for index, airs_path in enumerate(airs_paths):\n    # loads signal\n    airs_sig=np.load(airs_path)\n    # gets planet's id\n    planet_id=airs_path.stem.split('_')[-1]\n    # loads the mask for the signal\n    airs_mask=np.load(airs_dir/f'airs_mask_{planet_id}.npy')\n    # get ths four peaks\n    peaks = get_peaks(airs_mask=airs_mask, airs_sig=airs_sig, inner_shift=3, outter_shift=3)\n    # get the masked signal for plotting\n    masked_signal=~airs_mask*airs_sig\n    # get the smoothed signal for plotting\n    smoothed_signal=savgol_filter((masked_signal).sum(axis=(1,2)), 10, 3)\n    # The if statement is for making a holoviews overlay using the first plot\n    if index==0:\n        plots=get_plots(planet_id=planet_id, smoothed_signal=smoothed_signal, masked_signal=masked_signal, \n                  in_transition_peaks=peaks['inner_peaks'], \n                  out_of_transition_peaks=peaks['outter_peaks'])\n    else:\n        plots+=get_plots(planet_id=planet_id, smoothed_signal=smoothed_signal, masked_signal=masked_signal, \n                  in_transition_peaks=peaks['inner_peaks'], \n                  out_of_transition_peaks=peaks['outter_peaks'])","metadata":{"execution":{"iopub.status.busy":"2024-08-22T13:31:03.963066Z","iopub.execute_input":"2024-08-22T13:31:03.964106Z","iopub.status.idle":"2024-08-22T13:31:04.460875Z","shell.execute_reply.started":"2024-08-22T13:31:03.964054Z","shell.execute_reply":"2024-08-22T13:31:04.459576Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plots.opts(shared_axes=False).cols(2)","metadata":{"execution":{"iopub.status.busy":"2024-08-22T13:31:04.463253Z","iopub.execute_input":"2024-08-22T13:31:04.463682Z","iopub.status.idle":"2024-08-22T13:31:07.636639Z","shell.execute_reply.started":"2024-08-22T13:31:04.46364Z","shell.execute_reply":"2024-08-22T13:31:07.63496Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}