- 积分
- 25
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2012-2-16
- 最后登录
- 1970-1-1
|
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看,没有帐号?立即注册
x
网上目前并没有python读取及绘制 CALIOP L1B total attenuated backscatter signal的脚本, 因此与大家分享,希望能帮助到有需要的朋友们,有什么建议可以改进脚本也欢迎大家一起讨论!
以下脚本调的colorbar 及数值范围已与官网一致,
也画了一个官网展示的例子,也就是2010-04-15。
同时这个脚本也做了对时间的循环,方便大家批量出图。
同时还可对区域进行限定,也方便大家选取自己关注的区域。
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
import matplotlib.dates as mdates
from matplotlib.dates import DateFormatter
import matplotlib.ticker as ticker
import matplotlib.gridspec as gridspec
import matplotlib.ticker as mticker
from matplotlib.ticker import AutoMinorLocator
from matplotlib.colors import BoundaryNorm, ListedColormap
import glob
import sys
import xarray as xr
import h5py
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
from calendar import monthrange
def number_of_days_in_month(year=2020, month=3):
return monthrange(year, month)[1]
def find_nearest(array, value):
array = np.asarray(array)
idx = (np.abs(array - value)).argmin()
return array[idx],idx
#!/usr/bin/env python
from pyhdf.SD import SD, SDC # To read a HDF file (SD data)
from pyhdf.HDF import *
from pyhdf.VS import *
import datetime
def calc_dayofYear(year=2015,month=4,day=24):
day_of_year = datetime.datetime(year, month, day).strftime('%j')
return day_of_year
import pandas as pd
import matplotlib.colors as mcolors
cmap2 = ListedColormap(["navy","navy","royalblue","royalblue","royalblue","royalblue","royalblue","lime","green",\
"limegreen","yellow","yellow","gold","orange","darkorange","tomato","red","deeppink","hotpink","pink",\
"dimgrey","grey","darkgrey","silver","lightgrey","gainsboro","whitesmoke","white","snow",\
"snow","snow","snow","snow","snow"])
bounds = [1.0e-4,2.0e-4,3.0e-4,4.0e-4,5.0e-4,6.0e-4,7.0e-4,8.0e-4,9.0e-4,\
1.0e-3,1.5e-3,2.0e-3,2.5e-3,3.0e-3,3.5e-3,4.0e-3,4.5e-3,5.0e-3,5.5e-3,6.0e-3,6.5e-3,7.0e-3,\
7.5e-3,8.0e-3,\
1.0e-2,2.0e-2,3.0e-2,4.0e-2,5.0e-2,6.0e-2,7.0e-2,8.0e-2,9.0e-2,1.0e-1]
norm = BoundaryNorm(bounds, cmap2.N)
year_start, year_stop =2010, 2010
month_start,month_stop= 4, 4
day_start,day_stop = 15, 15
year_list = np.arange(start=year_start,stop=year_stop+1,step=1)
month_list= np.arange(start=month_start,stop=month_stop+1,step=1)
customize_day = True
# select zonal data
latn, lats = -40, -60
lonw, lone = -180, 180
for indyr in year_list:
for indmonth in month_list:
if customize_day:
for indday in np.arange(start=day_start,stop=day_stop+1,step=1):
tmp_file = sorted(glob.glob('/your_directory/CAL_LID_L1-Standard-V4-10.'+\
str(indyr).strip().zfill(4)+'-'+str(indmonth).strip().zfill(2)+'-'+str(indday).strip().zfill(2)+'*.hdf'))
if np.size(tmp_file) > 0:
for ilin, FILE_NAME in enumerate(tmp_file):
try:
# open the hdf file for reading
# extract altitude infomation from metadata
hdf_interface = HDF(FILE_NAME)
vs_interface = hdf_interface.vstart()
meta = vs_interface.attach("metadata")
field_infos = meta.fieldinfo()
all_data = meta.read(meta._nrecs)[0]
meta.detach()
data_dictionary = {}
field_name_index = 0
for field_info, data in zip(field_infos, all_data):
data_dictionary[field_info[field_name_index]] = data
lidar_altitudes = data_dictionary["Lidar_Data_Altitudes"]
vs_interface.end()
hdf_interface.close()
calipso_lidar_hgt = np.array(lidar_altitudes)
#---------- Read HDF Files ----------#
hdff = SD(FILE_NAME, SDC.READ)
tbc = hdff.select('Total_Attenuated_Backscatter_532')
lat = hdff.select('Latitude')
lon = hdff.select('Longitude')
profile_time = hdff.select('Profile_Time')
#---------- retrieve data ----------#
tbc_data = tbc.get()
lat_data = lat.get()
lon_data = lon.get()
ptime = profile_time.get()
lat_data = lat_data.flatten()
lon_data = lon_data.flatten()
ptime = ptime.flatten()
human_read_timetbcs = datetime.datetime(1993,1,1,0,0,0) + datetime.timedelta(seconds=ptime[0])
human_read_timetbce = datetime.datetime(1993,1,1,0,0,0) + datetime.timedelta(seconds=ptime[-1])
tmp_idx = np.where((lat_data <= latn) & (lat_data >= lats))[0]
ma_tbc_data = np.where(((tbc_data == -9999) | (tbc_data == -333)), np.nan, tbc_data)
xidx = np.arange(start=0,stop=len(tmp_idx),step=1)
if np.size(xidx > 0):
# plot ----------------------------------------------------------------------
dpi = 100
fig = plt.figure(figsize=(2000/dpi,1000/dpi), dpi=dpi) ##,constrained_layout=True)
spec= gridspec.GridSpec(ncols=1, nrows=1,figure=fig,left=0.1, bottom=0.3,right=0.9,top=0.8)
spec.update(wspace=0.18,hspace=0.0)
ax = fig.add_subplot(spec[0, 0])
cf = ax.pcolormesh(lat_data[tmp_idx], calipso_lidar_hgt,ma_tbc_data.T[:,tmp_idx],
cmap=cmap2,norm=norm,shading='auto')
ax.set_ylim(-2,30)
ax.set_xlim(lat_data[tmp_idx[0]], lat_data[tmp_idx[-1]])
ax.yaxis.set_minor_locator(AutoMinorLocator(5))
ax.set_ylabel('Height (km)',fontsize=14)
ax.xaxis.set_ticks(lat_data[tmp_idx][::1000])
ax.tick_params(axis='x', rotation=0)
xticklabels = []
lon_formatter = LongitudeFormatter()
lat_formatter = LatitudeFormatter()
for ii in range(0,len(lat_data[tmp_idx][:]),1000):
lon_str = lon_formatter(round(lon_data[tmp_idx][ii],1))
lat_str = lat_formatter(round(lat_data[tmp_idx][ii],1))
xticklabels.append(lat_str + '\n' + lon_str)
ax.set_xticklabels(xticklabels,)
ax.tick_params(labelsize=5)
title_time1 = human_read_timetbcs
title_time2 = human_read_timetbce
ax.set_title('Total Backscatter @ 532nm '+ title_time1.strftime("%Y-%m-%d %H:%M")+' - '+title_time2.strftime("%Y-%m-%d %H:%M UTC") ,loc='center',fontsize=16) #,fontweight='bold')
cbar = fig.colorbar(cf, ax=[ax],orientation='vertical', pad=0.02, shrink=0.85,aspect=25,
ticks=bounds[::1])
yticklabels = []
for ii in bounds[::1]:
if (ii == 1.0e-4) or (ii == 1.0e-3) or (ii == 1.0e-2) or (ii == 1.0e-1):
yticklabels.append("{:.1e}".format(ii))
elif (ii <= 9.0e-4) and (ii > 1.0e-4):
yticklabels.append(round(ii*1.0e4,1))
elif (ii <= 8.0e-3) and (ii > 1.0e-3):
yticklabels.append(round(ii*1.0e3,1))
elif (ii > 1.0e-2) and (ii <= 9.0e-2):
yticklabels.append(round(ii*1.0e2,1))
cbar.ax.set_yticklabels(yticklabels)
cbar.set_label('km$^{-1}$ sr$^{-1}$',size=8) #-- add colorbar title string
cbar.ax.tick_params(axis='y',direction='in',labelsize=8,rotation=0)
plt.show()
plt.close(fig)
except Exception as e:
print(e)
sys.exit()
else:
pass
for indday in np.arange(start=1,stop=number_of_days_in_month(year=indyr, month=indmonth)+1,step=1):
pass
|
|