爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

搜索
查看: 6896|回复: 1

读取及绘制CALIOP L1B total attenuated backscatter signal

[复制链接]
发表于 2022-2-7 12:30:49 | 显示全部楼层 |阅读模式

登录后查看更多精彩内容~

您需要 登录 才可以下载或查看,没有帐号?立即注册

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




密码修改失败请联系微信:mofangbao
发表于 2022-2-8 11:03:09 | 显示全部楼层
大神,新年吉祥!
感谢分享,加入了对时间的循环,批量出图的功能,酷毙了。
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

快速回复 返回顶部 返回列表