Ben October 08, 2019

Ongoing research ...

Goal: To develop an algorithm (in python) based on Gaussian processes to restore the Aqua MODIS Band 6 missing pixels:

### Get and plot the data

#!/usr/bin/env python

from pyhdf.SD import SD, SDC 
from pylab import figure, cm

import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
import matplotlib.cm as cm
import pprint

file_name = 'MYD021KM.A2007219.2010.006.2012077031102.hdf'

file = SD(file_name, SDC.READ)

datasets = file.datasets()

#for idx,sds in enumerate(datasets.keys()):
#    print( idx,sds )

data_selected = file.select('EV_500_Aggr1km_RefSB')

data = data_selected.get()

#print( data_selected.info() )

data_attributes = data_selected.attributes()

#pprint.pprint(data_attributes )

_FillValue = data_attributes['_FillValue']
_FillValue = 65528 # warning wrong _FillValue stored in attributes

reflectance_scales = data_attributes['reflectance_scales']

#----------------------------------------------------------------------------------------#
# Plot Modis band 06

data_band_6 = data[3,:,:]

data_band_6[ data_band_6 == _FillValue ] = 0.0

data_band_6 = data_band_6 * reflectance_scales[3]

data_band_6_vmin = np.min(data_band_6)
data_band_6_vmax = np.max(data_band_6)

data_shape = data_band_6.shape

print(data_shape)

cmap = [(0.0,0.0,0.0)] + [(cm.jet(i)) for i in range(1,256)] 
cmap = mpl.colors.ListedColormap(cmap)

img = plt.imshow(np.fliplr(data_band_6), cmap=cmap, 
                 interpolation='none', origin='lower')

plt.xticks([0,250,500,750,1000,1250], 
           ['0','250','500', '750', '1000','1250'], rotation=0, fontsize=8 )
plt.yticks([0,250,500,750,1000,1250,1500,1750,2000], 
           ['0','250','500', '750', '1000','1250','1500','1750','2000'], 
           rotation=0, fontsize=8 )

cbar = plt.colorbar()
cbar.ax.tick_params(labelsize=8)

plt.title('AQUA MODIS L1 Band 06', fontsize=8)

plt.savefig('modis_band6_missing_data_restoration_01.png', bbox_inches='tight')

#plt.show()

plt.close()


### Using only the spatial information

nb_missing_rows = np.where(data_band_6[:,0]==0.0)[0].shape[0]

band_06_missing_data_idx = np.where(data_band_6==0.0)

#print('data_shape --> ', data_shape)

#print('nb_missing_rows --> ', nb_missing_rows)

#print('number of missing data --> ', nb_missing_rows * data_shape[1])

#print('Fraction of missing data --> ', 100.0 * nb_missing_rows * data_shape[1] / ( data_shape[0] * data_shape[1] ))

#for i in range(100):
#   print(i,data_band_6[i,0])

nbx = 20

x = np.arange(nbx)

plt.scatter(x,data_band_6[0:nbx,0])

plt.savefig('modis_band6_missing_data_restoration_03.png', bbox_inches='tight')

#plt.show()

plt.close()

i_start = 0 #0
j_start = 0 #0

cmap = [(0.0,0.0,0.0)] + [(cm.jet(i)) for i in range(1,256)] 
cmap = mpl.colors.ListedColormap(cmap)

#img = plt.imshow(np.fliplr(data_band_6[0:12,0:12]), cmap=cmap, 
#                 interpolation='none', origin='lower')

img = plt.imshow(data_band_6[i_start:i_start+12,j_start:j_start+12], cmap=cmap, 
                 interpolation='none', origin='lower')

cbar = plt.colorbar()
cbar.ax.tick_params(labelsize=8)

plt.title('Band 06', fontsize=10)

plt.savefig('modis_band6_missing_data_restoration_05.png', bbox_inches='tight')

#plt.show()

plt.close()