CAL_LID_L2_333mCLay-ValStage1-V3-01.2006-10-30T21-41-43ZD.hdf
CAL_LID_L2_333mCLay-ValStage1-V3-01.2006-10-30T21-41-43ZD.h501_MODIS_L1B.ipynb.01_MODIS_L1B.ipynb and Section 3 for more information.
In [0]:
__author__ = 'ATSC-301 UBC'
CALIPSO is a joint NASA (USA) and CNES (France) environmental satellite, built in the Cannes Mandelieu Space Center, which was launched atop a Delta II rocket on April 28, 2006. Its name stands for Cloud-Aerosol Lidar and Infrared Pathfinder Satellite Observations.
Passive and active remote sensing Instruments on board the CALIPSO satellite monitor aerosols and clouds 24 hours a day. CALIPSO is part of the "A Train", flying in formation with several other satellites (Aqua, Aura and CloudSat).
You can findout more about CALIPSO's news and applications via **NASA's Satellite Mission Page**.
CALIPSO datasets are available in NASA's **Atmospheric Science Data Center**.
In order to get the datasets, we do the following steps:
We choose Cloud Profiles and Layers Product in CALIPSO L2 Lidar product family, it contains the following parameters:
There are three different choice of cloud layer resolution, we select **1/3km version**.
After you downloaded the CALIPSO files via FTP client. You need to convert HDF4 files to HDF5.
We have discussed this part in 01_MODIS_L1B.ipynb. But here we have a large number of files, a recommended way is to use h4toh5convert in batch sctipt.
_share/. copy *.bat to the direcctory of your HDF4 files (here say _data/CloudSat_L2), open the Terminal, cd to that path and type:dos
DOS_batch_h4toh5.batbash BASH_batch_h4toh5.sh
Bash works for Mac <--------------- !!! Except for what we will do in the following sections, **ccplot** is another choice for you:
ccplot is an open source command-line program for plotting profile, layer and earth view data sets from CloudSat, CALIPSO and Aqua MODIS products.
ccplot works for Windows, Linux and Mac. Visualising Data from CloudSat and CALIPSO Satellites - Peter Kuma.pdf in \_docs is its documentation. It also has a **Github Page**.
There is a CALIPSO_tools provided by Paul Cottle in _lib/, it organzied based on Python pandas module and can findout the cloud top height in Vertical Feature Mask, VFM product in CALIPSO L2 Lidar product family. Currenly, the problem it is a bit slow, and actually we do not use VFM product in this notebook, so we provide the tool, but we will not use it.
If you are interested in this tool, here is the example codes for how to use it:
# import the module
import sys
sys.path.insert(0, '_libs/CALIPSO_tools')
# search the file
hdf5_VFM=glob.glob('$YOUR_PATH/*VFM*.h5')
Cal_obj = Calipso.Calipso(h5file = hdf5_VFM[0], raw = True)
# Segment the orbit
lats = [45, 50]; lons = []; alts = [0,40000]
Cal_obj.window_select(latrange = lats,lonrange = lons,altrange = alts, inplace = True)
# Calculate Cloud top
Cal_masked = Cal_obj.feature_mask('cloud')
Cal_cloudtops = Calipso.findtops(Cal_masked)
# A simple visulation
plt.plot(Cal_cloudtops.VFM.index.get_level_values('Lats').values, Cal_cloudtops.VFM.values)
In [0]:
import h5py
import glob
import math
import scipy.io
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
% matplotlib inline
In order to read the time of orbit, here we need Python module datetime
In [0]:
import datetime
Searching all *.h5 files in _data/CALIPSO_L2/ directory, Make sure that you have the required HDF5 files exist in it.
In [0]:
hdf5_calipso=glob.glob('_data/CALIPSO_L2/*LID_L2_333mCLay*.h5')
print("CALIPSO 333m Cloud Layer files found:\n{}".format(hdf5_calipso))
In [0]:
calipso_obj=h5py.File(hdf5_calipso[0], 'r')
We have introduced how to play with HDF5 in 01_MODIS_L1B.ipynb, so here we skip that part.
Geolocation field
In [0]:
lon=calipso_obj['Longitude'][:]
lat=calipso_obj['Latitude'][:]
We need to make sure that longitude is monotonous decrease, because satellite is possible to cross the International Date Line (180E/180W line)
In [0]:
for id in range(0, len(lon)-1):
if lon[id+1] > lon[id]:
lon[id+1] = lon[id+1]-360
Currently NOT sure about height in the cell below, Lidar_data_Altitude in metadata has a shape of 583L, but the size of datasets are (55926L, 5L). I guess Lidar_data_Altitude represents the top height of every 55920/583 = 95 data points, but I don't how to make it corresponds to 5L in vertical dimension. <----------- !!!
In [0]:
height=calipso_obj['metadata'][0][-5] # <----------------- Not Sure, use the code below instead
height=np.array([0, 1, 2, 3, 4])
surf_elev=calipso_obj['DEM_Surface_Elevation'][:]
Lidar Cloud top Height
In [0]:
cloud_top_height=calipso_obj['Layer_Top_Altitude'][:];
cloud_top_height[cloud_top_height == -9999]=np.nan
Orbit time
In [0]:
profile_time=calipso_obj['Profile_Time'][:]
Now we try to plot the orbit track of the CALIPSO file based on lon, lat and profile_time
We need to know when the orbit start, and it is the 2nd element of metadata (You can check it via .keys() or other tools)
In [0]:
tai_start_string=calipso_obj['metadata'][0][1]
Then we convert the date string to datetime object.
In [0]:
orbitStart=datetime.datetime.strptime(tai_start_string, '%Y-%m-%dT%H:%M:%S.%fZ')
Create another datetime object which includes all time info. of data points.
In [0]:
time_vals=[]
#python datetime objects in utc
for the_time in profile_time:
date_time=orbitStart + datetime.timedelta(seconds=float(the_time))
time_vals.append(date_time)
date_day=time_vals
Finally, we plot the orbit with time on a map. There are more than 50000 data points in the file, so we have to skip some of them.
The code is a bit complex, here is the key idea:
time_step=1500; # plot every 1500 data points
label_step=2; # label every 2 ploted data points
for i in np.arange(0, len(date_day), time_step):
# every timestep
proj.plot(x[i], y[i], 'bo', markersize=8) # plot the data points
# labelling
if(i % (time_step*label_step) == 0): # just like "mod" in MATLAB
time_string=date_day[i].strftime('%H:%M UCT') # .strftime() is a call of output format for datetime object
axis.text(x_str[i], y_str[i], '{0:s}'.format(time_string), \ # axis.text(x, y, 'string', **keywords)
fontsize=8, fontweight='bold', ha='left', \
bbox={'edgecolor':'w', 'facecolor':'w', 'alpha':0.875, 'pad':0})
In [0]:
time_step=1500; label_step=2;
Orbit_number=calipso_obj['metadata'][0][10] # <---- Granule Number when orbit start
# locations of time label
lat_str=lat-1; lon_str=lon+4 # lat/lon for labels
# figure
fig=plt.figure(figsize=(15, 15))
axis=plt.gca()
lon_mid=0
proj = Basemap(projection='mill', resolution='c', \
llcrnrlat=-90, urcrnrlat=90, \
llcrnrlon=-240, urcrnrlon=120)
proj.drawcoastlines(color=[0.375, 0.375, 0.375])
# draw parallels and meridians.
proj.drawparallels(np.arange(-60, 90, 30), labels=[1, 0, 0, 0])
proj.drawmeridians(np.arange(0, 360, 60), labels=[0, 0, 0, 1])
# compute native x,y coordinates of grid.
x, y=proj(lon, lat)
x_str, y_str=proj(lon_str, lat_str)
# line of track
proj.plot(x, y, linewidth=3.5, color='k', linestyle='-')
# =================================== main part ================================= #
for i in np.arange(0, len(date_day), time_step):
# every timestep
proj.plot(x[i], y[i], 'bo', markersize=8)
# labelling
if(i % (time_step*label_step) == 0):
time_string=date_day[i].strftime('%H:%M UCT')
axis.text(x_str[i], y_str[i], '{0:s}'.format(time_string), \
fontsize=8, fontweight='bold', ha='left', \
bbox={'edgecolor':'w', 'facecolor':'w', 'alpha':0.875, 'pad':0})
# ================================================================================ #
# starting point in red
proj.plot(x[0], y[0], 'ro', markersize=8)
# bold axis
[i.set_linewidth(2) for i in axis.spines.itervalues()]
# title
title_str='CALIPSO Track\nGranule Number: 0' + str(Orbit_number)
axis.set_title(title_str, fontsize=14, fontweight='bold')
# savefig
plt.savefig('_figures/02_CALIPSO_track.png', dpi=450, facecolor='w', edgecolor='w',
orientation='portrait', papertype='a4', format='png',
transparent=True, bbox_inches='tight', pad_inches=0,
frameon=None)
plt.show()
In [0]:
Read the data we saved in 01_MODIS_L1B.ipynb.
Remembering that we use
scipy.io.savemat('_share/01_MODIS_L1B_TBright', {'longitude': longitude, 'latitude': latitude, 'cloud_Tbright': cloud_Tbright})
to save the file.
In [0]:
MODIS_L1B=scipy.io.loadmat('_share/01_MODIS_L1B_TBright.mat')
cloud_Tbright=MODIS_L1B['cloud_Tbright']
lon_L1B=MODIS_L1B['longitude']
lat_L1B=MODIS_L1B['latitude']
Mask the NaN value in Tbring
In [0]:
cloud_Tbright_masked=np.ma.masked_where(np.isnan(cloud_Tbright), cloud_Tbright)
We plot the Brightness Temperature to represents the clouds, and plot the orbit track. The code is similar to the last part of 01_MODIS_L1B.ipynb.
The orbit is a line, but here I plot it with super big linewidth to make it clear.
In [0]:
lonlim=[-155, -127]
latlim=[30, 60]
# Colormap
CMap=plt.cm.gray_r#gist_heat # hot, afmhot, gnuplot
#
proj=Basemap(resolution='l', projection='lcc', \
lat_1=30, lat_2=60, lat_0=45, lon_0=-140, \
llcrnrlon=-155, llcrnrlat=30, \
urcrnrlon=-110, urcrnrlat=56)
# create figure, add axes
fig=plt.figure(figsize=(12, 12))
ax=plt.gca()
## parallels and meridians.
parallels=np.arange(-90, 90, 5)
meridians=np.arange(0, 360, 5)
proj.drawparallels(parallels, labels=[1, 0, 0, 0],\
fontsize=10, latmax=90)
proj.drawmeridians(meridians, labels=[0, 0, 0, 1],\
fontsize=10, latmax=90)
# draw coast
proj.drawcoastlines(linewidth=1.5, linestyle='solid', color=[0.25, 0.25, 0.25])
# compute native x,y coordinates of grid.
x, y=proj(lon_L1B, lat_L1B)
x_track, y_track=proj(lon, lat)
# pcolor plot
CS=proj.pcolor(x, y, cloud_Tbright_masked, cmap=CMap, vmin=210, vmax=275)
# plot the track
proj.plot(x_track, y_track, linewidth=15, color='r', linestyle='-', alpha=0.625)
# title
ax.set_title('Cloud and Orbit Track',\
fontweight='bold', fontsize=14)
# Save figure
plt.savefig('_figures/02_Cloud_and_Orbit_Track.png', dpi=450, facecolor='w', edgecolor='w',
orientation='portrait', papertype='a4', format='png',
transparent=True, bbox_inches='tight', pad_inches=0,
frameon=None)
# Show
plt.show()
Based on the figure above, we can segment the orbit by latitude:
$35^\circ N < Latitude < 52^\circ N$
We use `numpy.searchsorted` to segment. This function has been discussed in 01_MODIS_L1B.ipynb.
Here, latitude is the "bins", and latitude=35, latitude=50 is the value.
In [0]:
minlat_indices=np.searchsorted(lat.flat, 35, 'right')-1
maxlat_indices=np.searchsorted(lat.flat, 52, 'right')-1
# Cut
cloud_top_height=cloud_top_height[minlat_indices:maxlat_indices]
lon_seg=lon[minlat_indices:maxlat_indices]
lat_seg=lat[minlat_indices:maxlat_indices]
surf_elev_seg=surf_elev[minlat_indices:maxlat_indices]
date_day_seg=date_day[minlat_indices:maxlat_indices]
Since the longitude of the orbit is monotone decreasing according to our discussion.
But the output of np.searchsorted is naturally in a monotone increasing sense, so here we change the "direction" of lon_seg.
In [0]:
lon_seg=lon_seg[::-1, :];
#cloud_top_height=cloud_top_height[::-1, :];
We save the result as *.mat
In [0]:
time_str=[0] * len(date_day_seg)
for i in range(len(date_day_seg)):
time_str[i]=date_day_seg[i].strftime('%H:%M UCT')
scipy.io.savemat('_share/02_CALIPSO_L2', {'longitude': lon_seg, 'latitude': lat_seg, 'time': time_str, 'surf_elev': surf_elev_seg, 'Cloud_Top_Height': cloud_top_height})
We plot the CALIPSO lidar cloud top height on the segmented orbit. Essentially, it is:
plt.plot(lat_seg, cloud_top_height)
Mask the NaN value.
In [0]:
cloud_top_height_masked=np.ma.masked_where(np.isnan(cloud_top_height), cloud_top_height)
In [0]:
from mpl_toolkits.axes_grid1 import host_subplot
import mpl_toolkits.axisartist as AA
# figure
fig=plt.figure(figsize=(10, 5))
# host_subplot
host = host_subplot(111, axes_class=AA.Axes)
plt.subplots_adjust(bottom=0.25)
# two more axis
par2 = host.twiny()
par3 = host.twiny()
# Still not so sure about how this works, currently looks good : )
offset = -25
new_fixed_axis = par2.get_grid_helper().new_fixed_axis
par2.axis["top"] = new_fixed_axis(loc="bottom", axes=par2, offset=(0, offset))
offset = -62.5
new_fixed_axis = par3.get_grid_helper().new_fixed_axis
par3.axis["top"] = new_fixed_axis(loc="bottom", axes=par3, offset=(0, offset))
# axis label <---- don't know why fontsize=14, fontweight='bold' doesn't work
host.set_ylabel('Lidar Cloud top Height (km)', fontsize=14, fontweight='bold')
host.set_xlabel('\n\nGeolocation (longitude/latitude)')
par3.set_xlabel('Orbit time')
# axis lim
par2.set_xlim(lat_seg.min(), lat_seg.max())
par3.set_xlim(0, 1)
host.set_xlim(lon_seg.min(), lon_seg.max())
host.set_ylim(-0.5, 10)
# plot
# fill the surface elevation
baseline=-0.5*np.ones(surf_elev_seg.shape).flat[:]
host.fill_between(lon_seg.flat[:], surf_elev_seg.flat[:], baseline, \
where=surf_elev_seg.flat[:]>=baseline, facecolor=[0.5, 0.5, 0.5], interpolate=False)
host.plot(lon_seg, cloud_top_height_masked[:, 0], linewidth=1.5, linestyle='-', color='k')
# xtick and ytick
host.set_yticks([1, 3, 5, 7, 9])
host.set_xticks(np.linspace(lon_seg.min(), lon_seg.max(), 6))
par2.set_xticks(np.linspace(lat_seg.min(), lat_seg.max(), 6))
par2.invert_xaxis()
count=0; labels_str=[0] * 6
for time_id in np.linspace(0, len(date_day_seg)-1, 6).astype(int):
labels_str[count]=date_day_seg[time_id].strftime('%H:%M UCT')
count+=1
par3.set_xticklabels(labels_str)
host.grid()
host.set_title('CALIPSO Lidar Cloud top Height', fontsize=16, fontweight='bold')
# savefig
plt.savefig('_figures/02_CALIPSO_Lidar_Cloud_top_Height.png', dpi=450, facecolor='w', edgecolor='w',
orientation='portrait', papertype='a4', format='png',
transparent=True, bbox_inches='tight', pad_inches=0,
frameon=None)
plt.show()
In [0]:
In [0]: