-
Notifications
You must be signed in to change notification settings - Fork 3
/
Image_PowerCorrection.py
76 lines (65 loc) · 3.33 KB
/
Image_PowerCorrection.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
import xarray as xr
import numpy as np
def Image_PowerCorrection(radar_image,flight_elev,depth_axis,attenuation_val,attenuation_type):
"""
% (C) Nick Holschuh - Amherst College - 2024 ([email protected])
% This function, applied to a depth_shifted image, applies a depth_power
% correction that allows you to interpret relative reflectivity within the ice column
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The inputs are as follows:
% radar_image --
% flight_elev --
% depth_axis --
% attenuation_val --
% attenuation_type -- 0: 1 wawy attenuation rate, 1: two way attenuation rate, 2: nx2 array with depths and attenuation rate values
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The outputs are as follows:
% power_corrected_radarim -- radar image (in dB) corrected for spreading and attenuation
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% e.g.:
radar_image = depth_data['new_data']
flight_elev = radar_data['Elevation']
depth_axis = depth_data['depth_axis']
attenuation_val = 8.0
attenuation_type = 0
"""
################ This is the import statement required to reference scripts within the package
import os,sys,glob
ndh_tools_path_opts = [
'/mnt/data01/Code/',
'/home/common/HolschuhLab/Code/'
]
for i in ndh_tools_path_opts:
if os.path.isdir(i): sys.path.append(i); correction_root_dir=i;
################################################################################################
############## Load in the pre-calculated spreading corrections
##### Built from Matlab 'Generate_SpreadingMatrix.m'
##### Converted to NC from 'Develop_spreadingcorrection.ipynb'
spreading_correction_vals = xr.open_dataset(correction_root_dir+'NDH_Tools/SpreadingCorrection.nc')
############## Construct a synthetic twtt for the depth_image
cair = 299792458
cice = 1.68e8
dz = np.mean(np.diff(depth_axis))
imdims = radar_image.shape
depth_image = np.cumsum(np.ones(imdims),axis=0)*dz
twtt_image = depth_image*2/cice+flight_elev*2/cair
elev_image = np.zeros(imdims)+flight_elev
x_search = xr.DataArray(elev_image.ravel(),dims=['vector_index'])
y_search = xr.DataArray(twtt_image.ravel(),dims=['vector_index'])
spreading_losses = spreading_correction_vals.interp(flight_elev=x_search,twtt=y_search)
spreading_losses = spreading_losses['spreadingloss_raytracing'].values.reshape(imdims)
######### One way attenuation Rate
if attenuation_type == 0:
attenuation_losses = depth_image*2*np.abs(attenuation_val)/1000
######### Two way attenuation Rate
elif attenuation_type == 1:
attenuation_losses = depth_image*np.abs(attenuation_val)/1000
######### Attenuation profile (assumes two way)
elif attenuation_type == 2:
attenuation_val = np.interp1(attenuation_val[:,0],np.abs(attenuation_val[:,1]),depth_axis)
attenuation_val = np.zeros(imdims)+attenuation_val
attenuation_losses = np.cumsum(attenuation_val,axis=0)*1000*2
power_corrected_radarim = 10*np.log10(radar_image)-spreading_losses+attenuation_losses
return power_corrected_radarim