How to plot MODIS C6 MYD06 Cloud Phase Optical Properties using python 3 ?

Examples of how to plot MODIS C6 MYD06 Cloud Phase Optical Properties using python 3:

Plot MODIS C6 MYD06 Cloud Phase Optical Properties

from pyhdf.SD import SD, SDC 
from matplotlib.pyplot import figure

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

modis_file = SD('../cer__tau_retrieval_module/inputs/MYD06_L2.A2008008.1420.006.2013342120951.hdf', SDC.READ)

data_selected_id = modis_file.select('Cloud_Phase_Optical_Properties')

data = data_selected_id.get()

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

cmap =  [(0.0,0.0,0.0)] + [(32./255., 107./255., 164./255.)] + \
       [(187./255., 217./255., 238./255.)] + [(255./255, 51./255, 51./255.)] 
cmap = mpl.colors.ListedColormap(cmap)

bounds = [-0.5, 1.5, 2.5, 3.5, 4.5]

norm = mpl.colors.BoundaryNorm(bounds, cmap.N)

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

cbar_bounds = [1.5, 2.5, 3.5, 4.5]
cbar_ticks = [2.0,3.0,4.0]  
cbar_labels = ['Liquid','Ice','Und.']

cbar = plt.colorbar(img, cmap=cmap, norm=norm, boundaries=cbar_bounds, ticks=cbar_ticks)
cbar.ax.set_yticklabels(cbar_labels, fontsize=10)

plt.title('Cloud Thermodynamic Phase 1km \n MYD06 C6 (2008-01-15; 14h35)', fontsize=10)

plt.savefig("modis_myd06_cloud_phase_optical_properties_python3_01.png", bbox_inches='tight')

plt.show()

Plot MODIS C6 MYD06 Cloud Phase Optical Properties with orthographic projection

To plot MODIS C6 MYD06 Cloud Phase Optical Properties with orthographic projection, a solution is to use basemap. Note: To fix the error message KeyError: 'PROJ_LIB':

conda install -c conda-forge proj4
conda install basemap

and add in the script

import os

os.environ['PROJ_LIB'] = '{{ path to anaconda }}/anaconda3/envs/worklab/share/proj'

replace worklab by the name of your environment (if you don't work in a specific environment: {{ path to anaconda }}/anaconda3/share/proj)

Code python to lot MODIS C6 MYD06 Cloud Phase Optical Properties:

#!/usr/bin/env python

from mpl_toolkits.basemap import Basemap, cm 
from pyhdf.SD import SD, SDC 
from scipy import interpolate
from scipy.interpolate import griddata

from matplotlib import colors
from pylab import figure, cm

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

# fix nearest interpolation create mask points too far
THRESHOLD = 5000.0

from scipy.interpolate.interpnd import _ndim_coords_from_arrays
from scipy.spatial import cKDTree

#----------------------------------------------------------------------------------------#
# Inputs

myd03_file_name = '../cer__tau_retrieval_module/inputs/MYD03.A2008008.1420.006.2012066144951.hdf'
myd06_file_name = '../cer__tau_retrieval_module/inputs/MYD06_L2.A2008008.1420.006.2013342120951.hdf'

#----------------------------------------------------------------------------------------#
# Read HDF Files

myd03 = SD(myd03_file_name, SDC.READ)
myd06 = SD(myd06_file_name, SDC.READ)

myd03_Latitude = myd03.select('Latitude')
myd03_Longitude = myd03.select('Longitude')
data_selected_id = myd06.select('Cloud_Phase_Optical_Properties')

myd03_Latitude_data = myd03_Latitude.get()
myd03_Longitude_data = myd03_Longitude.get()

#----------------------------------------------------------------------------------------#
# Find lat_0 and lon_0 of the MODIS granule

#print myd03_Latitude_data.shape
#print myd03_Longitude_data.shape

ArrayShape = myd03_Latitude_data.shape

#print ArrayShape[0], ArrayShape[1]

as0 = ArrayShape[0] - 1
as1 = ArrayShape[1] - 1

latmin = myd03_Latitude_data[0,0]
latmax = myd03_Latitude_data[as0,as1]

lat_0 = latmin + (latmax - latmin) / 2.

#print "lat_0:", lat_0

#----------------------------------------------------------------------------------------#

tmp_01 = myd03_Longitude_data[0,0]
tmp_02 = myd03_Longitude_data[as0,as1]

lonmin = min(myd03_Longitude_data[0,0],myd03_Longitude_data[as0,as1])
lonmax = max(myd03_Longitude_data[0,0],myd03_Longitude_data[as0,as1])

lon_0 = lonmin + (lonmax - lonmin) / 2.

if lon_0 > 180:
    lon_0 = - ( 360 - lon_0 )

#print "lon_0:", lon_0

myd03_Longitude_data[0,0] = tmp_01
myd03_Longitude_data[as0,as1] = tmp_02

#print latmin, latmax
#print lonmin, lonmax
#print "lat_0 and lon_0 done"

#----------------------------------------------------------------------------------------#
# Orthographic Map

fig = plt.figure(figsize=(6, 6), dpi=100)

ax = fig.add_subplot(111)
#ax.patch.set_facecolor('black')
ax.patch.set_facecolor((0.75,0.75,0.75))

m1 = Basemap(projection='ortho',lon_0=lon_0,lat_0=lat_0,resolution=None)

xpt0, ypt0 = m1(lon_0,lat_0)

#print 'xpt0, ypt0',xpt0, ypt0

xpt1, ypt1 = m1(myd03_Longitude_data[0,0],myd03_Latitude_data[0,0]) 
xpt2, ypt2 = m1(myd03_Longitude_data[0,as1],myd03_Latitude_data[0,as1]) 
xpt3, ypt3 = m1(myd03_Longitude_data[as0,as1],myd03_Latitude_data[as0,as1])
xpt4, ypt4 = m1(myd03_Longitude_data[as0,0],myd03_Latitude_data[as0,0])

llx = min(xpt1,xpt2,xpt3,xpt4) - xpt0  # lower left
lly = min(ypt1,ypt2,ypt3,ypt4) - ypt0

urx = max(xpt1,xpt2,xpt3,xpt4) - xpt0  # upper right
ury = max(ypt1,ypt2,ypt3,ypt4) - ypt0

m = Basemap(projection='ortho',lon_0=lon_0,lat_0=lat_0,resolution='l',\
    llcrnrx=llx,llcrnry=lly,urcrnrx=urx,urcrnry=ury)

#print "Orthographic map done"

#----------------------------------------------------------------------------------------#
# Plot MODIS data: Using Scipy Griddata

long = myd03_Longitude.get()
lat = myd03_Latitude.get()

x_igrid, y_igrid = m(long,lat)

x_igrid = x_igrid - xpt0
y_igrid = y_igrid - ypt0

z_igrid = data_selected_id.get()

attributes = data_selected_id.attributes()

z_igrid = z_igrid * attributes['scale_factor']

z_igrid_shape = z_igrid.shape

x1_igrid = x_igrid.ravel()
y1_igrid = y_igrid.ravel()
z1_igrid = z_igrid.ravel()

#print "z1_igrid.min(),z1_igrid.max()", z1_igrid.min(),z1_igrid.max()

xy1_igrid = np.vstack((x1_igrid, y1_igrid)).T

xi, yi = np.mgrid[llx:urx:1000j, lly:ury:1000j]

z = griddata(xy1_igrid, z1_igrid, (xi, yi), method='nearest')

#----------------------------------------------------------------------------------------#
# Construct kd-tree, functionality copied from scipy.interpolate

tree = cKDTree(xy1_igrid)
arr_x = _ndim_coords_from_arrays((xi, yi))
dists, indexes = tree.query(arr_x)

z = z[:]
z[dists > THRESHOLD] = np.nan

cmap =  [(0.0,0.0,0.0)] + [(32./255., 107./255., 164./255.)] + \
       [(187./255., 217./255., 238./255.)] + [(255./255, 51./255, 51./255.)] 
cmap = colors.ListedColormap(cmap)

bounds = [-0.5, 1.5, 2.5, 3.5, 4.5]
norm = mpl.colors.BoundaryNorm(bounds, cmap.N)

img = m.imshow(z.T, origin='lower', cmap=cmap, norm=norm )

cbar_bounds = [1.5, 2.5, 3.5, 4.5]
cbar_ticks = [2.0,3.0,4.0]  
cbar_labels = ['Liquid','Ice','Und.']

cbar = m.colorbar(img, cmap=cmap, norm=norm, boundaries=cbar_bounds, ticks=cbar_ticks)
cbar.ax.set_yticklabels(cbar_labels, fontsize=10)

m.drawcoastlines()

#m.fillcontinents(color='coral',lake_color='white')
#m.drawmapboundary(fill_color='white')

m.drawparallels(np.arange(-90.,120.,5.), color='k', labels=[True,False,False,False])
m.drawmeridians(np.arange(0.,420.,10.), color='k', labels=[False,False,False,True])

ax.set_xlabel("", fontsize=7)
ax.set_ylabel("", fontsize=7)

plt.title('Cloud Phase Optical Properties \n MYD06 C6', fontsize=10)

plt.savefig("modis_myd06_cloud_phase_optical_properties_python3.png", bbox_inches='tight')

plt.show()
plt.close()

References

• How can I avoid PROJ_LIB error in importing basemap?
• KeyError: 'PROJ_LIB' while installing basemap using conda install
• Basemap import error in PyCharm — KeyError: 'PROJ_LIB'