<jupyter_start><jupyter_text>A cloud effective radius (CER) and optical thickness (COT) retrieval framework using python; scikit-learn and TensorFlow: Application to the Moderate Resolution Imaging Spectroradiometer (MODIS) on board NASA's Terra & Aqua satellites.Author: <NAME> <EMAIL> <EMAIL><jupyter_code>import tensorflow as tf
import matplotlib.pyplot as plt
import numpy as np
import math
import matplotlib.patches as mpatches
import matplotlib.cm as cm
import pandas as pd
import matplotlib as mpl
import tensorflow_docs as tfdocs
import tensorflow_docs.plots
import tensorflow_docs.modeling
import os
import urllib.request
import tarfile
import warnings

from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF
from sklearn.gaussian_process.kernels import DotProduct, ConstantKernel as C
from sklearn import preprocessing
from pyhdf.SD import SD, SDC 
from matplotlib.pyplot import figure
from sklearn.preprocessing import PolynomialFeatures
from sklearn.linear_model import LinearRegression
from sklearn.pipeline import Pipeline
from sklearn.pipeline import make_pipeline
from sklearn.linear_model import Ridge
from pylab import figure, cm
from scipy import misc
from joblib import dump, load
from scipy.optimize import minimize, rosen, rosen_der    
from tensorflow import keras
from tensorflow.keras import layers

warnings.filterwarnings('ignore')
from IPython.core.display import HTML
HTML("""
<style>
.output_png {
    display: table-cell;
    text-align: center;
    vertical-align: middle;
}
</style>
""")<jupyter_output><empty_output><jupyter_text>Download ice and liquid reflectance libraries[Look Up Tables (LUTs)](https://atmosphere-imager.gsfc.nasa.gov/products/cloud/luts)<jupyter_code>%%time

if not os.path.exists('MODIS_C6_LUTS'):

    # Download LUT libraries
    
    url = 'https://atmosphere-imager.gsfc.nasa.gov/sites/default/files/ModAtmo/resources/modis_c6_luts.tar.gz'
    
    downloaded_filename = 'modis_c6_luts.tar.gz'
    
    urllib.request.urlretrieve(url, downloaded_filename)
    
    # Unzip .tar.gz
    # Ref: https://stackoverflow.com/questions/30887979/i-want-to-create-a-script-for-unzip-tar-gz-file-via-python
 
    fname = 'modis_c6_luts.tar.gz'

    if fname.endswith("tar.gz"):
   
        tar = tarfile.open(fname, "r:gz")
        tar.extractall()
        tar.close()
reflectance_ice_library = SD('./MODIS_C6_LUTS/modis_ice/reflectance_ice_library.hdf', SDC.READ)

print( reflectance_ice_library.info() )
datasets_dic = reflectance_ice_library.datasets()

for idx,sds in enumerate(datasets_dic.keys()):
    #print( idx,sds )
    sds_obj = reflectance_ice_library.select(sds)
    print( sds_obj.info() )
MultiScatBDReflectance = reflectance_ice_library.select('MultiScatBDReflectance')

MultiScatBDReflectance = MultiScatBDReflectance.get()

print(MultiScatBDReflectance.shape)

fig = figure(num=None, figsize=(12, 10), dpi=80, facecolor='w', edgecolor='k')

ax = fig.add_subplot(111)
    
plt.imshow(MultiScatBDReflectance[8,3,13,:,4,:], cmap='hot', interpolation='nearest',origin='lower')

plt.colorbar()

ax.set_aspect('auto')

plt.show()
OpticalThickness = reflectance_ice_library.select('OpticalThickness')
OpticalThickness = OpticalThickness.get()

print(OpticalThickness.shape)

print(OpticalThickness)
ParticleRadius = reflectance_ice_library.select('ParticleRadius')
ParticleRadius = ParticleRadius.get()

print(ParticleRadius.shape)

print(ParticleRadius)
y = MultiScatBDReflectance[8,3,13,:,4,:]

print(y.shape)
print( y.max() )
print( y.min() )
x1, x2 = np.meshgrid(np.arange(0,12, 1), np.arange(0,34, 1))

print(34*12)

x1 = x1.flatten()
x2 = x2.flatten()

x1 = ParticleRadius[x1]
x2 = OpticalThickness[x2]

y = y.flatten()
X = np.stack((x1,x2),axis=-1)
X = pd.DataFrame(X,columns=['tau','cer'])
X_stats = X.describe()
X_stats = X_stats.transpose()
X_stats
def data_scaling_function(x):
  return (x - X_stats['mean']) / X_stats['std']
X = data_scaling_function(X)

test_stats = X.describe()
test_stats = test_stats.transpose()
test_stats
y = pd.DataFrame(y,columns=['MultiScatBDReflectance'])
t = pd.concat((X,y),axis=1)
t = t.sample(frac = 1) 

t
X_train = t[['tau','cer']]
y_train = t[['MultiScatBDReflectance']]
def build_regression_model(intput_x):
    model = keras.Sequential([
    layers.Dense(10, activation='tanh', input_shape=[intput_x.shape[1]]),
    layers.Dense(10, activation='tanh'),
    layers.Dense(1)
    ])
    
    opt = tf.keras.optimizers.RMSprop(0.01)
    
    model.compile(loss='mse',
                optimizer='adam',
                metrics=['mse'])

    return model

reg_model = build_regression_model(X_train)
reg_model.fit(X_train, y_train, epochs=3000,
          validation_split = 0.2, verbose=0, callbacks=[tfdocs.modeling.EpochDots()])
X_test = np.stack((x1,x2),axis=-1)

X_test = pd.DataFrame(X_test,columns=['tau','cer'])

X_test = data_scaling_function(X_test)
yhat = reg_model.predict(X_test)
print(yhat.min())
print(yhat.max())
fig = figure(num=None, figsize=(12, 10), dpi=80, facecolor='w', edgecolor='k')

ax = fig.add_subplot(111)
    
#plt.imshow(yhat.reshape(34,12), cmap='hot', interpolation='nearest',origin='lower',vmin=0.0005173092,vmax=0.4885463)
plt.imshow(yhat.reshape(34,12), cmap='hot', interpolation='nearest',origin='lower',vmin=0.0005173092,vmax=0.4885463)

plt.colorbar()

ax.set_aspect('auto')

plt.show()
y_test = y.to_numpy()

y_test = y_test.reshape(-1,1)

diff = np.sqrt( (yhat - y_test) * (yhat - y_test) )

for i in range(20):
        print(yhat[i,0],y_test[i,0],diff[i,0])
fig = figure(num=None, figsize=(12, 10), dpi=80, facecolor='w', edgecolor='k')

ax = fig.add_subplot(111)
    
plt.imshow(diff.reshape(34,12), cmap='hot', interpolation='nearest',origin='lower')


plt.colorbar()

ax.set_aspect('auto')

plt.show()
x1, x2 = np.meshgrid(np.arange(0,60, 0.1), np.arange(0,160, 0.1))

x1 = x1.flatten()
x2 = x2.flatten()

X_test = np.stack((x1,x2),axis=-1)

X_test = pd.DataFrame(X_test,columns=['tau','cer'])

X_test = data_scaling_function(X_test)

yhat = reg_model.predict(X_test)

fig = figure(num=None, figsize=(12, 10), dpi=80, facecolor='w', edgecolor='k')

ax = fig.add_subplot(111)
    
plt.imshow(yhat.reshape(1600,600), cmap='hot', interpolation='nearest',origin='lower',vmin=0.0005173092,vmax=0.4885463)

plt.colorbar()

ax.set_aspect('auto')

plt.show()
y1 = MultiScatBDReflectance[8,3,13,:,0,:]
y2 = MultiScatBDReflectance[8,3,13,:,4,:]

y1 = y1.flatten()
y2 = y2.flatten()

Y = np.stack((y1,y2),axis=-1)

Y = pd.DataFrame(Y,columns=['band 1', 'band 6'])

Y
t = pd.concat((X,Y),axis=1)

t = t.sample(frac = 1) 

t
X_train = t[['tau','cer']]
y_train = t[['band 1', 'band 6']]
def build_multiple_regression_model(intput_x):
    model = keras.Sequential([
    layers.Dense(60, activation='tanh', input_shape=[intput_x.shape[1]]),
    layers.Dense(60, activation='tanh'),
    layers.Dense(60, activation='tanh'),
    layers.Dense(2)
    ])
    
    opt = tf.keras.optimizers.RMSprop(0.01)
    
    model.compile(loss='mse',
                optimizer='adam',
                metrics=['mse'])

    return model

mult_reg_model = build_multiple_regression_model(X_train)
mult_reg_model.fit(X_train, y_train, epochs=20000,
          validation_split = 0.2, verbose=0, callbacks=[tfdocs.modeling.EpochDots()])
x1, x2 = np.meshgrid(np.arange(0,60, 1), np.arange(0,160, 1))

x1 = x1.flatten()
x2 = x2.flatten()

X_test = np.stack((x1,x2),axis=-1)

X_test = pd.DataFrame(X_test,columns=['tau','cer'])

X_test = data_scaling_function(X_test)

yhat = mult_reg_model.predict(X_test)

yhat

fig = figure(num=None, figsize=(12, 10), dpi=80, facecolor='w', edgecolor='k')

ax = fig.add_subplot(111)
    
plt.imshow(yhat[:,0].reshape(160,60), cmap='hot', interpolation='nearest',origin='lower')

plt.colorbar()

ax.set_aspect('auto')

plt.show()
fig = figure(num=None, figsize=(12, 10), dpi=80, facecolor='w', edgecolor='k')

ax = fig.add_subplot(111)
    
plt.imshow(yhat[:,1].reshape(160,60), cmap='hot', interpolation='nearest',origin='lower')

plt.colorbar()

ax.set_aspect('auto')

plt.show()
file = SD('./MODIS_C6_LUTS/modis_ice/reflectance_ice_library.hdf', SDC.READ)

x_band_idx = 0
y_band_idx = 4


ice_lut_shape = MultiScatBDReflectance.shape

figure(num=None, figsize=(12, 10), dpi=80, facecolor='w', edgecolor='k')

x = [MultiScatBDReflectance[8,3,13,:,x_band_idx,index_re] for index_re in np.arange(ice_lut_shape[5])]
y = [MultiScatBDReflectance[8,3,13,:,y_band_idx,index_re] for index_re in np.arange(ice_lut_shape[5])]

plt.plot(x,y, 'steelblue',linewidth=0.5)

x = [MultiScatBDReflectance[8,3,13,index_tau,x_band_idx,:] for index_tau in np.arange(ice_lut_shape[3])]
y = [MultiScatBDReflectance[8,3,13,index_tau,y_band_idx,:] for index_tau in np.arange(ice_lut_shape[3])]

plt.plot(x,y, 'coral',linewidth=0.5)

pop_a = mpatches.Patch(color='coral', label=r'Cloud Effective Radius ($\mu m$)')
pop_b = mpatches.Patch(color='steelblue', label='Cloud Optical Thickness')

plt.legend(handles=[pop_a,pop_b],fontsize=8)


plt.show()
plt.close()


x1, x2 = np.meshgrid(np.arange(0,12, 1), np.arange(0,34, 1))

print(34*12)

x1 = x1.flatten()
x2 = x2.flatten()

x1 = ParticleRadius[x1]
x2 = OpticalThickness[x2]

X_test = np.stack((x1,x2),axis=-1)

X_test = pd.DataFrame(X_test,columns=['tau','cer'])

X_test = data_scaling_function(X_test)

yhat = mult_reg_model.predict(X_test)

yhat.shape

yhat = yhat.reshape(34,12,2)
figure(num=None, figsize=(12, 10), dpi=80, facecolor='w', edgecolor='k')



x = [yhat[:,index_re,0] for index_re in np.arange(ice_lut_shape[5])]
y = [yhat[:,index_re,1] for index_re in np.arange(ice_lut_shape[5])]

plt.plot(x,y, 'steelblue',linewidth=0.5)


x = [yhat[index_tau,:,0] for index_tau in np.arange(ice_lut_shape[3])]
y = [yhat[index_tau,:,1] for index_tau in np.arange(ice_lut_shape[3])]

plt.plot(x,y, 'coral',linewidth=0.5)

pop_a = mpatches.Patch(color='steelblue', label=r'Cloud Effective Radius ($\mu m$)')
pop_b = mpatches.Patch(color='coral', label='Cloud Optical Thickness')

plt.legend(handles=[pop_a,pop_b],fontsize=8)


plt.show()
plt.close()


figure(num=None, figsize=(12, 10), dpi=80, facecolor='w', edgecolor='k')

x = [MultiScatBDReflectance[8,3,13,:,x_band_idx,index_re] for index_re in np.arange(ice_lut_shape[5])]
y = [MultiScatBDReflectance[8,3,13,:,y_band_idx,index_re] for index_re in np.arange(ice_lut_shape[5])]

plt.plot(x,y, 'steelblue',linewidth=0.5)


x = [yhat[:,index_re,0] for index_re in np.arange(ice_lut_shape[5])]
y = [yhat[:,index_re,1] for index_re in np.arange(ice_lut_shape[5])]

plt.plot(x,y, 'steelblue',linewidth=0.5)

x = [MultiScatBDReflectance[8,3,13,index_tau,x_band_idx,:] for index_tau in np.arange(ice_lut_shape[3])]
y = [MultiScatBDReflectance[8,3,13,index_tau,y_band_idx,:] for index_tau in np.arange(ice_lut_shape[3])]

plt.plot(x,y, 'coral',linewidth=0.5)

x = [yhat[index_tau,:,0] for index_tau in np.arange(ice_lut_shape[3])]
y = [yhat[index_tau,:,1] for index_tau in np.arange(ice_lut_shape[3])]

plt.plot(x,y, 'coral',linewidth=0.5)

pop_a = mpatches.Patch(color='coral', label=r'Cloud Effective Radius ($\mu m$)')
pop_b = mpatches.Patch(color='steelblue', label='Cloud Optical Thickness')

plt.legend(handles=[pop_a,pop_b],fontsize=8)

plt.show()
plt.close()
from graphviz import Digraph
from graphviz import Source

temp = '''
digraph G {


     graph[ fontname = "Helvetica-Oblique",
            fontsize = 12,
            label = "",
            size = "7.75,10.25" ];

    rankdir = LR;
    splines=false;
    edge[style=invis];
    ranksep= 1.4;
    {
    node [shape=circle, color=chartreuse, style=filled, fillcolor=chartreuse];
    x1 [label=<CER>];
    x2 [label=<TAU>]; 
}
{
    node [shape=circle, color=dodgerblue, style=filled, fillcolor=dodgerblue];
    a12 [label=<a<sub>1</sub><sup>(2)</sup>>];
    a22 [label=<a<sub>2</sub><sup>(2)</sup>>];
    a32 [label=<a<sub>3</sub><sup>(2)</sup>>];
    a42 [label=<a<sub>4</sub><sup>(2)</sup>>];
    a52 [label=<a<sub>5</sub><sup>(2)</sup>>];
    a13 [label=<a<sub>1</sub><sup>(3)</sup>>];
    a23 [label=<a<sub>2</sub><sup>(3)</sup>>];
    a33 [label=<a<sub>3</sub><sup>(3)</sup>>];
    a43 [label=<a<sub>4</sub><sup>(3)</sup>>];
    a53 [label=<a<sub>5</sub><sup>(3)</sup>>];
}
{
    node [shape=circle, color=coral1, style=filled, fillcolor=coral1];
    O1 [label=<Band 1>];
    O2 [label=<Band 6>]; 
    O3 [label=<Band 7>]; 
}
    {
        rank=same;
        x1->x2;
    }
    {
        rank=same;
        a12->a22->a32->a42->a52;
    }
    {
        rank=same;
        a13->a23->a33->a43->a53;
    }
    {
        rank=same;
        O1->O2->O3;
    }
    l0 [shape=plaintext, label="layer 1 (input layer)"];
    l0->x1;
    {rank=same; l0;x1};
    l1 [shape=plaintext, label="layer 2 (hidden layer)"];
    l1->a12;
    {rank=same; l1;a12};
    l2 [shape=plaintext, label="layer 3 (hidden layer)"];
    l2->a13;
    {rank=same; l2;a13};
    l3 [shape=plaintext, label="layer 4 (output layer)"];
    l3->O1;
    {rank=same; l3;O1};
    edge[style=solid, tailport=e, headport=w];
    {x1; x2} -> {a12;a22;a32;a42;a52};
    {a12;a22;a32;a42;a52} -> {a13;a23;a33;a43;a53};
    {a13;a23;a33;a43;a53} -> {O1,O2,O3};
}'''

s = Source(temp)

s
from graphviz import Digraph
from graphviz import Source

temp = '''
digraph G {


     graph[ fontname = "Helvetica-Oblique",
            fontsize = 12,
            label = "",
            size = "7.75,10.25" ];

    rankdir = LR;
    splines=false;
    edge[style=invis];
    ranksep= 1.4;
    {
    node [shape=circle, color=chartreuse, style=filled, fillcolor=chartreuse];
    x1 [label=<x1>];
    x2 [label=<x2>]; 
}
{
    node [shape=circle, color=dodgerblue, style=filled, fillcolor=dodgerblue];
    a12 [label=<a<sub>1</sub><sup>(2)</sup>>];
    a22 [label=<a<sub>2</sub><sup>(2)</sup>>];
    a32 [label=<a<sub>3</sub><sup>(2)</sup>>];
    a42 [label=<a<sub>4</sub><sup>(2)</sup>>];
    a52 [label=<a<sub>5</sub><sup>(2)</sup>>];
    a13 [label=<a<sub>1</sub><sup>(3)</sup>>];
    a23 [label=<a<sub>2</sub><sup>(3)</sup>>];
    a33 [label=<a<sub>3</sub><sup>(3)</sup>>];
    a43 [label=<a<sub>4</sub><sup>(3)</sup>>];
    a53 [label=<a<sub>5</sub><sup>(3)</sup>>];
}
{
    node [shape=circle, color=coral1, style=filled, fillcolor=coral1];
    O1 [label=<y1>];
    O2 [label=<y2>]; 
    O3 [label=<y3>]; 
}
    {
        rank=same;
        x1->x2;
    }
    {
        rank=same;
        a12->a22->a32->a42->a52;
    }
    {
        rank=same;
        a13->a23->a33->a43->a53;
    }
    {
        rank=same;
        O1->O2->O3;
    }
    l0 [shape=plaintext, label="layer 1 (input layer)"];
    l0->x1;
    {rank=same; l0;x1};
    l1 [shape=plaintext, label="layer 2 (hidden layer)"];
    l1->a12;
    {rank=same; l1;a12};
    l2 [shape=plaintext, label="layer 3 (hidden layer)"];
    l2->a13;
    {rank=same; l2;a13};
    l3 [shape=plaintext, label="layer 4 (output layer)"];
    l3->O1;
    {rank=same; l3;O1};
    edge[style=solid, tailport=e, headport=w];
    {x1; x2} -> {a12;a22;a32;a42;a52};
    {a12;a22;a32;a42;a52} -> {a13;a23;a33;a43;a53};
    {a13;a23;a33;a43;a53} -> {O1,O2,O3};
}'''



s = Source(temp)

s<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.18091714021408747, 'fraction_numerical': 0.04929298268798731, 'mean_word_length': 4.730783794017418, 'pattern_counts': {'":': 0, '<': 124, '<?xml version=': 0, '>': 160, 'https://': 3, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '2', 'dir': 'jupyter-structured-clean-dedup', 'id': '504279', 'license': 'MIT', 'n_tokens_mistral': 6720, 'n_tokens_neox': 6321, 'n_words': 1000, 'path': 'docs/train_deep_learning_model_from_modis_myd06_lut.ipynb', 'repo_name': 'benjamin-hg-marchant/deep-inversion', 'size': '15185'}