Read radiance data from NetCDF files

import numpy as np
import netCDF4
import matplotlib.pyplot as plt
import scipy.ndimage as ndimage
import glob
from PIL import Image

# open data files
dataf = sorted(glob.glob("../data/" + "uvi_2018111[7-9]*l3b_*.nc"))   # all images during 2018-11-17 -- 19
nf = len(dataf)   # number of files
print("number of images = ", nf)

# extract information from the first file
nc = netCDF4.Dataset(dataf[0], 'r')
rad = nc.variables["radiance"][0, :, :]   # radiance
nx = rad.shape[1]   # number of longitudinal grids
ny = rad.shape[0]   # number of latitudinal grids
lon = nc.variables["longitude"][:]   # longitude
lat = nc.variables["latitude"][:]   # latitude
nc.close()

number of images =  35

Average photometrically-corrected images and plot the averaged image. Photometric correction has been done by dividing original image by cosine of incidence angle.

# averaging photometrically-corrected images

# prepare a list for the averaged radiance
averaged_rad = np.zeros([ny, nx])

# prepare lists for each processed image
model = np.zeros([ny, nx])
corrected_rad = np.zeros([ny, nx])
smoothed_rad = np.zeros([ny, nx])
highpass_rad = np.zeros([ny, nx])

# ---- [practice 5] ---- #
n = 35       # not greater than 35

for i in range(0, n):   # read all images
    print(i, end=" ")

    # read data
    nc = netCDF4.Dataset(dataf[i], 'r')
    rad = nc.variables["radiance"][0, :, :]   # radiance
    inang = nc.variables["inangle"][0, :, :]   # incidence angle
    nc.close()

    # photometric correction
    model = np.cos(np.pi/180.*inang)
    model = np.where(model > 0., model, 0.)
    corrected_rad = rad/model

    # avaraging
    averaged_rad = averaged_rad + corrected_rad / float(n)

# plot the averaged image
plt.figure(1, (12, 6))
plt.title("averaged image")
plt.pcolormesh(lon, lat, averaged_rad, shading='auto', cmap="gray", vmin=0, vmax=1e8)
plt.show()

0 

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

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 

Apply high-pass filter to the averaged image and plot it to find feature fixed to the latitude-longitude points, i.e., the geography of Venus.

# high-pass filtering of the averaged image

# ---- [practice 6] ---- #
sigma = 20     # e-folding half width in pixels for gaussian smoothing

smoothed_rad = ndimage.gaussian_filter(averaged_rad, [sigma, sigma])
highpass_rad = averaged_rad - smoothed_rad

# plot high-passed image
plt.figure(2, (12, 6))
plt.title("averaged & high-pass filtered image")
plt.pcolormesh(lon, lat, highpass_rad, shading='auto', cmap="gray", vmin=-1e6, vmax=1e6)
plt.show()