OB.DAAC Logo
NASA Logo
Ocean Color Science Software

ocssw V2022
dbocean.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 math
7 import time
8 import argparse
9 import numpy as np
10 import lmfit as lm
11 import xarray as xr
12 import matplotlib.pyplot as plt
13 import scipy.interpolate as trp
14 
15 W470 = 0
16 W550 = 1
17 W659 = 2
18 W860 = 3
19 W124 = 4
20 W164 = 5
21 W213 = 6
22 NWL = 7
23 NSZA = 22
24 NVZA = 20
25 NRAA = 21
26 NAOT = 14
27 NFMF = 14
28 NF1 = 5
29 NF2 = 9
30 NWS = 6
31 NCHL = 4
32 D2R = np.pi/180.0
33 
34 class dbocean(object):
35 
36  def __init__(self, lut_filepath, mode):
37  print ("Reading Deepblue LUT: " + lut_filepath)
38  self.mode = mode
39  self.chlc = -2.0
40  self.wndc = 1.0
41  self.type = 0
42  self.lut = np.zeros((NWL,NCHL,NWS,NFMF,NAOT,NRAA,NVZA,NSZA))
43  fmf_pts = np.zeros((NFMF))
44  try:
45  self.lut[W470][:,:,0:NF1,:,:,:,:]=xr.load_dataset(lut_filepath,group='/OCEAN_AEROSOL_DUST')['IoverF_m03']
46  self.lut[W470][:,:,NF1-1:NF2-1,:,:,:,:]=xr.load_dataset(lut_filepath,group='/OCEAN_AEROSOL_MIXED')['IoverF_m03']
47  self.lut[W470][:,:,NF2-2:NFMF,:,:,:,:]=xr.load_dataset(lut_filepath,group='/OCEAN_AEROSOL_FINE')['IoverF_m03']
48  self.lut[W550][:,:,0:NF1,:,:,:,:]=xr.load_dataset(lut_filepath,group='/OCEAN_AEROSOL_DUST')['IoverF_m04']
49  self.lut[W550][:,:,NF1-1:NF2-1,:,:,:,:]=xr.load_dataset(lut_filepath,group='/OCEAN_AEROSOL_MIXED')['IoverF_m04']
50  self.lut[W550][:,:,NF2-2:NFMF,:,:,:,:]=xr.load_dataset(lut_filepath,group='/OCEAN_AEROSOL_FINE')['IoverF_m04']
51  self.lut[W659][:,:,0:NF1,:,:,:,:]=xr.load_dataset(lut_filepath,group='/OCEAN_AEROSOL_DUST')['IoverF_m05']
52  self.lut[W659][:,:,NF1-1:NF2-1,:,:,:,:]=xr.load_dataset(lut_filepath,group='/OCEAN_AEROSOL_MIXED')['IoverF_m05']
53  self.lut[W659][:,:,NF2-2:NFMF,:,:,:,:]=xr.load_dataset(lut_filepath,group='/OCEAN_AEROSOL_FINE')['IoverF_m05']
54  self.lut[W860][:,:,0:NF1,:,:,:,:]=xr.load_dataset(lut_filepath,group='/OCEAN_AEROSOL_DUST')['IoverF_m07']
55  self.lut[W860][:,:,NF1-1:NF2-1,:,:,:,:]=xr.load_dataset(lut_filepath,group='/OCEAN_AEROSOL_MIXED')['IoverF_m07']
56  self.lut[W860][:,:,NF2-2:NFMF,:,:,:,:]=xr.load_dataset(lut_filepath,group='/OCEAN_AEROSOL_FINE')['IoverF_m07']
57  self.lut[W124][:,:,0:NF1,:,:,:,:]=xr.load_dataset(lut_filepath,group='/OCEAN_AEROSOL_DUST')['IoverF_m08']
58  self.lut[W124][:,:,NF1-1:NF2-1,:,:,:,:]=xr.load_dataset(lut_filepath,group='/OCEAN_AEROSOL_MIXED')['IoverF_m08']
59  self.lut[W124][:,:,NF2-2:NFMF,:,:,:,:]=xr.load_dataset(lut_filepath,group='/OCEAN_AEROSOL_FINE')['IoverF_m08']
60  self.lut[W164][:,:,0:NF1,:,:,:,:]=xr.load_dataset(lut_filepath,group='/OCEAN_AEROSOL_DUST')['IoverF_m10']
61  self.lut[W164][:,:,NF1-1:NF2-1,:,:,:,:]=xr.load_dataset(lut_filepath,group='/OCEAN_AEROSOL_MIXED')['IoverF_m10']
62  self.lut[W164][:,:,NF2-2:NFMF,:,:,:,:]=xr.load_dataset(lut_filepath,group='/OCEAN_AEROSOL_FINE')['IoverF_m10']
63  self.lut[W213][:,:,0:NF1,:,:,:,:]=xr.load_dataset(lut_filepath,group='/OCEAN_AEROSOL_DUST')['IoverF_m11']
64  self.lut[W213][:,:,NF1-1:NF2-1,:,:,:,:]=xr.load_dataset(lut_filepath,group='/OCEAN_AEROSOL_MIXED')['IoverF_m11']
65  self.lut[W213][:,:,NF2-2:NFMF,:,:,:,:]=xr.load_dataset(lut_filepath,group='/OCEAN_AEROSOL_FINE')['IoverF_m11']
66  self.chl_lut = xr.load_dataset(lut_filepath,group='/LOG_CHL')['LOG_CHL']
67  raa_pts = xr.load_dataset(lut_filepath,group='/OCEAN_AEROSOL_MIXED')['Relative_Azimuth_Angle']
68  vza_pts = xr.load_dataset(lut_filepath,group='/OCEAN_AEROSOL_MIXED')['View_Zenith_Angle']
69  sza_pts = xr.load_dataset(lut_filepath,group='/OCEAN_AEROSOL_MIXED')['Solar_Zenith_Angle']
70  wnd_pts = xr.load_dataset(lut_filepath,group='/OCEAN_AEROSOL_MIXED')['Wind_Speed']
71  chl_pts = xr.load_dataset(lut_filepath,group='/OCEAN_AEROSOL_MIXED')['Chl_Conc']
72  aot_pts = xr.load_dataset(lut_filepath,group='/OCEAN_AEROSOL_MIXED')['Aerosol_Optical_Depth_550']
73  fmf_pts[0:NF1] = xr.load_dataset(lut_filepath,group='/OCEAN_AEROSOL_DUST')['Fine_Mode_Fraction_550']
74  fmf_pts[NF1-1:NF2-1] = xr.load_dataset(lut_filepath,group='/OCEAN_AEROSOL_MIXED')['Fine_Mode_Fraction_550']
75  fmf_pts[NF2-2:NFMF] = xr.load_dataset(lut_filepath,group='/OCEAN_AEROSOL_FINE')['Fine_Mode_Fraction_550']
76  self.wl_pts = xr.load_dataset(lut_filepath,group='/OCEAN_AEROSOL_MIXED')['Band_Central_Wavelength']
77  self.coef = np.array([0.06,0.06,0.04,0.04,0.07,0.06,0.1])
78  self.pars = lm.Parameters()
79  self.pars.add(name='fmf', value=0.5, min=0, max=1)
80  self.pars.add(name='aot', value=1.0, min=0, max=100)
81  if(self.mode==0 or self.mode==1):
82  self.rpts = (chl_pts, wnd_pts, raa_pts, vza_pts, sza_pts)
83  self.mpts = (fmf_pts, aot_pts)
84  self.lut = np.transpose(self.lut, (1,2,5,6,7,3,4,0))
85  elif(self.mode==2):
86  self.rpts = (wnd_pts, raa_pts, vza_pts, sza_pts)
87  self.mpts = (fmf_pts, aot_pts, chl_pts)
88  self.lut = np.transpose(self.lut, (2,5,6,7,3,4,1,0))
89  self.pars.add(name='chl', value=-2.0, min=-10.0, max=2.0)
90  elif(self.mode==3):
91  self.rpts = (chl_pts, raa_pts, vza_pts, sza_pts)
92  self.mpts = (fmf_pts, aot_pts, wnd_pts)
93  self.lut = np.transpose(self.lut, (1,5,6,7,3,4,2,0))
94  self.pars.add(name='wnd', value=1.0, min=0, max=100)
95  elif(self.mode==4):
96  self.rpts = (raa_pts, vza_pts, sza_pts)
97  self.mpts = (fmf_pts, aot_pts, chl_pts, wnd_pts)
98  self.lut = np.transpose(self.lut, (5,6,7,3,4,1,2,0))
99  self.pars.add(name='chl', value=-2.0, min=-10.0, max=2.0)
100  self.pars.add(name='wnd', value=1.0, min=0, max=100)
101  except Exception as inst:
102  print(type(inst))
103  print(inst)
104  print ("Unable to read LUT file ... exiting")
105  sys.exit()
106 
107  def minfun(self, pars, data, scale, rlut):
108  if(self.mode==0 or self.mode==1):
109  rxi = np.stack((pars['fmf'], pars['aot']))
110  elif(self.mode==2):
111  rxi = np.stack((pars['fmf'], pars['aot'], pars['chl']))
112  elif(self.mode==3):
113  rxi = np.stack((pars['fmf'], pars['aot'], pars['wnd']))
114  elif(self.mode==4):
115  rxi = np.stack((pars['fmf'], pars['aot'], pars['chl'], pars['wnd']))
116  model = trp.interpn(self.mpts,rlut,rxi,
117  bounds_error=False,fill_value=None )[0]
118  return (model - data)*scale
119 
120  def process(self,rfl,sza,vza,raa,wnd,chl):
121  cosz = np.cos(sza*D2R)
122  rfl *= cosz/np.pi
123  if(self.mode==0 or self.mode==1):
124  if(self.mode==1):
125  chl, wnd = self.chlc, self.wndc
126  xi = np.stack((chl, wnd, raa, vza, sza))
127  elif(self.mode==2):
128  wnd = self.wndc
129  xi = np.stack((wnd, raa, vza, sza))
130  elif(self.mode==3):
131  chl = self.chlc
132  xi = np.stack((chl, raa, vza, sza))
133  elif(self.mode==4):
134  xi = np.stack((raa, vza, sza))
135  tlut = trp.interpn(self.rpts, self.lut, xi)[0]
136  scale = 1.0/(self.coef*(rfl + 0.00001))
137  mfit = lm.minimize(self.minfun, self.pars, args=(rfl, scale, tlut))
138  self.fmf = mfit.params['fmf'].value
139  self.aot = mfit.params['aot'].value
140  self.sse = mfit.redchi
141  if(self.mode==0 or self.mode==1):
142  self.chl = float(chl)
143  self.wnd = float(wnd)
144  mxi = np.stack((self.fmf,self.aot))
145  elif(self.mode==2):
146  self.wnd = float(wnd)
147  self.chl = mfit.params['chl'].value
148  mxi = np.stack((self.fmf,self.aot,self.chl))
149  elif(self.mode==3):
150  self.chl = float(chl)
151  self.wnd = mfit.params['wnd'].value
152  mxi = np.stack((self.fmf,self.aot,self.wnd))
153  elif(self.mode==4):
154  self.chl = mfit.params['chl'].value
155  self.wnd = mfit.params['wnd'].value
156  mxi = np.stack((self.fmf,self.aot,self.chl,self.wnd))
157  mrfl = trp.interpn(self.mpts, tlut, mxi,
158  bounds_error=False, fill_value=None )[0]
159 # convert return values to unnormalized units
160  self.rfl = rfl*np.pi/cosz
161  self.mrfl = mrfl*np.pi/cosz
162  self.rsd = self.mrfl - self.rfl
163  self.sse = np.dot(self.rsd,self.rsd)/(NWL-2)
164  return self.fmf, self.aot, self.chl, self.wnd, self.sse, self.type
165 
166  def plot(self, iy, ix):
167  plt.clf()
168  plt.grid(True)
169  plt.plot(self.wl_pts, self.rfl, marker='.', color='b', label='measured')
170  plt.plot(self.wl_pts, self.mrfl, marker='.', color='g', label='modeled')
171  plt.plot(self.wl_pts, self.rsd, marker='.', color='r', label='residual')
172  plt.xlabel('wavelength (nm)')
173  plt.ylabel('reflectance')
174  tstr = "dbocean mode {3:d} -- y={4:}, x={5:} aot: {0:.3f} fmf: {1:.3f} sse: {2:.3}"
175  plt.title(tstr.format(self.aot, self.fmf, self.sse, self.mode, iy, ix))
176  plt.legend(loc='upper right')
177  plt.show()
178 
179 
180 class input(object):
181 
182  def __init__(self, l1b_filepath):
183  self.ifile = l1b_filepath
184  print ("Reading sensor data: " + self.ifile)
185  try:
186  self.rfl[W470] = xr.load_dataset(l1b_filepath,group='/observations')['rhot_490']
187  self.rfl[W550] = xr.load_dataset(l1b_filepath,group='/observations')['rhot_550']
188  self.rfl[W659] = xr.load_dataset(l1b_filepath,group='/observations')['rhot_670']
189  self.rfl[W860] = xr.load_dataset(l1b_filepath,group='/observations')['rhot_865']
190  self.rfl[W124] = xr.load_dataset(l1b_filepath,group='/observations')['rhot_1240']
191  self.rfl[W164] = xr.load_dataset(l1b_filepath,group='/observations')['rhot_1610']
192  self.rfl[W213] = xr.load_dataset(l1b_filepath,group='/observations')['rhot_2250']
193  self.sza = xr.load_dataset(l1b_filepath,group='/geolocation')['solar_zenith'].values
194  self.vza = xr.load_dataset(l1b_filepath,group='/geolocation')['sensor_zenith'].values
195  self.raa = xr.load_dataset(l1b_filepath,group='/geolocation')['relative_azimuth'].values
196  self.lat = xr.load_dataset(l1b_filepath,group='/navigation_data')['latitude'].values
197  self.lon = xr.load_dataset(l1b_filepath,group='/navigation_data')['longitude'].values
198  self.cld = xr.load_dataset(l1b_filepath,group='/ancillary')['cloud_mask'].values
199  self.wnd = xr.load_dataset(l1b_filepath,group='/ancillary')['wind_speed'].values
200  self.chl = xr.load_dataset(l1b_filepath,group='/ancillary')['chlorophyll'].values
201  self.rfl = np.transpose(self.rfl, (1,2,0))
202  self.shape = self.sza.shape
203  except Exception as inst:
204  print(type(inst))
205  print(inst)
206  print ("Unable to read from file ... exiting")
207  sys.exit()
208 
209 class output(object):
210 
211  def __init__(self, out_filepath, ydim, xdim):
212  self.ofile = out_filepath
213  self.ydim = ydim
214  self.xdim = xdim
215  try:
216  rfl = xr.DataArray(np.zeros((NWL,ydim,xdim)),dims=('wl','y', 'x'))
217  lat = xr.DataArray(np.zeros((ydim,xdim)),dims=('y', 'x'))
218  lon = xr.DataArray(np.zeros((ydim,xdim)),dims=('y', 'x'))
219  sza = xr.DataArray(np.zeros((ydim,xdim)),dims=('y', 'x'))
220  vza = xr.DataArray(np.zeros((ydim,xdim)),dims=('y', 'x'))
221  raa = xr.DataArray(np.zeros((ydim,xdim)),dims=('y', 'x'))
222  wnd = xr.DataArray(np.zeros((ydim,xdim)),dims=('y', 'x'))
223  chl = xr.DataArray(np.zeros((ydim,xdim)),dims=('y', 'x'))
224  aot = xr.DataArray(np.zeros((ydim,xdim)),dims=('y', 'x'))
225  fmf = xr.DataArray(np.zeros((ydim,xdim)),dims=('y', 'x'))
226  sse = xr.DataArray(np.zeros((ydim,xdim)),dims=('y', 'x'))
227  self.ds = xr.Dataset({'rfl': rfl, 'lat': lat, 'lon': lon,
228  'sza': sza, 'vza': vza, 'raa': raa,
229  'wnd': wnd, 'chl': chl, 'aot': aot,
230  'fmf': fmf, 'sse': sse, })
231  except Exception as inst:
232  print(type(inst))
233  print(inst)
234  print ("Unable to initialize output file ... exiting")
235  sys.exit()
236 
237  def write(self):
238  print ("Writing to file: " + self.ofile)
239  self.ds.to_netcdf(self.ofile)
240 
241 def main():
242 
243  parser = argparse.ArgumentParser()
244  parser.add_argument("-if", "--ifile", type=argparse.FileType('r'),
245  help="input file", required=True)
246  parser.add_argument("-of", "--ofile", type=argparse.FileType('w'),
247  help="output file", required=False)
248  parser.add_argument("-lf", "--lut", type=argparse.FileType('r'),
249  help="lookup table file", required=True)
250  parser.add_argument("-pf", "--plot", type=bool, default=False,
251  help="plot pixel data", required=False)
252  parser.add_argument("-mf", "--mode", type=int, default=0,
253  help="mode option", required=False)
254  parser.add_argument("y", type=int, help="start line")
255  parser.add_argument("x", type=int, help="start pixel")
256  parser.add_argument("z", type=int, help="square side ")
257 
258  args = parser.parse_args()
259  vin = input(args.ifile.name)
260  dimx = dimy = args.z
261  if (args.z <= 0):
262  args.y = args.x = 0
263  dimy = vin.shape[0]
264  dimx = vin.shape[1]
265 
266  vout = output(args.ofile.name, dimy, dimx)
267 
268  dbo = dbocean(args.lut.name, args.mode)
269 
270  print ("Processing mode {0}".format(args.mode))
271  for iy in range(args.y, args.y+dimy):
272  tic = time.perf_counter()
273  for ix in range(args.x, args.x+dimx):
274 
275  if(vin.cld[iy,ix] or vin.chl[iy,ix]<-999):
276  fmf,aot,chl,wnd,sse = -999.9, -999.9, -999.9, -999.9, -999.9
277  else:
278  try:
279  fmf,aot,chl,wnd,sse = dbo.process(vin.rfl[iy,ix],vin.sza[iy,ix],
280  vin.vza[iy,ix],vin.raa[iy,ix],
281  vin.wnd[iy,ix],vin.chl[iy,ix])
282  except Exception as inst:
283  print(type(inst))
284  print(inst)
285  print ("processing error at pixel", iy, ix)
286 
287  if(args.plot):
288  dbo.plot(iy,ix)
289 
290  vout.ds.rfl[:,iy-args.y,ix-args.x] = vin.rfl[iy,ix]
291  vout.ds.lat[iy-args.y,ix-args.x] = vin.lat[iy,ix]
292  vout.ds.lon[iy-args.y,ix-args.x] = vin.lon[iy,ix]
293  vout.ds.sza[iy-args.y,ix-args.x] = vin.sza[iy,ix]
294  vout.ds.vza[iy-args.y,ix-args.x] = vin.vza[iy,ix]
295  vout.ds.raa[iy-args.y,ix-args.x] = vin.raa[iy,ix]
296  vout.ds.wnd[iy-args.y,ix-args.x] = wnd
297  vout.ds.chl[iy-args.y,ix-args.x] = chl
298  vout.ds.aot[iy-args.y,ix-args.x] = aot
299  vout.ds.fmf[iy-args.y,ix-args.x] = fmf
300  vout.ds.sse[iy-args.y,ix-args.x] = sse
301 
302  toc = time.perf_counter()
303  tpp = (toc-tic)/dimx
304 # print('.', end = '', flush=True)
305  print(iy, tpp, flush=True)
306  print(' Done!')
307  vout.write()
308 
309 if __name__ == "__main__":
310 
311  main()
dbocean.input.__init__
def __init__(self, l1b_filepath)
Definition: dbocean.py:182
dbocean.dbocean.plot
def plot(self, iy, ix)
Definition: dbocean.py:166
dbocean.output.__init__
def __init__(self, out_filepath, ydim, xdim)
Definition: dbocean.py:211
dbocean.main
def main()
Definition: dbocean.py:241
names::if
character(len=1000) if
Definition: names.f90:13
dbocean.dbocean.__init__
def __init__(self, lut_filepath, mode)
Definition: dbocean.py:36
dbocean.dbocean.minfun
def minfun(self, pars, data, scale, rlut)
Definition: dbocean.py:107
color_dtdb.range
range
Definition: color_dtdb.py:331
dbocean.dbocean.process
def process(self, rfl, sza, vza, raa, wnd, chl)
Definition: dbocean.py:120
dbocean.output.write
def write(self)
Definition: dbocean.py:237


ABOUT

Our Team

Contact Us

DATA

Get Started

Find Data

Use Data

TOOLS

SeaBASS

SeaDAS

Giovanni Data Visualizer

Worldview Mapping Interface

Level 3 Time Series Plotter

Request a SeaHawk Image

Mission Quality Monitor

HELP

Level 1 & 2 Browser Help

Level 3 & 4 Browser Help

File Search Help

Search & Download Help

Earthdata Forum

Join Mailing List

How to Cite Data

CONNECT WITH US
COMPLIANCE

Web Privacy Policy

Data & Information Policy

Communications Policy

Accessibility Policy

Freedom of Information Act

USA.gov

Responsible NASA Official: Sean Bailey
Curator: OceanColor Webmaster