windsim/snr_check.py at master · solatale/windsim · GitHub

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

# Usage: python3 snr_check.py Cssos_magsim_SNR_tile_077_424.txt

import astropy.io.ascii as ascii
import numpy as np
from scipy import stats
import sys
import matplotlib.pyplot as plt
from scipy.stats import sigmaclip


def sigmaclip(ndarray, nsigma):
    mean = np.median(ndarray)
    std = np.std(ndarray)
    idx = np.where((ndarray>(mean-std*nsigma)) & (ndarray<(mean+std*nsigma)))
    return idx

def clip(ndarray, min, max):
    idx = np.where((ndarray>=min) & (ndarray<=max))
    return idx


snr = 10

data = ascii.read(sys.argv[1])
# data = ascii.read('Cssos_magsim_SNR_tile_077_424.txt')


slope, intercept, r_value, p_value, std_err = stats.linregress(data['MOD_r'],data['MagSim_r'])

print(slope,intercept)

plt.scatter(data['MOD_r'], data['MagSim_r'], s=1)
plt.plot(np.array([18,30]), slope*np.array([18,30])+intercept, color='red')
plt.xlim(18,28)
plt.ylim(18,28)
plt.show()


data['Resid_r'] = data['MagSim_r'] - (slope*data['MOD_r']+intercept)

data = data[sigmaclip(data['Resid_r'], 3)]
data = data[sigmaclip(data['Resid_r'], 3)]

dataclipsnr = data[clip(data['SNR_r'], snr-0.3, snr+0.3)]

plt.scatter(dataclipsnr['MOD_r'], dataclipsnr['MagSim_r'], s=1)
plt.plot(np.array([18,30]), slope*np.array([18,30])+intercept, color='red')
plt.xlim(min(dataclipsnr['MOD_r']),max(dataclipsnr['MOD_r']))
plt.ylim(min(dataclipsnr['MOD_r']),max(dataclipsnr['MOD_r']))
plt.show()

plt.scatter(dataclipsnr['MOD_r'], dataclipsnr['Resid_r'], s=1)
plt.show()

std_resid = np.std(dataclipsnr['Resid_r'])

print('SNR =', snr)
print('STD =', std_resid)
print(snr*std_resid)