-
Notifications
You must be signed in to change notification settings - Fork 0
/
plot_gems_snr_pix.pro
154 lines (122 loc) · 4.47 KB
/
plot_gems_snr_pix.pro
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
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
pro plot_gems_snr_pix, filename, xidx, yidx
sXidx = string(xidx, format='(i03)')
sYidx = string(yidx, format='(i03)')
;------------------------------------------------------------------------------
; Configuration
;------------------------------------------------------------------------------
pos = [0.12, 0.1, 0.9, 0.85]
leg_pos = [0.85, 0.8]
winliminitpos = strpos(filename, 'winliminit')
fitrange = strmid(filename, winliminitpos+10, 3) + '340'
fit_range = strmid(fitrange, 0, 3) + '-' + strmid(fitrange, 3, 3)
L1Cmaker = 'EOSRL'
path ='/data1/gems/o3p/ds/GEMS_O3P_Yonsei/'
outputpath = path + 'out/'
projectpath = outputpath ;+ 'softcal/' + fitrange + '/'
filelist = file_search(projectpath + '*.nc4')
runtime = strmid(filelist, 13, 10, /reverse)
runtimeidx = sort(runtime)
recentrunfile = filelist[runtimeidx[-1]]
;fn = file_basename(recentrunfile)
;fn = file_basename(filefullpath)
datepos = stregex(filename,'[0-9]{8}_[0-9]{4}', length=len)
date = strmid(filename, datepos, len)
; xidx=89, yidx=43
scp_dest = '[email protected]:/home/soodal/works/plot/'
; read GEMS_L2_O3P
data = ds_read_gems_l2_o3p(filename, $
varlist=['EffectiveCloudFractionUV', $
'ResidualsOfFit', 'Latitude', 'Longitude', $
'SimulatedRadiances', 'FitWeights', 'SignalToNoiseRatio', 'WavelengthsWholeRange'])
basename = file_basename(filename)
; get variable from GEMS_L2_O3P
ecf = data.EffectiveCloudFractionUV
ecfnanidx = where(ecf lt -990, /null)
ecf[ecfnanidx] = !values.f_nan
lat = data.latitude
latnanidx = where(lat lt -990, /null)
lat[latnanidx] = !values.f_nan
lon = data.longitude
lonnanidx = where(lon lt -990, /null)
lon[lonnanidx] = !values.f_nan
fitres = data.ResidualsOfFit
nanidx = where(fitres lt -990)
fitres[nanidx] = !values.f_nan
zeroidx = where(fitres eq 0)
fitres[zeroidx] = !values.f_nan
fitresdims = size(fitres, /dim)
simrad = data.SimulatedRadiances
nanidx = where(simrad lt -990)
simrad[nanidx] = !values.f_nan
zeroidx = where(simrad eq 0)
simrad[zeroidx] = !values.f_nan
simraddims = size(simrad, /dim)
fitweights = data.FitWeights
snr = data.SignalToNoiseRatio
wav = data.WavelengthsWholeRange
;wavpath = strmid(fn, 13, 10, /reverse) + '/'
;wavpath = ''
;wavfn = projectpath + wavpath + 'waves_rad_310_340.txt'
;readcol, wavfn, wav
;==============================================================================
; Plot
;------------------------------------------------------------------------------
; for latitude < 20
;------------------------------------------------------------------------------
latclr_idx = where(lat lt 20 and ecf lt 0.2)
;fitres = reform(fitres, fitresdims[0], 174L*512)
;simrad = reform(simrad, simraddims[0], 174L*512)
residuals = fltarr(fitresdims[0])
residuals[*] = !values.f_nan
stdres = fltarr(fitresdims[0])
stdres[*] = !values.f_nan
npix = n_elements(latclr_idx)
; pick 10 pixel from 1/20, 3/20, 5/20, ..., 19/20 from the whole idx
;selected_idx = (indgen(10) * 2 + 1) * npix/20
;for in = 0, n_elements(selected_idx)-1 do begin
;for iwav = 0, fitresdims[0]-1 do begin
;_fitres = reform(fitres[iwav, latclr_idx[in]])
;_simrad = reform(simrad[iwav, latclr_idx[in]])
;residuals[iwav] = mean(_fitres, /nan)/mean(_simrad, /nan)*100
;stdres[iwav] = stddev(_fitres, /nan)/mean(_simrad, /nan)*100
;endfor
residuals = fitres[*, xidx, yidx]$
/simrad[*, xidx, yidx] * 100
wav0idx = where(wav eq 0, /null)
wav[wav0idx] = !values.f_nan
residuals[wav0idx] = !values.f_nan
mm = minmax(residuals)
yrange = [-4, 4]
;IF mm[1] GE 0.3 THEN BEGIN
;yrange[1] = mm[1]
;ENDIf
;IF mm[0] LE -0.3 THEN BEGIN
;yrange[0] = mm[0]
;ENDIF
p1 = plot(wav[0:203], 1/snr[0:203, yidx]*100, /buffer, $
xtitle='Wavelength[nm]', $
ytitle='Measurement Error [1/SNR] [%]', $
title='GEMS L2 O3P ' + fit_range + 'nm Measurement Error from SNR', $
symbol='s', $
yrange=[0, 1], $
xrange=[300, 340], $
;name='\c Residuals', $
position=pos)
p_ = plot(indgen(500), fltarr(500), 'gray', /buffer, /over)
;p2 = plot(wav[0:fitresdims[0]-1], stdres, 'r', /buffer, $
;/over, name='Stddev')
;t1 = text(0.25, 0.75, 'Effective Cloud Fraction < 0.2', /norm)
;t2 = text(0.25, 0.70, 'Latitude < 20', /norm)
;leg = legend(target=[p0, p1, p2, p3, pap], pos=leg_pos)
pngfile = './plot/' +basename + '_snr_x' +sXidx +'_y' +sYidx +'.png'
p1.errorbar_color='indian_red'
p1.errorbar_capsize=0.1
p1.title.font_size=16
p1.save, pngfile
p1.close
; send image to pc
spawn, 'scp -P18742 -p ' + pngfile + ' ' + scp_dest
;endfor
print, 'Done.'
stop
end