OB.DAAC Logo
NASA Logo
Ocean Color Science Software

ocssw V2022
dtland.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 argparse
8 import numpy as np
9 import xarray as xr
10 import matplotlib.pyplot as plt
11 import scipy.interpolate as interp
12 from matplotlib import *
13 from matplotlib.ticker import MaxNLocator
14 from mpl_toolkits.mplot3d import Axes3D
15 from datetime import datetime
16 import Tasmanian
17 
18 W470 = 0
19 W550 = 1
20 W659 = 2
21 W860 = 3
22 W124 = 4
23 W164 = 5
24 W213 = 6
25 
26 D213 = 3
27 slope_644 = 0.559
28 yint_644 = 0.0
29 slope_466 = 0.645
30 yint_466 = 0.0
31 
32 ISMALL = 0
33 IBIG = 1
34 NSEASONS = 4
35 NLATS = 180
36 NLONS = 360
37 NLWL = 4
38 NWL =7
39 NSZA = 11
40 NVZA = 15
41 NRAA = 16
42 NTAU = 7
43 NAOT = 6
44 NTAB = 5
45 PRESSURE_P0 = 1013.0
46 D2R = np.pi/180.0
47 
48 class dtland(object):
49 
50  def __init__(self, lut_filepath, mode):
51  print ("Reading Darktarget LUT: " + lut_filepath)
52  self.all = np.zeros((NSEASONS,NLATS,NLONS))
53  self.int = np.zeros((NTAB,NLWL,NTAU,NSZA,NVZA,NRAA))
54  self.t = np.zeros((NTAB,NLWL,NTAU,NSZA,NVZA))
55  self.fd = np.zeros((NTAB,NLWL,NTAU,NSZA))
56  self.sbar = np.zeros((NTAB,NLWL,NTAU,NSZA))
57  self.opth = np.zeros((NTAB,NLWL,NTAU))
58  self.massc = np.zeros((NTAB,NLWL,NTAU))
59  self.extn = np.zeros((NTAB,NLWL,NTAU))
60  self.ssa = np.zeros((NTAB,NLWL,NTAU))
61  self.qext = np.zeros((NTAB,NLWL,NTAU))
62  self.bext = np.zeros((NTAB,NLWL,NTAU))
63  self.vext = np.zeros((NTAB,NLWL,NTAU))
64  self.mext = np.zeros((NTAB,NLWL,NTAU))
65  self.opth = np.zeros((NTAB,NLWL,NTAU))
66  try:
67  self.wl_pts=xr.load_dataset(lut_filepath,group='/LAND AEROSOL')['WAV_NL'].values
68  self.raa_pts=xr.load_dataset(lut_filepath,group='/LAND AEROSOL')['PHI_NL'].values
69  self.vza_pts=xr.load_dataset(lut_filepath,group='/LAND AEROSOL')['THE_NL'].values
70  self.sza_pts=xr.load_dataset(lut_filepath,group='/LAND AEROSOL')['THET0_NL'].values
71  self.mu0_pts=xr.load_dataset(lut_filepath,group='/LAND AEROSOL')['MU0_NL'].values
72  self.all=xr.load_dataset(lut_filepath,group='/LAND AEROSOL')['AEROSOL_ALL'].values
73  self.opth=xr.load_dataset(lut_filepath,group='/LAND AEROSOL')['OPTH_NL0'].values
74  self.int=xr.load_dataset(lut_filepath,group='/LAND AEROSOL')['INT_NL0'].values
75  self.t=xr.load_dataset(lut_filepath,group='/LAND AEROSOL')['T_NL0'].values
76  self.fd=xr.load_dataset(lut_filepath,group='/LAND AEROSOL')['Fd_NL0'].values
77  self.sbar=xr.load_dataset(lut_filepath,group='/LAND AEROSOL')['SBAR_NL0'].values
78  self.massc=xr.load_dataset(lut_filepath,group='/LAND AEROSOL')['MASSCOEF_NL0'].values
79  self.extn=xr.load_dataset(lut_filepath,group='/LAND AEROSOL')['EXTNORM_NL0'].values
80  self.ssa=xr.load_dataset(lut_filepath,group='/LAND AEROSOL')['SSA_NL0'].values
81  self.qext=xr.load_dataset(lut_filepath,group='/LAND AEROSOL')['QEXT_NL0'].values
82  self.bext=xr.load_dataset(lut_filepath,group='/LAND AEROSOL')['BEXT_NL0'].values
83  self.vext=xr.load_dataset(lut_filepath,group='/LAND AEROSOL')['VEXT_NL0'].values
84  self.mext=xr.load_dataset(lut_filepath,group='/LAND AEROSOL')['MEXT_NL0'].values
85  except:
86  print ("Unable to read LUT file ... exiting")
87  sys.exit()
88 
89  self.int = np.transpose(self.int, (3,4,5,1,0,2))
90  self.t = np.transpose(self.t, (3,4,1,0,2))
91  self.fd = np.transpose(self.fd, (3,1,0,2))
92  self.sbar = np.transpose(self.sbar, (3,1,0,2))
93  self.pts3 = (self.sza_pts, self.vza_pts, self.raa_pts)
94  self.pts2 = (self.sza_pts, self.vza_pts)
95  self.pts1 = (self.sza_pts,)
96 
97  def process(self,rfl,sza,vza,raa,height,mtab):
98  xit3 = np.stack((sza, vza, raa))
99  xit2 = np.stack((sza, vza))
100  xit1 = [[sza]]
101 # interpolate angles
102  int_nl = interp.interpn(self.pts3, self.int, xit3)[0]
103  t_nl = interp.interpn(self.pts2, self.t, xit2)[0]
104  fd_nl = interp.interpn(self.pts1, self.fd, xit1)[0]
105  sbar_nl = interp.interpn(self.pts1, self.sbar, xit1)[0]
106  fdt_nl = np.multiply(t_nl,fd_nl)
107 # interpolate elevation
108  eqwav_nl = np.zeros((NLWL))
109  pressure = np.exp(-(height/7.5))
110  for iw in range(3):
111  wl = self.wl_pts[iw]
112  expfactor = -4.15 + (0.2*wl)
113  rod_pres = 0.0088*pressure*np.power(wl, expfactor)
114  lambda1 = 0.1
115  lambda2 = 2.0
116  diff0 = 99.0
117  while diff0 > 0.00001:
118  lambda0 = (lambda1 + lambda2) / 2.0
119  ftau0 = 0.0088*np.power(lambda0,-4.15+0.2*lambda0)
120  ftau1 = 0.0088*np.power(lambda1,-4.15+0.2*lambda1)
121  ftau2 = 0.0088*np.power(lambda2,-4.15+0.2*lambda2)
122  if ((ftau1 > rod_pres) and (ftau2 < rod_pres)):
123  if (ftau0 > rod_pres):
124  lambda1 = (lambda1 + lambda2)/2.0
125  else:
126  lambda2 = (lambda1 + lambda2)/2.0
127  diff0 = np.abs(ftau0 - rod_pres)
128  eqwav_nl[iw] = np.log(lambda0);
129  logwl = (np.log(self.wl_pts),)
130  int_nl[:-1] = np.exp(interp.interpn(logwl, np.log(int_nl), \
131  (eqwav_nl[:-1],), bounds_error=False, fill_value=None))
132  t_nl[:-1] = np.exp(interp.interpn(logwl, np.log(t_nl), \
133  (eqwav_nl[:-1],), bounds_error=False, fill_value=None))
134  fd_nl[:-1] = np.exp(interp.interpn(logwl, np.log(fd_nl), \
135  (eqwav_nl[:-1],), bounds_error=False, fill_value=None))
136  sbar_nl[:-1] = np.exp(interp.interpn(logwl, np.log(sbar_nl), \
137  (eqwav_nl[:-1],), bounds_error=False, fill_value=None))
138  fdt_nl[:-1] = np.exp(interp.interpn(logwl, np.log(fdt_nl), \
139  (eqwav_nl[:-1],), bounds_error=False, fill_value=None))
140 # simulate toa
141  rho = np.zeros((NLWL,NTAB,NTAU))
142  rm213 = np.ones((NTAB,NTAU))*rfl[W213]
143  rho[D213] = np.divide(int_nl[D213]-rm213, \
144  np.multiply(sbar_nl[D213],(int_nl[D213]-rm213))-fdt_nl[D213])
145  rho[W659] = slope_644*rho[D213] + yint_644
146  rho[W470] = slope_466*rho[W659] + yint_466
147  rho[W550].fill(0)
148  rho_star = int_nl + np.divide(np.multiply(fdt_nl,rho), \
149  (np.ones_like(rho)-np.multiply(sbar_nl,rho)))
150 
151  tauxs = np.zeros(NSMALL)
152  tauxb = np.zeros(NBIG)
153  for ism in range(NSMALL):
154  tauxs[ism] = np.interp(rfl[W860], mrfls[W860,ism], self.taus[W550,ism])
155  for ibm in range(NBIG):
156  tauxb[ibm] = np.interp(rfl[W860], mrflb[W860,ibm], self.taub[W550,ibm])
157 
158  trfls = np.zeros((NWL,NSMALL))
159  trflb = np.zeros((NWL,NBIG))
160  for iwl in range(NWL):
161  for ism in range(NSMALL):
162  trfls[iwl,ism] = np.interp(tauxs[ism],self.taus[W550,ism],mrfls[iwl,ism])
163  for ibm in range(NBIG):
164  trflb[iwl,ibm] = np.interp(tauxb[ibm],self.taub[W550,ibm],mrflb[iwl,ibm])
165 
166  aot = np.zeros((NSMALL,NBIG))
167  fmf = np.zeros((NSMALL,NBIG))
168  sse = np.zeros((NSMALL,NBIG))
169  for ism in range(NSMALL):
170  for ibm in range(NBIG):
171  denom = rfl - mrflr + 0.01
172  rbdif = np.divide((rfl - trflb[:,ibm]),denom)
173  sbdif = np.divide((trfls[:,ism] - trflb[:,ibm]),denom)
174  xm = np.dot(rbdif,sbdif)/np.dot(sbdif,sbdif)
175  xm = np.max((np.min((xm,1.0)),0.0))
176  mrfl = xm*trfls[:,ism] + (1.0-xm)*trflb[:,ibm]
177  err = np.divide((rfl - mrfl),denom)
178  sse[ism,ibm] = np.dot(err,err)
179  aot[ism,ibm] = xm*tauxs[ism] + (1.0-xm)*tauxb[ibm]
180  fmf[ism,ibm] = xm
181  im = divmod(sse.argmin(), sse.shape[1])
182  self.sse = sse[im[0],im[1]]
183  self.aot = aot[im[0],im[1]]
184  self.fmf = fmf[im[0],im[1]]
185  self.rfl = rfl
186  self.mrfl = self.fmf*trfls[:,im[0]] + (1.0-self.fmf)*trflb[:,im[1]]
187  self.rsd = self.mrfl - self.rfl
188  self.sse = np.dot(self.rsd,self.rsd)/(NWL-2)
189  self.rayl = mrflr
190 
191  return self.fmf, self.aot, self.sse
192 
193  def plot(self, iy, ix):
194  plt.clf()
195  plt.grid(True)
196  plt.plot(self.wl_pts, self.rfl, marker='.', color='b', label='measured')
197  plt.plot(self.wl_pts, self.mrfl, marker='.', color='g', label='modeled')
198  plt.plot(self.wl_pts, self.rsd, marker='.', color='r', label='residual')
199  plt.xlabel('wavelength (nm)')
200  plt.ylabel('reflectance')
201  tstr = "dtland -- y={3:}, x={4:} aot: {0:.3f} fmf: {1:.3f} sse: {2:.3}"
202  plt.title(tstr.format(self.aot, self.fmf, self.sse, iy, ix))
203  plt.legend(loc='upper right')
204  plt.show()
205 
206  def plot_points(self):
207  fig = plt.figure()
208  ax = fig.add_subplot(111, projection='3d')
209  for s in small:
210  xs = toast[W550,s,:]
211  ys = toast[W659,s,:]
212  zs = toast[W860,s,:]
213  ax.plot(xs, ys, zs, c="b", label="small model")
214  ax.scatter(xs, ys, zs, c="b", marker="o")
215  for b in big:
216  xb = toabt[W550,b,:]
217  yb = toabt[W659,b,:]
218  zb = toabt[W860,b,:]
219  ax.plot(xb, yb, zb, c="r", label="big model")
220  ax.scatter(xb, yb, zb, c="r", marker="o")
221  for s in small:
222  for b in big:
223  xp = bspair[:,W550,s,b]
224  yp = bspair[:,W659,s,b]
225  zp = bspair[:,W860,s,b]
226  ax.plot(xp, yp, zp, c="g", label="big/small continuum")
227 
228  xs = trfls[W550,:]
229  ys = trfls[W659,:]
230  zs = trfls[W860,:]
231  ax.scatter(xs, ys, zs, c="b", marker="o", s=50)
232  xb = trflb[W550,:]
233  yb = trflb[W659,:]
234  zb = trflb[W860,:]
235  ax.scatter(xb, yb, zb, c="r", marker="o", s=50)
236  xm = rfl[W550,row,col]
237  ym = rfl[W659,row,col]
238  zm = rfl[W860,row,col]
239  ax.scatter(xm, ym, zm, c="k", marker="o", s=50)
240  ax.set_xlabel('W550')
241  ax.set_ylabel('W659')
242  ax.set_zlabel('W860')
243  #plt.show()
244 
245  fig = plt.figure()
246  ax = fig.add_subplot(111, projection='3d')
247  ws = np.arange(NWL)
248  ts = np.arange(NAOT+1)
249  ws, ts = np.meshgrid(ws, ts)
250  ls = np.zeros((NWL,NAOT+1))
251  lb = np.zeros((NWL,NAOT+1))
252  ls[:,:-1] = toast[:,s,:]
253  ls[:,NAOT] = ls[:,NAOT-1]
254  ax.plot_wireframe(ws, ts, ls, color='b', label="small model")
255  lb[:,:-1] = toabt[:,b,:]
256  lb[:,NAOT] = lb[:,NAOT-1]
257  ax.plot_wireframe(ws, ts, lb, color='r', label="big model")
258  plt.show()
259 
260 class input(object):
261 
262  def __init__(self, l1b_filepath):
263  self.ifile = l1b_filepath
264  print ("Reading VIIRS Data: " + self.ifile)
265  try:
266  self.rfl = np.append(xr.load_dataset(l1b_filepath,group='/reflectance')['toa_reflectance'][1:6,:,:].values,
267  xr.load_dataset(l1b_filepath,group='/reflectance')['toa_reflectance'][7:,:,:].values,axis=0)
268  self.sza = xr.load_dataset(l1b_filepath,group='/geolocation')['solar_zenith'].values
269  self.vza = xr.load_dataset(l1b_filepath,group='/geolocation')['sensor_zenith'].values
270  self.raa = xr.load_dataset(l1b_filepath,group='/geolocation')['relative_azimuth'].values
271  self.height = xr.load_dataset(l1b_filepath,group='/geolocation')['elevation'].values
272  self.lat = xr.load_dataset(l1b_filepath,group='/navigation_data')['latitude'].values
273  self.lon = xr.load_dataset(l1b_filepath,group='/navigation_data')['longitude'].values
274  self.cld = xr.load_dataset(l1b_filepath,group='/ancillary')['cloud_mask'].values
275  self.wnd = xr.load_dataset(l1b_filepath,group='/ancillary')['wind_speed'].values
276  self.rfl = np.transpose(self.rfl, (1,2,0))
277  except:
278  print ("Unable to read from input file ... exiting")
279  sys.exit()
280 
281 
282 class output(object):
283 
284  def __init__(self, out_filepath, ydim, xdim):
285  self.ofile = out_filepath
286  self.ydim = ydim
287  self.xdim = xdim
288  try:
289  self.rfl = xr.DataArray(np.zeros((NWL,ydim,xdim)),dims=('wl','y', 'x'))
290  self.lat = xr.DataArray(np.zeros((ydim,xdim)),dims=('y', 'x'))
291  self.lon = xr.DataArray(np.zeros((ydim,xdim)),dims=('y', 'x'))
292  self.sza = xr.DataArray(np.zeros((ydim,xdim)),dims=('y', 'x'))
293  self.vza = xr.DataArray(np.zeros((ydim,xdim)),dims=('y', 'x'))
294  self.raa = xr.DataArray(np.zeros((ydim,xdim)),dims=('y', 'x'))
295  self.wnd = xr.DataArray(np.zeros((ydim,xdim)),dims=('y', 'x'))
296  self.aot = xr.DataArray(np.zeros((ydim,xdim)),dims=('y', 'x'))
297  self.fmf = xr.DataArray(np.zeros((ydim,xdim)),dims=('y', 'x'))
298  self.sse = xr.DataArray(np.zeros((ydim,xdim)),dims=('y', 'x'))
299  self.ds = xr.Dataset({'rfl': self.rfl, 'lat': self.lat, 'lon': self.lon,
300  'sza': self.sza, 'vza': self.vza, 'raa': self.raa,
301  'wnd': self.wnd, 'aot': self.aot,
302  'fmf': self.fmf, 'sse': self.sse })
303  except:
304  print ("Unable to initialize output file ... exiting")
305  sys.exit()
306 
307  def write(self):
308  print ("Writing to file: " + self.ofile)
309  self.ds.to_netcdf(self.ofile)
310 
311 def main():
312 
313  parser = argparse.ArgumentParser()
314  parser.add_argument("-i", "--ifile", type=argparse.FileType('r'),
315  help="input file", required=True)
316  parser.add_argument("-o", "--ofile", type=argparse.FileType('w'),
317  help="output file", required=False)
318  parser.add_argument("-l", "--lut", type=argparse.FileType('r'),
319  help="lookup table file", required=True)
320  parser.add_argument("-p", "--plot", type=bool, default=False,
321  help="plot pixel data", required=False)
322  parser.add_argument("-m", "--mode", type=int, default=0,
323  help="mode option", required=False)
324  parser.add_argument("y", type=int, help="start line")
325  parser.add_argument("x", type=int, help="start pixel")
326  parser.add_argument("z", type=int, help="square side ")
327 
328  args = parser.parse_args()
329  vin = input(args.ifile.name)
330  vout = output(args.ofile.name, args.z, args.z)
331  dtl = dtland(args.lut.name, args.mode)
332 
333  print ("Processing mode {0}".format(args.mode))
334  for iy in range(args.y, args.y+args.z):
335  for ix in range(args.x, args.x+args.z):
336 
337  ilon = 180 + np.round(vin.lon + 0.5)
338  ilat = 90 - np.round(vin.lat + 0.5)
339 # mtab = np.array([int(dtl.all[1,ilat,ilon]+1), 4])
340  mtab = np.array([1, 4])
341 
342  if(vin.cld[iy,ix]):
343  fmf,aot,sse = -999.9, -999.9, -999.9
344  else:
345  fmf,aot,sse = dtl.process(vin.rfl[iy,ix,:],vin.sza[iy,ix],
346  vin.vza[iy,ix],vin.raa[iy,ix],
347  vin.height[iy,ix]/1000.0,mtab)
348  if(args.plot):
349  dtl.plot(iy,ix)
350 
351  vout.ds.rfl[:,iy-args.y,ix-args.x] = vin.rfl[iy,ix]
352  vout.ds.lat[iy-args.y,ix-args.x] = vin.lat[iy,ix]
353  vout.ds.lon[iy-args.y,ix-args.x] = vin.lon[iy,ix]
354  vout.ds.sza[iy-args.y,ix-args.x] = vin.sza[iy,ix]
355  vout.ds.vza[iy-args.y,ix-args.x] = vin.vza[iy,ix]
356  vout.ds.raa[iy-args.y,ix-args.x] = vin.raa[iy,ix]
357  vout.ds.wnd[iy-args.y,ix-args.x] = vin.wnd[iy,ix]
358  vout.ds.aot[iy-args.y,ix-args.x] = aot
359  vout.ds.fmf[iy-args.y,ix-args.x] = fmf
360  vout.ds.sse[iy-args.y,ix-args.x] = sse
361 
362  print('.', end = '')
363  print(' Done!')
364  vout.write()
365 
366 if __name__ == "__main__":
367 
368  main()
dtland.dtland.plot
def plot(self, iy, ix)
Definition: dtland.py:193
dtland.dtland.__init__
def __init__(self, lut_filepath, mode)
Definition: dtland.py:50
dtland.main
def main()
Definition: dtland.py:311
dtland.input.__init__
def __init__(self, l1b_filepath)
Definition: dtland.py:262
names::if
character(len=1000) if
Definition: names.f90:13
dtland.output.__init__
def __init__(self, out_filepath, ydim, xdim)
Definition: dtland.py:284
dtland.dtland.process
def process(self, rfl, sza, vza, raa, height, mtab)
Definition: dtland.py:97
dtland.output.write
def write(self)
Definition: dtland.py:307
dtland.dtland.plot_points
def plot_points(self)
Definition: dtland.py:206
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

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