ocssw V2020
mie_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.stats import lognorm
9 from scipy.integrate import quad
10 from scipy import integrate as int
11 from scipy import signal
12 from mpl_toolkits.mplot3d import Axes3D
13 from matplotlib import colors as mcolors
14 import matplotlib.pyplot as plt
15 from itertools import permutations
16 from random import sample
17 import cmath
18 import bhmie as bh
19 import mie_kern as mkrn
20 from pytmatrix import tmat
21 
22 
23 color_lvl = 8
24 rgb = list(permutations(range(0,256,color_lvl),3))
25 lnclrs = np.array(sample(rgb,120))
26 lnclrs = lnclrs/256.0
27 
28 start = 200
29 stop = 20000
30 m = 1.33+1.0E-10j
31 
32 mk = mkrn.mie_kern("_test_", 'oci', m, start, stop, 80)
33 mk.load('/accounts/ssander2/test/oci/kernel/KERN_137i10_80_071819.txt')
34 mk.kern[:,0] = 48*(mk.wl[0]/mk.wl[:])**4
35 mk.plot("kern")
36 
37 # Generate arbitrary size distribution
38 sd = np.zeros(mk.nsz)
39 for isz in range(mk.nsz):
40  sd[isz] = (1+np.cos(np.pi/100 + 2*np.pi*isz/50))/2
41 plt.clf()
42 #plt.plot(mk.sz,sd)
43 # plt.show()
44 # Generate a vector of reflectance values from the size distribution
45 refl = np.dot(mk.kern,sd)
46 plt.clf()
47 plt.plot(mk.wl,refl)
48 #plt.show()
49 """
50 rp0 = 1e-8
51 for j in range(10):
52  rp0 /= 10.0
53  rp1 = 1e-8
54  for q in range(10):
55  rp1 /= 10.0
56  rp2 = 1e-8
57  for k in range(2):
58  rp1 /= 1.0
59  rp2 /= 10.0
60  rp3 = 1e-4
61  for l in range(1):
62  rp1 /= 1.0
63  rp2 /= 1.0
64  rp3 /= 10.0
65  rikern = la.inv(np.dot(mk.kern.T,mk.kern) + rp0*mk.I + rp1*mk.H1 + rp2*mk.H2 + rp3*mk.H3)
66  psd = np.tensordot(np.dot(rikern,mk.kern.T),refl,(((1),(0))))
67  irefl = np.tensordot(mk.kern,psd,(((1),(0))))
68  err = la.norm(refl-irefl, ord=2, axis=0)/(la.norm(refl, ord=2, axis=0)+.001)
69 
70  plt.clf()
71  fig, (ax1,ax2) = plt.subplots(2,figsize=(40, 30))
72  fig.suptitle('Reflectance as a function of wavelength and particle size: ' + \
73  "{: .0e}".format(rp0) + " " + "{: .0e}".format(rp1) + " " + \
74  "{: .0e}".format(rp2) + " " + "{: .0e}".format(rp3))
75  ax1.set(xlabel='Wavelength (Nm)', ylabel='Reflectance',
76  title='Reflectance as a function of Wavelength')
77  ax1.grid(b=True, which='both', axis='both')
78  ax2.set(xlabel='Particle Size (Nm)', ylabel='Inversion',
79  title='Representation as a function of Particle Size')
80  ax2.grid(b=True, which='both', axis='both')
81 
82  ax1.plot(mk.wl,refl)
83  ax1.plot(mk.wl,irefl)
84  ax2.semilogx(mk.sz,psd)
85  ax2.semilogx(mk.sz,sd)
86  plt.show()
87 # plt.close(fig)
88 
89 """
90 
91 
92 # open simulated oci file and read reflectances into matrix refl
93 import tables
94 fin = tables.open_file('OCI2020084005000.L1B_PACE.nc', 'r')
95 oci = np.append(fin.get_node("/", "observation_data/Lt_blue")[8:59], \
96  fin.get_node("/", "observation_data/Lt_red")[1:,:,:],axis=0)
97 fin.close()
98 if not os.path.isdir("/accounts/ssander2/test/oci/figures"):
99  os.mkdir("/accounts/ssander2/test/oci/figures")
100 opath = "/accounts/ssander2/test/oci/figures"
101 
102 rp0 = 1e-7
103 for j in range(15):
104  rp0 /= 10.0
105  rp1 = 1e-0
106  for q in range(15):
107  rp1 /= 10.0
108  rp2 = 1e-0
109  for k in range(6):
110  rp1 /= 1.0
111  rp2 /= 10.0
112  rp3 = 1e-0
113  for l in range(10):
114  rp1 /= 1.0
115  rp2 /= 1.0
116  rp3 /= 10.0
117  rikern = la.inv(np.dot(mk.kern.T,mk.kern) + rp0*mk.I + rp1*mk.H1 + rp2*mk.H2 + rp3*mk.H3)
118  # psd = np.tensordot(np.dot(rikern,mk.kern.T),oci,(((1),(0)))) + \
119  # np.tensordot(rikern,rp1*np.multiply(np.ones_like(oci).T,guess).T,(((1),(0))))
120  psd = np.tensordot(np.dot(rikern,mk.kern.T),oci,(((1),(0))))
121  roci = np.tensordot(mk.kern,psd,(((1),(0))))
122 
123  print( "{: .0e}".format(rp0), "{: .0e}".format(rp1), "{: .0e}".format(rp2), "{: .0e}".format(rp3))
124  fname = opath + "/SZ_" + "{: .0e}".format(rp0) + "_" + "{: .0e}".format(rp1) + \
125  "_" + "{: .0e}".format(rp2) + "_" + "{: .0e}".format(rp3) + ".png"
126 
127  plt.clf()
128  fig, (ax1,ax2) = plt.subplots(2,figsize=(40, 30))
129  fig.suptitle('Reflectance as a function of wavelength and particle size: ' + \
130  "{: .0e}".format(rp0) + " " + "{: .0e}".format(rp1) + " " + \
131  "{: .0e}".format(rp2) + " " + "{: .0e}".format(rp3))
132  ax1.set(xlabel='Wavelength (Nm)', ylabel='Reflectance',
133  title='Reflectance as a function of Wavelength')
134  ax1.grid(b=True, which='both', axis='both')
135  ax2.set(xlabel='Particle Size (Nm)', ylabel='Inversion',
136  title='Representation as a function of Particle Size')
137  ax2.grid(b=True, which='both', axis='both')
138 
139  for i in range(0,1700,20):
140  ax1.plot(mk.wl,oci[:,i,600])
141  ax1.plot(mk.wl,roci[:,i,600], color=lnclrs[(13*i)%120])
142  ax2.semilogx(mk.sz,psd[:,i,600], color=lnclrs[(13*i)%120])
143  # plt.show(block=False)
144  plt.savefig(fname)
145  plt.close(fig)
146 
list(APPEND LIBS ${NETCDF_LIBRARIES}) list(APPEND LIBS $
Definition: CMakeLists.txt:9