ocssw V2020
oci_kernal.py
Go to the documentation of this file.
1 #!/usr/bin/env python
2 # encoding: utf-8
3 
4 import os
5 import sys
6 import numpy as np
7 from numpy import linalg as la
8 from scipy import signal
9 import matplotlib.pyplot as plt
10 from itertools import permutations
11 from random import sample
12 import cmath
13 import bhmie as bh
14 from pytmatrix import tmat
15 
16 def normalized(a, order=2, axis=-1 ):
17  l2 = np.atleast_1d(np.linalg.norm(a, order, axis))
18  l2[l2==0] = 1
19  return a / np.expand_dims(l2, axis)
20 
21 color_lvl = 8
22 rgb = list(permutations(range(0,256,color_lvl),3))
23 colors = np.array(sample(rgb,120))
24 colors = colors/256.0
25 
26 awl = [310, 315, 320, 325, 330, 335, 340, 345, 350, 355, 360, 365, 370, 375, 380, 385,
27  390, 395, 400, 405, 410, 415, 420, 425, 430, 435, 440, 445, 450, 455, 460, 465,
28  470, 475, 480, 485, 490, 495, 500, 505, 510, 515, 520, 525, 530, 535, 540, 545,
29  550, 555, 560, 565, 570, 575, 580, 585, 590, 595, 600, 605, 610, 615, 620, 625,
30  630, 635, 640, 645, 650, 655, 660, 665, 670, 675, 680, 685, 690, 695, 700, 705,
31  710, 715, 720, 725, 730, 735, 740, 745, 750, 755, 760, 765, 770, 775, 780, 785,
32  790, 795, 800, 805, 810, 815, 820, 825, 830, 835, 840, 845, 850, 855, 860, 865,
33  870, 875, 880, 885, 890, 895, 940, 1040, 1250, 1378, 1615, 2130, 2260 ]
34 
35 start = 200
36 stop = 20000
37 wl = np.copy(awl[8:-7])
38 nwl = len(wl)
39 x = np.arange(start,stop,1, )
40 nx = len(x)
41 lsz = np.linspace(np.log(start),np.log(stop),nwl)
42 sz = np.exp(lsz)
43 nsz = len(sz)
44 # regularization matrices
45 I = np.identity(nsz, dtype=np.float64)
46 # first difference matrix
47 H0 = -np.identity(nsz, dtype=np.float64)
48 H0[0,0] = 0
49 for i in range(1,nsz):
50  for j in range(nsz):
51  H0[i,j] = 1 if (j==i-1) else H0[i,j]
52 # sum of first differences squared
53 H1 = np.dot(H0.T,H0)
54 # sum of second differences squared
55 H2 = H1.copy()
56 H2[0,:] = 0
57 H2 = np.dot(H2.T,H2)
58 # sum of third differences squared
59 H3 = np.zeros([nsz,nwl])
60 for i in range(3,nsz):
61  for j in range(nwl):
62  if (j==i):
63  H3[i,j] = -1
64  elif (j==i-1):
65  H3[i,j] = 3
66  elif (j==i-2):
67  H3[i,j] = -3
68  elif (j==i-3):
69  H3[i,j] = 1
70 H3 = np.dot(H3.T,H3)
71 # The standard deviations of the particle size distributions for each size bin
72 #lsig = np.zeros(nsz)
73 #lsig[0] = np.log((sz[1]-sz[0])/2)
74 sig = np.zeros(nsz)
75 sig[0] = (sz[1]-sz[0])/2
76 for i in range(1,nsz-1):
77 # lsig[i] = (lsz[i+1]-lsz[i-1])/7.0
78  sig[i] = (sz[i+1]-sz[i-1])/4.5
79 sig[nsz-1] = (sz[nsz-1]-sz[nsz-2])/2
80 pdf = np.zeros((nx, nsz))
81 for i in range(1,nsz-1):
82 # pdf[:,i] = (np.exp(-(np.log(x[:])-lsz[i])**2/(2*lsig[i])**2))/(x[:]*lsig[i]*np.sqrt(2*np.pi))
83  pdf[:,i] = np.exp(-(x[:]-sz[i])**2/(2*sig[i]**2))/(sig[i]*np.sqrt(2*np.pi))
84 spdf = np.sum(pdf,axis=1)
85 plt.clf()
86 plt.semilogx(x, pdf)
87 plt.semilogx(x, spdf)
88 #plt.show()
89 # Index of refraction
90 m = 1.5+0.00001j
91 nang = 1
92 # T-Matrix code
93 #tm = tmat.tmat()
94 #tm.generate_tables()
95 #ktm = np.zeros((nwl, nsz), dtype=np.float64)
96 #for i in range(nwl):
97 # for j in range(int(nsz/4)):
98 # r, v, q, ktm[i,j], al, s, f = tm.mie(wl[i], m, sz[j], np.exp(lsig[j]), nang)
99 # if (j%10 == 0):
100 # print(str(j))
101 ng = 7
102 std = 1
103 of = np.linspace(-(ng-1)/2,(ng-1)/2, ng)
104 gsz = np.zeros(ng)
105 kbh = np.zeros((nwl, nsz), dtype=np.float64)
106 kbtz = np.zeros((nwl, nsz), dtype=np.float64)
107 kbhg = np.ones((ng, nwl, nsz), dtype=np.float64)
108 kbht = np.ones((nwl, nsz, ng), dtype=np.float64)
109 wnd = np.array(signal.gaussian(7, std=std))/(std * np.sqrt(2*np.pi))
110 # compute kernel based on Mie scattering
111 for i in range(nwl):
112  for j in range(nsz):
113  gsz[:] = sz[j] + of[:]*sig[j]
114  for k in range(ng):
115  s1,s2,qx,kbht[i,j,k],qb,gs = bh.bhmie(2*np.pi*gsz[k]/wl[i],m,nang)
116  kbh = np.tensordot(wnd,kbht,axes=([0],[2]))
117  kbtz = np.trapz(np.multiply(kbht,wnd),gsz,axis=2)
118 # Add Rayleigh scattering component to the first column
119 kbh[:,0] = 1*(wl[0]/wl[:])**4
120 plt.clf()
121 plt.plot(wl,kbh[:,0])
122 #plt.show()
123 # Normalize the kernal
124 kbh = normalized(kbh, axis=0)
125 # Generate a plausible vector of measurements by wavelength
126 R = 0.5
127 spect = np.ones(nwl)/5 + R*kbh[:,0]
128 plt.plot(wl,spect)
129 #plt.show()
130 # Generate a plausible size distribution
131 guess = np.ones(nsz)/10
132 guess[0] = 1
133 # Generate arbitrary size distribution
134 sd = np.zeros(nsz)
135 for i in range(nsz):
136  sd[i] = (1+np.cos(np.pi + 2*np.pi*i/50))/2
137 plt.clf()
138 plt.plot(sz,sd)
139 # plt.show()
140 # Generate a vector of reflectance values from the size distribution
141 aot = np.dot(kbh,sd)
142 plt.clf()
143 plt.plot(wl,aot)
144 # plt.show()
145 # Attempt to find an inversion matrix that correctly recovers the model
146 fig, (ax1, ax2) = plt.subplots(2,figsize=(20, 30))
147 fig.suptitle('Reflectance as a function of wavelength and particle size')
148 ax1.set(xlabel='Wavelength (Nm)', ylabel='Reflectance',
149  title='Reflectance as a function of Wavelength')
150 ax1.grid(b=True, which='both', axis='both')
151 ax2.set(xlabel='Particle Size (Nm)', ylabel='Reflectance',
152  title='Reflectance as a function of Particle Size')
153 ax2.grid(b=True, which='both', axis='both')
154 rp0 = 1e-6
155 rp1 = 0.001
156 rp2 = 0.001
157 rp3 = 0.001
158 for i in range(1):
159  rp0 /= 1
160  rp1n = rp1
161  rp2n = rp2
162  rp3n = rp3
163  plt.cla()
164  for j in range(10):
165  rp1n /= 10
166  rp2n = rp2
167  rp3n = rp3
168  ax2.set(xlabel='Particle Size (Nm)', ylabel='Reflectance',
169  title='Reflectance as a function of Particle Size')
170  ax2.grid(b=True, which='both', axis='both')
171  for k in range(10):
172  rp2n /= 10
173  rp3n = rp3
174  for l in range(10):
175  rp3n /= 10
176  rikern = la.inv(np.dot(kbh.T,kbh) + rp0*I + rp1n*H1 + rp2n*H2 + rp3n*H3)
177  # psd = np.tensordot(np.dot(rikern,kbh.T),oci,(((1),(0)))) + \
178  # np.tensordot(rikern,rp1*np.multiply(np.ones_like(oci).T,guess).T,(((1),(0))))
179  psd = np.tensordot(np.dot(rikern,kbh.T),spect,(((1),(0))))
180  rguess = np.tensordot(kbh,psd,(((1),(0))))
181 
182  print( rp0, rp1n, rp2n, rp3n )
183  ax1.plot(wl,spect, color=colors[(i+k+j) + 29])
184  ax1.plot(wl,rguess, color=colors[(i+k+j)])
185  ax2.semilogx(sz,psd, color=colors[(i+k+j)])
186 # plt.show(block=False)
187 
188 
189 # open simulated oci file and read reflectances into matrix refl
190 import tables
191 fin = tables.open_file('OCI2020084005000.L1B_PACE.nc', 'r')
192 oci = np.append(fin.get_node("/", "observation_data/Lt_blue")[8:59], \
193  fin.get_node("/", "observation_data/Lt_red")[1:,:,:],axis=0)
194 fin.close()
195 if not os.path.isdir("figures"):
196  os.mkdir("figures")
197 opath = "figures"
198 
199 rp0 = 0.00001
200 rp1 = 0.00001
201 rp2 = 0.00001
202 rp3 = 0.00001
203 
204 rikern = la.inv(np.dot(kbh.T,kbh) + rp0*I + rp1*H1 + rp2*H2 + rp3*H3)
205 #psd = np.tensordot(np.dot(rikern,kbh.T),oci,(((1),(0)))) + \
206 # np.tensordot(rikern,rpm*np.multiply(np.ones_like(oci).T,guess).T,(((1),(0))))
207 psd = np.tensordot(np.dot(rikern,kbh.T),oci,(((1),(0))))
208 roci = np.tensordot(kbh,psd,(((1),(0))))
209 
210 for i in range(0,1000,50):
211 
212  ax1.plot(wl,oci[:,500,i])
213  ax1.plot(wl,roci[:,500,i], color=colors[int(i/50)])
214  ax2.semilogx(sz,psd[:,500,i], color=colors[int(i/50)])
215 # plt.show(block=False)
216 
217 rp1 = 1e-5
218 for i in range(5):
219  rp0 = 1e-7
220  rp1 /= 10.0
221  rp2 = 1e-5
222  rp3 = 1e-5
223  for k in range(5):
224  rp1 /= 1.0
225  rp2 /= 10.0
226  rp3 = 1e-5
227  for l in range(5):
228  rp1 /= 1.0
229  rp2 /= 1.0
230  rp3 /= 10.0
231  rikern = la.inv(np.dot(kbh.T,kbh) + rp0*I + rp1*H1 + rp2*H2 + rp3*H3)
232  # psd = np.tensordot(np.dot(rikern,kbh.T),oci,(((1),(0)))) + \
233  # np.tensordot(rikern,rp1*np.multiply(np.ones_like(oci).T,guess).T,(((1),(0))))
234  psd = np.tensordot(np.dot(rikern,kbh.T),oci,(((1),(0))))
235  roci = np.tensordot(kbh,psd,(((1),(0))))
236 
237  print( "{: .0e}".format(rp0), "{: .0e}".format(rp1), "{: .0e}".format(rp2), "{: .0e}".format(rp3))
238  fname = opath + "/SZ_" + "{: .0e}".format(rp0) + "_" + "{: .0e}".format(rp1) + \
239  "_" + "{: .0e}".format(rp2) + "_" + "{: .0e}".format(rp3) + ".png"
240 
241  plt.clf()
242  fig, (ax1,ax2) = plt.subplots(2,figsize=(30, 30))
243  fig.suptitle('Reflectance as a function of wavelength and particle size: ' + \
244  "{: .0e}".format(rp0) + " " + "{: .0e}".format(rp1) + " " + \
245  "{: .0e}".format(rp2) + " " + "{: .0e}".format(rp3))
246  ax1.set(xlabel='Wavelength (Nm)', ylabel='Reflectance',
247  title='Reflectance as a function of Wavelength')
248  ax1.grid(b=True, which='both', axis='both')
249  ax2.set(xlabel='Particle Size (Nm)', ylabel='Reflectance',
250  title='Reflectance as a function of Particle Size')
251  ax2.grid(b=True, which='both', axis='both')
252 
253  for i in range(0,1000,20):
254  ax1.plot(wl,oci[:,500,i])
255  ax1.plot(wl,roci[:,500,i], color=colors[int(13*i)%120])
256  ax2.semilogx(sz,psd[:,500,i], color=colors[int(13*i)%120])
257 # plt.show(block=False)
258  plt.savefig(fname)
259 
def normalized(a, order=2, axis=-1)
Definition: oci_kernal.py:16
list(APPEND LIBS ${NETCDF_LIBRARIES}) list(APPEND LIBS $
Definition: CMakeLists.txt:9