Due to the lapse in federal government funding, NASA is not updating this website. We sincerely regret this inconvenience.
A sphere with three ocean waves in differing shades of blue. Next to the sphere, there is the bolded text "Ocean Color". Under that, there are two acronyms, separated by a vertical pipe: OB.DAAC (Ocean Biology Distributed Active Archive Center) and OBPG (Ocean Biology Processing Group). Lorem ipsum dolor sit amet
NASA Logo
Ocean Color Science Software

ocssw V2022
dtocean.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 xarray as xr
11 import matplotlib.pyplot as plt
12 import scipy.interpolate as interp
13 from matplotlib import *
14 from matplotlib.ticker import MaxNLocator
15 from mpl_toolkits.mplot3d import Axes3D
16 from datetime import datetime
17 
18 D490 = 0
19 D550 = 1
20 D670 = 2
21 D860 = 3
22 D124 = 4
23 D164 = 5
24 D213 = 6
25 W860 = D860+2
26 NDWL = 7
27 NWS = 4
28 NSZA = 11
29 NVZA = 16
30 NRAA = 16
31 NAOT = 6
32 NBIG = 5
33 NSMALL = 4
34 D2R = np.pi/180.0
35 
36 class dtocean(object):
37 
38  def __init__(self, lut_filepath, mode):
39  print ("Reading Darktarget LUT: " + lut_filepath)
40  self.toas = np.zeros((NDWL,NSMALL,NAOT,NWS,NSZA,NVZA,NRAA))
41  self.toab = np.zeros((NDWL,NBIG,NAOT,NWS,NSZA,NVZA,NRAA))
42  self.toar = np.zeros((NDWL,NWS,NSZA,NVZA,NRAA))
43  self.taus = np.zeros((NDWL,NSMALL,NAOT))
44  self.taub = np.zeros((NDWL,NBIG,NAOT))
45  try:
46  self.toas=xr.load_dataset(lut_filepath,group='/OCEAN_AEROSOL')['AINTS'].values
47  self.toab=xr.load_dataset(lut_filepath,group='/OCEAN_AEROSOL')['AINTB'].values
48  self.toar=xr.load_dataset(lut_filepath,group='/OCEAN_AEROSOL')['REF_RAYALL'].values
49  self.taus=xr.load_dataset(lut_filepath,group='/OCEAN_AEROSOL')['TAUAS'].values
50  self.taub=xr.load_dataset(lut_filepath,group='/OCEAN_AEROSOL')['TAUAB'].values
51  self.wl_pts=xr.load_dataset(lut_filepath,group='/OCEAN_AEROSOL')['WAVE'].values
52  self.raa_pts=xr.load_dataset(lut_filepath,group='/OCEAN_AEROSOL')['PHC'].values
53  self.vza_pts=xr.load_dataset(lut_filepath,group='/OCEAN_AEROSOL')['THET'].values
54  self.sza_pts=xr.load_dataset(lut_filepath,group='/OCEAN_AEROSOL')['THET0'].values
55  self.taus = np.transpose(self.taus,(0,1,2))
56  self.taub = np.transpose(self.taub,(0,1,2))
57  self.toas = np.transpose(self.toas,(3,4,5,6,0,1,2))
58  self.toab = np.transpose(self.toab,(3,4,5,6,0,1,2))
59  self.toar = np.transpose(self.toar,(1,2,3,4,0))
60  self.wnd_pts = np.array([0.0, 6.0, 10.0, 16.0])
61  self.tau_pts = (self.taus[D550,0,:])
62  self.pts = (self.wnd_pts, self.sza_pts, self.vza_pts, self.raa_pts)
63  except Exception as inst:
64  print(type(inst))
65  print(inst)
66  print ("Unable to read LUT file ... exiting")
67  sys.exit()
68 
69  self.tslut = np.zeros((NAOT,NWS,NSZA,NVZA,NRAA,NSMALL))
70  self.tblut = np.zeros((NAOT,NWS,NSZA,NVZA,NRAA,NBIG))
71  self.rfl_pts = np.linspace(0.0,0.2,NAOT)
72  self.tpts = (self.rfl_pts, self.wnd_pts, self.sza_pts, self.vza_pts, self.raa_pts)
73 
74  try:
75  tic = time.perf_counter()
76  maxs = np.max(self.toas)
77  maxb = np.max(self.toab)
78 
79  for iwn in range(NWS):
80  for isz in range(NSZA):
81  for ith in range(NVZA):
82  for iph in range(NRAA):
83  for ism in range(NSMALL):
84  self.tslut[:,iwn,isz,ith,iph,ism] = \
85  np.interp(self.rfl_pts,self.toas[iwn,isz,ith,iph,D860,ism,:],self.tau_pts)
86  for ism in range(NBIG):
87  self.tblut[:,iwn,isz,ith,iph,ism] = \
88  np.interp(self.rfl_pts,self.toab[iwn,isz,ith,iph,D860,ism,:],self.tau_pts)
89 
90  toc = time.perf_counter()
91 
92  except Exception as inst:
93  print(type(inst))
94  print(inst)
95 
96  def proc_gran(self,rfl,sza,vza,raa,wnd):
97  xit = np.stack((wnd, sza, vza, raa),axis=-1)
98  units = np.pi/np.cos(sza*D2R)
99  shp = sza.shape
100 
101  tstart = time.perf_counter()
102  try:
103  tic = time.perf_counter()
104  mrfls = np.multiply(interp.interpn(self.pts, self.toas, xit),np.expand_dims(units,axis=(2,3,4)))
105  mrflb = np.multiply(interp.interpn(self.pts, self.toab, xit),np.expand_dims(units,axis=(2,3,4)))
106  mrflr = np.multiply(interp.interpn(self.pts, self.toar, xit),np.expand_dims(units,axis=2))
107  toc = time.perf_counter()
108  except Exception as inst:
109  print(type(inst))
110  print(inst)
111  print('\nGeometry interpolation: ', toc-tic, 'sec, ', (toc-tic)/shp[0]/shp[1], "per pixel")
112 
113  tauxs = np.zeros((shp[0],shp[1],NSMALL))
114  tauxb = np.zeros((shp[0],shp[1],NBIG))
115  trfls = np.zeros((shp[0],shp[1],NDWL,NSMALL))
116  trflb = np.zeros((shp[0],shp[1],NDWL,NBIG))
117  AA = np.zeros((shp[0],shp[1],NSMALL,NBIG,7,2))
118  tic = time.perf_counter()
119  for ism in range(NSMALL):
120  mrfls[:,:,:,ism,0] = mrflr
121  for ibm in range(NBIG):
122  mrflb[:,:,:,ibm,0] = mrflr
123  for iy in range(shp[0]):
124  for ix in range(shp[1]):
125 
126  for ism in range(NSMALL):
127  tauxs[iy,ix,ism] = np.interp(rfl[iy,ix,W860], mrfls[iy,ix,D860,ism], self.taus[D550,ism])
128  for ibm in range(NBIG):
129  tauxb[iy,ix,ibm] = np.interp(rfl[iy,ix,W860], mrflb[iy,ix,D860,ibm], self.taub[D550,ibm])
130 
131  for iwl in range(NDWL):
132  for ism in range(NSMALL):
133  trfls[iy,ix,iwl,ism] = np.interp(tauxs[iy,ix,ism],self.taus[D550,ism],mrfls[iy,ix,iwl,ism])
134  for ibm in range(NBIG):
135  trflb[iy,ix,iwl,ibm] = np.interp(tauxb[iy,ix,ibm],self.taub[D550,ibm],mrflb[iy,ix,iwl,ibm])
136 
137  for ism in range(NSMALL):
138  for ibm in range(NBIG):
139  AA[iy,ix,ism,ibm] = np.vstack((trfls[iy,ix,:,ism],trflb[iy,ix,:,ibm])).T
140 
141  toc = time.perf_counter()
142  print('Per-pixel interpolation: ', toc-tic, 'sec, ', (toc-tic)/shp[0]/shp[1], "per pixel")
143 
144  tic = time.perf_counter()
145  rfle = np.expand_dims(rfl[2:],axis=(2,3,5))
146  try:
147  xm = np.clip(np.matmul(np.linalg.pinv(AA),rfle),0,1)[:,:,:,:,0]
148  mrfl = np.multiply(xm[:,:,:,:],AA[:,:,:,:,:,0]) + np.multiply((1-xm[:,:,:,:]),AA[:,:,:,:,:,1])
149  err = np.squeeze(rfle,axis=5)-mrfl
150  sse = np.linalg.norm(err[:,:,:,:,1:],axis=4)**2/(NDWL-2)
151  except Exception as inst:
152  print(type(inst))
153  print(inst)
154  toc = time.perf_counter()
155  print('FMF calculation: ', toc-tic, 'sec, ', (toc-tic)/shp[0]/shp[1], "per pixel")
156  tstop = time.perf_counter()
157  print('Total processing time: ', shp[0], 'lines, ', tstop-tstart, 'sec, ', (tstop-tstart)/shp[0]/shp[1], "per pixel")
158 
159  rmin = sse.reshape(shp[0],shp[1],NSMALL*NBIG).argmin(axis=2)
160  im = divmod(rmin,NBIG)
161  self.type = (NBIG-1)*im[0] + im[1]
162  Y,X = np.ogrid[:shp[0],:shp[1]]
163  self.sse = sse.reshape(shp[0],shp[1],NSMALL*NBIG)[Y,X,rmin]
164  self.fmf = xm.reshape(shp[0],shp[1],NSMALL*NBIG)[Y,X,rmin]
165  self.aot = np.multiply(self.fmf,tauxs[Y,X,im[0]]) + np.multiply((1.0-self.fmf),tauxb[Y,X,im[1]])
166 
167  return self.fmf, self.aot, self.sse, self.type
168 
169  def process(self,rfl,sza,vza,raa,wnd):
170  xit = np.stack((wnd, sza, vza, raa))
171  units = np.pi/np.cos(sza*D2R)
172  rfld = rfl[2:]
173  mrfls = interp.interpn(self.pts, self.toas, xit, bounds_error=False, fill_value=None)[0]*units
174  mrflb = interp.interpn(self.pts, self.toab, xit, bounds_error=False, fill_value=None)[0]*units
175  mrflr = interp.interpn(self.pts, self.toar, xit, bounds_error=False, fill_value=None)[0]*units
176 
177  tauxs = np.zeros(NSMALL)
178  tauxb = np.zeros(NBIG)
179  for ism in range(NSMALL):
180  mrfls[:,ism,0] = mrflr
181  tauxs[ism] = interp.interpn((mrfls[D860,ism,:],), self.taus[D550,ism,:], [rfld[D860]], bounds_error=False, fill_value=None)
182  for ibm in range(NBIG):
183  mrflb[:,ibm,0] = mrflr
184  tauxb[ibm] = interp.interpn((mrflb[D860,ibm,:],), self.taub[D550,ibm,:], [rfld[D860]], bounds_error=False, fill_value=None)
185 
186  trfls = np.zeros((NDWL,NSMALL))
187  trflb = np.zeros((NDWL,NBIG))
188  for iwl in range(NDWL):
189  for ism in range(NSMALL):
190  trfls[iwl,ism] = interp.interpn((self.taus[D550,ism],), mrfls[iwl,ism], [tauxs[ism]], bounds_error=False, fill_value=None)
191  for ibm in range(NBIG):
192  trflb[iwl,ibm] = interp.interpn((self.taub[D550,ibm],), mrflb[iwl,ibm], [tauxb[ibm]], bounds_error=False, fill_value=None)
193 
194  AA = np.zeros((NSMALL,NBIG,7,2))
195  for ism in range(NSMALL):
196  for ibm in range(NBIG):
197  AA[ism,ibm] = np.vstack((trfls[:,ism],trflb[:,ibm])).T
198 
199  try:
200  xm = np.clip(np.dot(np.linalg.pinv(AA),rfld)[:,:,0],0,1)
201  xmd = np.expand_dims(xm,axis=2)
202  mrfl = xmd*AA[:,:,:,0] + (1-xmd)*AA[:,:,:,1]
203  err = rfld-mrfl
204  sse = np.linalg.norm(err[:,:,1:],axis=2)**2/(NDWL-2)
205  except Exception as inst:
206  print(type(inst))
207  print(inst)
208 
209  im = divmod(sse.argmin(),NBIG)
210  self.type = (NBIG-1)*im[0] + im[1]
211  self.sse = sse[im[0],im[1]]
212  self.fmf = xm[im[0],im[1]]
213  self.aot = self.fmf*tauxs[im[0]] + (1.0-self.fmf)*tauxb[im[1]]
214  self.trfls = trfls
215  self.trflb = trflb
216  self.srfl = trfls[:,im[0]]
217  self.brfl = trflb[:,im[1]]
218  self.mrfl = self.fmf*trfls[:,im[0]] + (1.0-self.fmf)*trflb[:,im[1]]
219  self.rsd = (rfld - self.mrfl)/rfld
220  self.rfl = rfld
221  self.sse = np.dot(self.rsd,self.rsd)/(NDWL-2)
222 
223  return self.fmf, self.aot, wnd, self.sse, self.type
224 
225  def plot(self, iy, ix):
226  plt.clf()
227  plt.grid(True)
228  plt.plot(self.wl_pts, self.srfl, marker='.', color='b', label='fine mode')
229  plt.plot(self.wl_pts, self.brfl, marker='.', color='k', label='coarse mode')
230  plt.plot(self.wl_pts, self.trfls, marker='', color='b', label='fine mode', linewidth=1)
231  plt.plot(self.wl_pts, self.trflb, marker='', color='k', label='coarse mode', linewidth=1)
232  plt.plot(self.wl_pts, self.rfl, marker='.', color='r', label='measured')
233  plt.plot(self.wl_pts, self.mrfl, marker='.', color='g', label='modeled')
234 # plt.plot(self.wl_pts, self.rsd, marker='.', color='r', label='residual')
235  plt.xlabel('wavelength (nm)')
236  plt.ylabel('reflectance')
237  tstr = "dtocean -- y={3:}, x={4:} aot: {0:.3f} fmf: {1:.3f} sse: {2:.3}"
238  plt.title(tstr.format(self.aot, self.fmf, self.sse, iy, ix))
239  plt.legend(loc='upper right')
240  plt.show()
241 
242  def plot_points(self):
243  fig = plt.figure()
244  ax = fig.add_subplot(111, projection='3d')
245  small = [0,1,2,3]
246  for s in small:
247  xs = toast[D550,s,:]
248  ys = toast[D670,s,:]
249  zs = toast[D860,s,:]
250  ax.plot(xs, ys, zs, c="b", label="small model")
251  ax.scatter(xs, ys, zs, c="b", marker="o")
252  big = [0,1,2,3,5]
253  for b in big:
254  xb = toabt[D550,b,:]
255  yb = toabt[D670,b,:]
256  zb = toabt[D860,b,:]
257  ax.plot(xb, yb, zb, c="r", label="big model")
258  ax.scatter(xb, yb, zb, c="r", marker="o")
259  for s in small:
260  for b in big:
261  xp = bspair[:,D550,s,b]
262  yp = bspair[:,D670,s,b]
263  zp = bspair[:,D860,s,b]
264  ax.plot(xp, yp, zp, c="g", label="big/small continuum")
265 
266  xs = trfls[D550,:]
267  ys = trfls[D670,:]
268  zs = trfls[D860,:]
269  ax.scatter(xs, ys, zs, c="b", marker="o", s=50)
270  xb = trflb[D550,:]
271  yb = trflb[D670,:]
272  zb = trflb[D860,:]
273  ax.scatter(xb, yb, zb, c="r", marker="o", s=50)
274  xm = rfl[W550,row,col]
275  ym = rfl[W670,row,col]
276  zm = rfl[W860,row,col]
277  ax.scatter(xm, ym, zm, c="k", marker="o", s=50)
278  ax.set_xlabel('D550')
279  ax.set_ylabel('D670')
280  ax.set_zlabel('D860')
281  #plt.show()
282 
283  fig = plt.figure()
284  ax = fig.add_subplot(111, projection='3d')
285  ws = np.arange(NDWL)
286  ts = np.arange(NAOT+1)
287  ws, ts = np.meshgrid(ws, ts)
288  ls = np.zeros((NDWL,NAOT+1))
289  lb = np.zeros((NDWL,NAOT+1))
290  ls[:,:-1] = toast[:,s,:]
291  ls[:,NAOT] = ls[:,NAOT-1]
292  ax.plot_wireframe(ws, ts, ls, color='b', label="small model")
293  lb[:,:-1] = toabt[:,b,:]
294  lb[:,NAOT] = lb[:,NAOT-1]
295  ax.plot_wireframe(ws, ts, lb, color='r', label="big model")
296  plt.show()
297 
dtocean.dtocean.plot
def plot(self, iy, ix)
Definition: dtocean.py:225
dtocean.dtocean.__init__
def __init__(self, lut_filepath, mode)
Definition: dtocean.py:38
dtocean.dtocean.plot_points
def plot_points(self)
Definition: dtocean.py:242
print
void print(std::ostream &stream, const char *format)
Definition: PrintDebug.hpp:38
dtocean.dtocean.proc_gran
def proc_gran(self, rfl, sza, vza, raa, wnd)
Definition: dtocean.py:96
dtocean.dtocean.process
def process(self, rfl, sza, vza, raa, wnd)
Definition: dtocean.py:169
color_dtdb.range
range
Definition: color_dtdb.py:331


ABOUT

Our Team

Contact Us

DATA

Get Started

Find Data

Use Data

TOOLS

SeaBASS

SeaDAS

Giovanni Data Visualizer

Worldview Mapping Interface

Mission Quality Monitor

HELP

How to Use Earthdata Search

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

Freedom of Information Act

USA.gov

Responsible NASA Official: Sean Bailey
Curator: OceanColor Webmaster