- 积分
- 25770
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2017-9-4
- 最后登录
- 1970-1-1
|
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看，没有帐号？立即注册
x
本帖最后由 灭火器 于 2020-10-20 08:04 编辑
HDF EOS 的官网上有用 Python 读取 CloudsSat 2B-GEOPROF 数据的实例脚本（https://hdfeos.org/zoo/index_openCDPC_Examples.php#CloudSAT）用该网站上提供的实例数据进行测试，能够得到正常的结果。图片如下
但是当我用自己下载的数据进行测试时，发现 CloudSat 的数据在 2019 年初升级为 R05 版本，而该脚本（包括 NCL 版）无法正常处理 R05 的数据。尝试后发现原因是 R05 中的变量不再附有自描述的属性，所以下面这样的代码将会失效：
- # Read attributes.
- attrs = dset.attributes(full=1)
- lna=attrs["long_name"]
- long_name = lna[0]
- sfa=attrs["factor"]
- scale_factor = sfa[0]
- vra=attrs["valid_range"]
- valid_min = vra[0][0]
- valid_max = vra[0][1]
- ua=attrs["units"]
- units = ua[0]
雷达反射率的剖面图中额外添上了地形（灰色部分）。结果与上面的图片一致，可以发现原图片中低空处的强回波是地形导致的。相关代码直接贴在后面：
- from pyhdf import HDF, VS, SD
- import numpy as np
- import numpy.ma as ma
- import pandas as pd
- class reader:
- '''用于读取CloudSat R05 2B产品的类.'''
- def __init__(self, fname):
- '''打开HDF,Vdata table,和scientific dataset.'''
- self.hdf = HDF.HDF(fname, HDF.HC.READ)
- self.vs = self.hdf.vstart()
- self.sd = SD.SD(fname, SD.SDC.READ)
- def attach_vdata(self, varname):
- '''从Vdata table中读取一个变量的所有记录.'''
- vdata = self.vs.attach(varname)
- data = vdata[:]
- vdata.detach()
- return data
-
- def scale_and_mask(self, data, varname):
- '''依据变量的factor进行放缩,并根据valid_range来mask掉填充值.'''
- factor = self.attach_vdata(f'{varname}.factor')[0][0]
- valid_min, valid_max = \
- self.attach_vdata(f'{varname}.valid_range')[0][0]
- invalid = np.logical_or(data <= valid_min, data >= valid_max)
- data = ma.array(data, mask=invalid)
- data = data / factor
- return data
-
- def read_geo(self, process=True):
- '''读取经纬度和地形高度.参数process指定是否放缩并mask地形高度.'''
- lon = np.array(self.attach_vdata('Longitude'))[:,0]
- lat = np.array(self.attach_vdata('Latitude'))[:,0]
- elv = np.array(self.attach_vdata('DEM_elevation'))[:,0]
- if process:
- elv = self.scale_and_mask(elv, 'DEM_elevation')
-
- return lon, lat, elv
-
- def read_time(self, datetime=True):
- '''读取每个数据点对应的时间.
-
- datetime=True: 返回所有数据点的日期时间组成的DatetimeIndex.
- datetime=False: 返回所有数据点相对于第一个点所经过的秒数组成的numpy数组.
- '''
- seconds = np.array(self.attach_vdata('Profile_time'))[:,0]
- if datetime:
- TAI = self.attach_vdata('TAI_start')[0][0]
- start = pd.to_datetime('1993-01-01') + pd.Timedelta(seconds=TAI)
- offsets = pd.to_timedelta(seconds, unit='s')
- time = pd.date_range(start=start, end=start, periods=offsets.size)
- time = time + offsets
- return time
- else:
- return seconds
-
- def read_sds(self, varname, process=True):
- '''读取scientific dataset.参数process指定是否放缩并mask.'''
- data = self.sd.select(varname)[:]
- if process:
- data = self.scale_and_mask(data, varname)
-
- return data
-
- def close(self):
- '''关闭HDF文件.'''
- self.vs.end()
- self.sd.end()
- self.hdf.close()
- # 测试数据.
- if __name__ == '__main__':
- fname = '2010128055614_21420_CS_2B-GEOPROF_GRANULE_P1_R05_E03_F00.hdf'
- f = reader(fname)
- f.close()
- #-----------------------------------------------------------
- # 2020/10/19
- # 读取CloudSat 2B-GEOPROF数据,画出卫星轨迹和雷达反射率因子的剖面图.
- # 参考:
- # [1] https://hdfeos.org/zoo/index_openCDPC_Examples.php#CloudSAT
- # [2] https://moonbooks.org/Articles/How-to-read-CloudSat-2B-GEOPROF-GRANULE-HDF4-file-using-python-and-pyhdf-/
- #-----------------------------------------------------------
- from read_CloudSat import reader
- import numpy as np
- import numpy.ma as ma
- import matplotlib.pyplot as plt
- from mpl_toolkits.axes_grid1 import make_axes_locatable
- import cartopy.crs as ccrs
- from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
- def set_map(ax):
- '''设置地图.'''
- proj = ccrs.PlateCarree()
- ax.coastlines(lw=0.5)
-
- xticks = np.arange(-180, 180 + 60, 60)
- yticks = np.arange(-90, 90 + 30, 30)
- ax.set_xticks(xticks, crs=proj)
- ax.set_yticks(yticks, crs=proj)
- ax.xaxis.set_major_formatter(LongitudeFormatter())
- ax.yaxis.set_major_formatter(LatitudeFormatter())
- ax.tick_params('both', labelsize='x-small')
- ax.set_global()
- def draw_track(ax, lon1D, lat1D):
- '''根据经纬度画出轨迹,并标识出起点.'''
- ax.plot(lon1D, lat1D, lw=2, color='b', transform=ccrs.Geodetic())
- ax.plot(lon1D[0], lat1D[0], 'ro', ms=3, transform=ccrs.PlateCarree())
- ax.text(
- lon1D[0] + 5, lat1D[0], 'start', color='r', fontsize='x-small',
- transform=ccrs.PlateCarree()
- )
- def draw_cross_section(ax, time, hgt, data):
- '''画出data的时间-高度剖面.'''
- im = ax.pcolormesh(time, hgt, data.T, cmap='jet', shading='nearest')
-
- ax.set_ylim(0, 20)
- ax.set_xlim(0, time.max())
- ax.tick_params(labelsize='x-small')
- ax.set_xlabel('Seconds since the start of the granule [s]', fontsize='x-small')
- ax.set_ylabel('Height [km]', fontsize='x-small')
- divider = make_axes_locatable(ax)
- cax = divider.append_axes('bottom', size='10%', pad=0.42)
- cbar = plt.colorbar(im, cax=cax, extend='both', orientation='horizontal')
- cbar.ax.tick_params(labelsize='x-small')
- cbar.set_label('dBZe', fontsize='x-small')
- def draw_elevation(ax, time, elv):
- '''画出地形抬升.'''
-
- ax.fill_between(time, elv, color='gray')
- if __name__ == '__main__':
- # 读取数据
- fname = '2010128055614_21420_CS_2B-GEOPROF_GRANULE_P1_R05_E03_F00.hdf'
- f = reader(fname)
- lon, lat, elv = f.read_geo()
- data = f.read_sds('Radar_Reflectivity')
- height = f.read_sds('Height')
- time = f.read_time(datetime=False)
- start_time, end_time = f.read_time(datetime=True)[[0,-1]]
- f.close()
- # 每个profile的高度有微小的差别,这里用平均值代替所有profile的高度.
- hgt = height.mean(axis=0).data
- # 高度单位全部转为km.
- elv /= 1000
- hgt /= 1000
- fig = plt.figure(dpi=200, figsize=(6, 6))
- ax1 = fig.add_axes([0.25, 0.45, 0.5, 0.5], projection=ccrs.PlateCarree())
- ax2 = fig.add_axes([0.25, 0.2, 0.5, 0.3])
- set_map(ax1)
- draw_track(ax1, lon, lat)
- # 把granule的起止时间标上.
- start_str = start_time.strftime('%Y-%m-%d %H:%M')
- end_str = end_time.strftime('%Y-%m-%d %H:%M')
- ax1.set_title(f'From {start_str} to {end_str}', fontsize='small')
- draw_cross_section(ax2, time, hgt, data)
- draw_elevation(ax2, time, elv)
- ax2.set_title('Radar Reflectivity Factor', fontsize='small')
- fig.savefig('test', bbox_inches='tight')
- plt.close(fig)
|
|
窥视卡
雷达卡
发表于 2020-10-19 15:10:09
提升卡
置顶卡
沉默卡
喧嚣卡
变色卡
抢沙发
千斤顶
显身卡
楼主
