Earthdata Forum

Code: Select all

import os
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
from pyhdf import HDF, SD, VS
from pyhdf.SD import SD, SDC

FILE_NAME = 'CAL_LID_L2_05kmAPro-Standard-V4-20.2015-05-14T23-44-05ZD.hdf'
DATAFIELD_NAME = 'Extinction_Coefficient_532'
sd = SD(FILE_NAME, SDC.READ)
hdf = HDF.HDF(FILE_NAME)
vs = hdf.vstart()
hdf = SD(FILE_NAME, SDC.READ)

#get altitude
xid = vs.find('metadata')
altid = vs.attach(xid)
altid.setfields('Lidar_Data_Altitudes')
nrecs, _, _, _, _ = altid.inquire()
altitude = altid.read(nRec=nrecs)
altid.detach()
alti = np.array(altitude[0][0])
alt = np.flip(alti)

# Read geolocation datasets.
# They are all same for every file so we read only once.
lat = hdf.select("Latitude")
latitude = np.squeeze(lat[:,0])

#subset for latitude (10-20N)
idx_lat = np.where((latitude >= 10) & (latitude <= 20))
idx_subset_lat = idx_lat[0]
subset_lat = np.array([latitude[i] for i in idx_subset_lat])

# Read dataset.
data2D = hdf.select(DATAFIELD_NAME)
data = data2D[:,:].astype(np.float64)
# Retrieve the attributes.
attrs = data2D.attributes(full=1)
# print(attrs)

fva=attrs["fillvalue"]
fill_value = fva[0]
ua=attrs["units"]
units = ua[0]

# Replace the missing values with NaN.        
data[data == float(fill_value)] = np.nan

#subset data to accomodate the subsetting of latitude
data_arrays = []
for x in range(len(data.T)):
   subset = np.take(data[:,x], idx_subset_lat)
   data_arrays.append(subset)
   if x >= len(data.T):
       break
subset_data = np.stack(data_arrays)


# Contour the data on a grid of latitude vs. altitude
lats, alts = np.meshgrid(subset_lat, alt)

#plot data
plt.figure(figsize=(7.20,3.60))
plt.contourf(lats, alts, subset_data, cmap=plt.get_cmap('jet'))
cb = plt.colorbar()
cb.set_label(units)
long_name = DATAFIELD_NAME 
plt.title("{0}\n{1}".format(FILE_NAME, long_name), fontsize=8)
plt.xlabel('Latitude (degrees north)')
plt.ylabel('Altitude (km)')
fig = plt.gcf()
pngfile = "{0}.v.py.png".format(FILE_NAME)
fig.savefig(pngfile, dpi=200)