import numpy as np
import netCDF4
import matplotlib.pyplot as plt
import scipy.ndimage as ndimage

# read data
dataf = "../data/uvi_20181126_150113_283_l3bx_v10.nc"

nc = netCDF4.Dataset(dataf, 'r')
rad = nc.variables["radiance"][0, :, :]   # radiance (unit: W/m2/str/m)
inang = nc.variables["inangle"][0, :, :]   # incidence angle (unit: radian)
x = nc.variables["axis1"][:]   # pixel coordinate -X
y = nc.variables["axis2"][:]   # pixel coordinate -Y

# plot the original image
plt.figure(1, (10, 10))
plt.title("original image")
plt.pcolormesh(x, y, rad, shading='auto', cmap="gray")
plt.show()

# model image created from cos(incidence angle)

# ---- [practice 1] ---- # (Use np.cos(), np.pi, inang)
model = np.cos(np.pi*inang/180.)

model = np.where(model > 0., model, 0.)  # Negative value is replaced by 0

# plot the model image
plt.figure(1, (10, 10))
plt.title("model image")
plt.pcolormesh(x, y, model, shading='auto', cmap="gray")
plt.show()

C:\ProgramData\Anaconda3\lib\site-packages\numpy\ma\core.py:1015: RuntimeWarning: overflow encountered in multiply
  result = self.f(da, db, *args, **kwargs)

# photometric correction (dividing the original image by the model image)

# ---- [practice 2] ---- # (Use rad, model)
corrected_rad = rad/model

# plot the image after photometric correction
plt.figure(1, (10, 10))
plt.title("corrected image")
plt.pcolormesh(x, y, corrected_rad, shading='auto', cmap="gray", vmin=0, vmax=5e7)
plt.show()

# high-pass filtering

# ---- [practice 3] ---- #
sigma = 8      # e-folding half width in pixels

smoothed_rad = ndimage.gaussian_filter(corrected_rad, [sigma, sigma])

# plot the smoothed image
plt.figure(1, (10, 10))
plt.title("smoothed image")
plt.pcolormesh(x, y, smoothed_rad, shading='auto', cmap="gray", vmin=0, vmax=5e7)
plt.show()

# high-pass filtering
highpass_rad = corrected_rad - smoothed_rad

# plot the highpass-filtered image (gray scale)
plt.figure(1, (10, 10))
plt.title("highpass-filtered image")
plt.pcolormesh(x, y, highpass_rad, shading='auto', cmap="gray", vmin=-6e5, vmax=6e5)
plt.show()

# plot the highpass-filtered image (color scale)
plt.figure(1, (10, 10))
plt.title("highpass-filtered image")
plt.pcolormesh(x, y, highpass_rad, shading='auto', cmap="jet", vmin=-6e5, vmax=6e5)
plt.show()