ocssw V2020
color_dtdb.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 datetime
7 import math
8 from PIL import Image as im
9 from PIL import ImageDraw
10 from PIL import ImageFilter
11 from PIL import ImageMath
12 from PIL import ImageEnhance
13 from PIL import ImageOps as imo
14 from PIL import ImageColor as imc
15 
16 import viirs_files as vf
17 import numpy as np
18 import xarray as xr
19 import matplotlib as mpl
20 import matplotlib.pyplot as plt
21 from matplotlib import cm
22 
23 Kb = 0.114
24 Kr = 0.299
25 
26 RGB2YCbCr = np.array([[0.257, 0.504, 0.098],[-0.148, -0.291, 0.439],[0.439, -0.368, -0.071]])
27 YCbCr2RGB = np.array([[1.0, 0.0, 1.596],[1.0, -0.392, -0.813],[1.0, 2.017, 0.0]])
28 
29 def rgb2ycbcr(r,g,b):
30  y = RGB2YCbCr[0,0]*r + RGB2YCbCr[0,1]*g + RGB2YCbCr[0,2]*b
31  u = RGB2YCbCr[1,0]*r + RGB2YCbCr[1,1]*g + RGB2YCbCr[1,2]*b
32  v = RGB2YCbCr[2,0]*r + RGB2YCbCr[2,1]*g + RGB2YCbCr[2,2]*b
33  return y,u,v
34 
35 def ycbcr2rgb(y,cb,cr):
36  r = YCbCr2RGB[0,0]*y + YCbCr2RGB[0,1]*cb + YCbCr2RGB[0,2]*cr
37  g = YCbCr2RGB[1,0]*y + YCbCr2RGB[1,1]*cb + YCbCr2RGB[1,2]*cr
38  b = YCbCr2RGB[2,0]*y + YCbCr2RGB[2,1]*cb + YCbCr2RGB[2,2]*cr
39  return r,g,b
40 
41 
42 RGB2YUV = np.array([[0.299, 0.587, 0.114],[-0.14713, -0.28886, 0.436],[0.615, -0.51499, -0.10001]])
43 YUV2RGB = np.array([[1.0, 0.0, 1.13983],[1.0, -0.39465, -0.58060],[1.0, 2.03211, 0.0]])
44 
45 def rgb2yuv(r,g,b):
46  y = (RGB2YUV[0,0]*r + RGB2YUV[0,1]*g + RGB2YUV[0,2]*b)
47  u = (RGB2YUV[1,0]*r + RGB2YUV[1,1]*g + RGB2YUV[1,2]*b)
48  v = (RGB2YUV[2,0]*r + RGB2YUV[2,1]*g + RGB2YUV[2,2]*b)
49  return y,u,v
50 
51 def yuv2rgb(y,u,v):
52  r = (YUV2RGB[0,0]*y + YUV2RGB[0,1]*u + YUV2RGB[0,2]*v)
53  g = (YUV2RGB[1,0]*y + YUV2RGB[1,1]*u + YUV2RGB[1,2]*v)
54  b = (YUV2RGB[2,0]*y + YUV2RGB[2,1]*u + YUV2RGB[2,2]*v)
55  return r,g,b
56 
57 
58 def rgb2ycbcr2(r,g,b):
59  y = Kr*r + (1-Kr-Kb)*g + Kb*b
60  pb = (b-y)/(1-Kb)/2
61  pr = (r-y)/(1-Kr)/2
62  return y,pb,pr
63 
64 def ycbcr2rgb2(y,pb,pr):
65  r = 2*pr*(1-Kr) + y
66  b = 2*pb*(1-Kb) + y
67  g = (y - Kr*r - Kb*b)/(1-Kr-Kb)
68  return r,g,b
69 
70 
71 def histeq(im,nbr_bins=256):
72 
73  #get image histogram
74  imhist,bins = np.histogram(im.flatten(),nbr_bins)
75 
76  cdf = imhist.cumsum() #cumulative distribution function
77  cdf = cdf/cdf[-1]
78 
79  #use linear interpolation of cdf to find new pixel values
80  im2 = np.interp(im.flatten(),bins[:-1],cdf)
81 
82  plt.plot(bins[:-1],imhist)
83  plt.plot(bins[:-1],cdf)
84  #plt.show()
85  plt.clf()
86 
87  plt.plot(np.divide(im2[0:10000],im.flatten()[0:10000]))
88  #plt.show()
89  plt.clf()
90 
91  return im2.reshape(im.shape), cdf
92 
93 
94 def pseudocolor(val, minval, maxval):
95  # Scale val to be in the range [0, 1]
96  val = (val - minval) / (maxval - minval)
97  # Return RGBA tuple from nipy_spectral colormap
98  return cm.nipy_spectral(val)
99 
100 
101 # Start program execution
102 
103 args = sys.argv
104 
105 dtdb_dirpath = args[1]
106 dataset = args[2]
107 
108 colors = "123"
109 if len(args) > 3:
110  colors = args[3]
111  bc = np.int(args[3][0])-1
112  gc = np.int(args[3][1])-1
113  rc = np.int(args[3][2])-1
114 else:
115  bc = 1
116  gc = 2
117  rc = 3
118 
119 logfile = ""
120 if len(args) > 4:
121  logfile = args[4]
122  command = "date > " + logfile
123  result = os.system(command)
124 
125 dtdb_dircontents = os.listdir(dtdb_dirpath)
126 dtdb_dircontents.sort()
127 
128 for x in dtdb_dircontents:
129 
130  if x[0] == ".":
131  continue
132 
133  filepath = dtdb_dirpath + "/" + x
134 
135  if not os.path.isfile(filepath):
136  continue
137 
138  try:
139  if (dataset == "reflectance"):
140  refl = xr.load_dataset(filepath,group='/reflectance')['toa_reflectance']
141  bluec = refl[bc,:,:]
142  greenc = refl[gc,:,:]
143  redc = refl[rc,:,:]
144  elif (dataset == "aot"):
145  ocean = xr.load_dataset(filepath,group='/geophysical_data')['AOT_ocean']
146  land = xr.load_dataset(filepath,group='/geophysical_data')['AOT_land']
147  land = xr.load_dataset(filepath,group='/geolocation')['land_water']
148  bluec = ocean[bc,:,:]
149  greenc = ocean[gc,:,:]
150  redc = ocean[rc,:,:]
151  np.putmask(bluec, lw>0, land[0,:,:])
152  np.putmask(greenc, lw>0, land[1,:,:])
153  np.putmask(redc, lw>0, land[2,:,:])
154  elif (dataset == "cloudmask"):
155  cldmsk = xr.load_dataset(filepath,group='/meteorology')['cloud_mask']
156  cldtst = xr.load_dataset(filepath,group='/quality')['cloud_test']
157  bluec = np.zeros_like(cldmsk, dtype=np.int16)
158  greenc = np.zeros_like(cldmsk, dtype=np.int16)
159  redc = np.zeros_like(cldmsk, dtype=np.int16)
160  nomask = np.zeros_like(cldmsk, dtype=np.int16)
161  bluec = cldmsk[:,:]
162  np.putmask(greenc[:,:], cldtst[:,:]<0, -cldtst[:,:])
163  np.putmask(bluec[:,:], cldtst[:,:]<0, nomask[:,:])
164  np.putmask(redc[:,:], cldtst[:,:]>0, cldtst[:,:])
165  np.putmask(bluec[:,:], cldtst[:,:]>0, nomask[:,:])
166  colors = "rgb"
167  elif (dataset == "retrievals"):
168  ae = xr.load_dataset(filepath,group='/geophysical_data')['ae'][0,:,:].values
169  aot = xr.load_dataset(filepath,group='/geophysical_data')['aot_550'][:,:].values
170  fmf = xr.load_dataset(filepath,group='/geophysical_data')['fmf'][:,:].values
171  err = xr.load_dataset(filepath,group='/geophysical_data')['residual_error'][:,:].values
172  except Exception as e:
173  print (" Error loading dataset ... exiting")
174  continue
175 
176  if (dataset == "reflectance" or dataset == "aot" or dataset == "cloudmask"):
177  ro = np.clip(np.nan_to_num(redc),0,10.0)
178  go = np.clip(np.nan_to_num(greenc),0,10.0)
179  bo = np.clip(np.nan_to_num(bluec),0,10.0)
180 
181  if np.isnan(np.sum(ro)):
182  print (np.max(ro), np.max(go), np.max(bo))
183  continue
184 
185  y, u, v = rgb2yuv(ro, go, bo)
186  y += 0.001
187  y = y/np.max(y)
188  yh, Y_cdf = histeq(y)
189 # u *= 1
190 # v *= 1
191 # rh, gh, bh = yuv2rgb(yh, u, v)
192 
193  mino_all = min(np.min(ro), np.min(go), np.min(bo))
194  # ro -= mino_all
195  # go -= mino_all
196  # bo -= mino_all
197  maxo_all = max(np.max(ro), np.max(go), np.max(bo)) + 0.001
198  ro /= maxo_all
199  go /= maxo_all
200  bo /= maxo_all
201 
202  gamma = 1.0
203  ga = np.ones_like(ro)*gamma
204  rg = np.power(ro, ga)
205  gg = np.power(go, ga)
206  bg = np.power(bo, ga)
207  yg = np.power(yh, ga)
208 
209  ys = np.divide(yg,y)
210  rg = np.multiply(ys,ro)
211  gg = np.multiply(ys,go)
212  bg = np.multiply(ys,bo)
213 
214  scale = 255
215  im_red = im.fromarray(np.uint8(np.clip(rg,0,1.0)*scale))
216  im_green = im.fromarray(np.uint8(np.clip(gg,0,1.0)*scale))
217  im_blue = im.fromarray(np.uint8(np.clip(bg,0,1.0)*scale))
218 
219  im_rgb = im.merge("RGB", (im_red, im_green, im_blue))
220 
221  # im_rgb.show()
222  set = dataset.partition("/")
223  outfilename = x + "_" + set[2] + "_" + colors + ".png"
224 
225  image_path = dtdb_dirpath + "/IMAGES"
226  if not os.path.isdir(image_path):
227  os.mkdir(image_path)
228  image_path = image_path + "/" + set[0]
229  if not os.path.isdir(image_path):
230  os.mkdir(image_path)
231  image_path = image_path + set[1] + set[2] + "_" + colors
232  if not os.path.isdir(image_path):
233  os.mkdir(image_path)
234 
235  outpath = image_path + "/" + outfilename
236  im_rgb.save( outpath )
237 
238  elif (dataset == "retrievals"):
239  image_path = dtdb_dirpath + "/IMAGES"
240  if not os.path.isdir(image_path):
241  os.mkdir(image_path)
242  image_path = image_path + "/geophysical_data"
243 
244  outfilename = x + "_" + "ae" + ".png"
245  path = image_path + "_ae"
246  if not os.path.isdir(path):
247  os.mkdir(path)
248  outpath = path + "/" + outfilename
249  maxa = 2.0
250  mina = -1.0
251  np.clip(ae, mina, maxa, ae)
252  aep = (ae - mina)/(maxa - mina)
253  ml = mpl.cm.ScalarMappable(norm=None, cmap='nipy_spectral')
254  ml.set_array(ae)
255  aec = ml.to_rgba(aep, alpha=None)
256  plt.clf()
257  fig = plt.figure(facecolor='k')
258  fig.figimage(aec,resize=True)
259  fig.set_figwidth(fig.get_figwidth()*1.1)
260  ax = fig.add_axes([0.94, 0.1, 0.012, 0.8], frameon=True)
261  cbar = fig.colorbar(ml, ax=ax, cax=ax)
262  cbar.ax.tick_params(colors='c', labelsize=10, length=10, width=2, pad=5, right=True)
263  cbar.ax.yaxis.set_major_locator(plt.MaxNLocator(12))
264  fig.suptitle(x +' Angstrom Exponent', color='c', size=20)
265 # plt.show()
266  plt.savefig(outpath, facecolor=fig.get_facecolor())
267  plt.close(fig)
268 
269  outfilename = x + "_" + "aot_550" + ".png"
270  path = image_path + "_aot_550"
271  if not os.path.isdir(path):
272  os.mkdir(path)
273  outpath = path + "/" + outfilename
274  maxa = 2.0
275  mina = 0.0
276  if (np.max(aot) < maxa):
277  maxa = np.max(aot)
278  np.clip(aot, mina, maxa, aot)
279  aotp = (aot - mina)/(maxa - mina)
280  ml = mpl.cm.ScalarMappable(norm=None, cmap='nipy_spectral')
281  ml.set_array(aot)
282  aotc = ml.to_rgba(aotp, alpha=None, bytes=True, norm=False)
283  plt.clf()
284  fig = plt.figure(facecolor='k')
285  fig.figimage(aotc,resize=True)
286  fig.set_figwidth(fig.get_figwidth()*1.1)
287  ax = fig.add_axes([0.94, 0.1, 0.012, 0.8], frameon=True)
288  cbar = fig.colorbar(ml, ax=ax, cax=ax)
289  cbar.ax.tick_params(colors='c', labelsize=10, length=10, width=2, pad=5, right=True)
290  cbar.ax.yaxis.set_major_locator(plt.MaxNLocator(12))
291  fig.suptitle(x +' Aerosol Optical Depth at 550 nm', color='c', size=20)
292  plt.savefig(outpath, facecolor=fig.get_facecolor())
293  plt.close(fig)
294 
295  outfilename = x + "_" + "fmf" + ".png"
296  path = image_path + "_fmf"
297  if not os.path.isdir(path):
298  os.mkdir(path)
299  outpath = path + "/" + outfilename
300  maxa = 1.0
301  mina = 0.0
302  np.clip(fmf, mina, maxa, fmf)
303  fmfp = (fmf - mina)/(maxa - mina)
304  ml = mpl.cm.ScalarMappable(norm=None, cmap='nipy_spectral')
305  ml.set_array(fmf)
306  fmfc = ml.to_rgba(fmfp, alpha=None, bytes=True, norm=False)
307  plt.clf()
308  fig = plt.figure(facecolor='k')
309  fig.figimage(fmfc,resize=True)
310  fig.set_figwidth(fig.get_figwidth()*1.1)
311  ax = fig.add_axes([0.94, 0.1, 0.012, 0.8], frameon=True)
312  cbar = fig.colorbar(ml, ax=ax, cax=ax)
313  cbar.ax.tick_params(colors='c', labelsize=10, length=10, width=2, pad=5, right=True)
314  cbar.ax.yaxis.set_major_locator(plt.MaxNLocator(12))
315  fig.suptitle(x +' Fine Mode Fraction', color='c', size=20)
316  plt.savefig(outpath, facecolor=fig.get_facecolor())
317  plt.close(fig)
318 
319  outfilename = x + "_" + "error" + ".png"
320  path = image_path + "_error"
321  if not os.path.isdir(path):
322  os.mkdir(path)
323  outpath = path + "/" + outfilename
324  maxa = 10.0
325  mina = 0.0
326  np.clip(err, mina, maxa, err)
327  errp = (err - mina)/(maxa - mina)
328  ml = mpl.cm.ScalarMappable(norm=None, cmap='nipy_spectral')
329  ml.set_array(err)
330  errc = ml.to_rgba(errp, alpha=None, bytes=True, norm=False)
331  plt.clf()
332  fig = plt.figure(facecolor='k')
333  fig.figimage(errc,resize=True)
334  fig.set_figwidth(fig.get_figwidth()*1.1)
335  ax = fig.add_axes([0.94, 0.1, 0.012, 0.8], frameon=True)
336  cbar = fig.colorbar(ml, ax=ax, cax=ax)
337  cbar.ax.tick_params(colors='c', labelsize=10, length=10, width=2, pad=5, right=True)
338  cbar.ax.yaxis.set_major_locator(plt.MaxNLocator(20))
339  fig.suptitle(x +' Sum of Squares Error', color='c', size=20)
340  plt.savefig(outpath, facecolor=fig.get_facecolor())
341  plt.close(fig)
342 
343  print (x)
344 
def ycbcr2rgb(y, cb, cr)
Definition: color_dtdb.py:35
def rgb2ycbcr2(r, g, b)
Definition: color_dtdb.py:58
def yuv2rgb(y, u, v)
Definition: color_dtdb.py:51
def rgb2ycbcr(r, g, b)
Definition: color_dtdb.py:29
def histeq(im, nbr_bins=256)
Definition: color_dtdb.py:71
def rgb2yuv(r, g, b)
Definition: color_dtdb.py:45
def pseudocolor(val, minval, maxval)
Definition: color_dtdb.py:94
def ycbcr2rgb2(y, pb, pr)
Definition: color_dtdb.py:64