【Python基础绘图】Python提取ICESat-2/ATL08并可视化

【Python基础绘图】Python提取ICESat-2/ATL08并可视化

01 引言：

利用Python借助h5py包从h5格式数据中提取关键字段并存为txt，并可视化，现记录分享给更多有需要的同学。

02 读取数据：

# -*- encoding: utf-8 -*-
'''
@File    :   height_data.py
@Time    :   2022/05/07 14:14:30
@Author  :   HMX 
@Version :   1.0
@Contact :   kzdhb8023@163.com
'''
​
# here put the import lib
import numpy as np
import os
import h5py
import time
​
def hdf2arr(item, objprop):
    dset_name = f'/{item}/land_segments/'+objprop
    datavar = f[dset_name]
    data = datavar[:]
    units = datavar.attrs['units']
    long_name = datavar.attrs['long_name']
    arr=np.array(data)
    return arr
​
def hdf2arr2(item, objprop):
    dset_name = f'/{item}/land_segments/canopy/'+objprop
    datavar = f[dset_name]
    data = datavar[:]
    units = datavar.attrs['units']
    long_name = datavar.attrs['long_name']
    arr=np.array(data)
    return arr
​
def hdf2arr3(item, objprop):
    dset_name = f'/orbit_info/'+objprop
    datavar = f[dset_name]
    data = datavar[:]
    units = datavar.attrs['units']
    long_name = datavar.attrs['long_name']
    # arr=np.array(data)
    return data
​
​
t1 = time.time()
​
f = open (r'.\ATL08_V5.txt','a')
# 读取的关键字段
headlist = [
'lon',
'lat',
'cloud_flag_atm',
'night_flag',
'urban_flag',
'segment_landcover',
'h_mean_canopy',
'h_max_canopy',
'sc_orient',
'track',
'date',
]
f.write(','.join(headlist)+'\n')
f.close()
h5_path = r'Z:\Dataset\Canopy_Height\ATL08_V5\demo'
flist = os.listdir(h5_path)
for file in  flist:
    fn = os.path.join(h5_path,file)
    for item in ['gt3r','gt2r','gt1r','gt3l','gt2l','gt1l']:#循环读取多个激光带
        with h5py.File(fn, mode='r') as f:
            try:
                latvar = f[f'{item}/land_segments/latitude']
                latitude = latvar[:]             
                lonvar = f[f'{item}/land_segments/longitude']
                longitude = lonvar[:]
                lon=np.array(longitude)
                lat=np.array(latitude)
                arr1 = hdf2arr(item, objprop = 'cloud_flag_atm')
                arr2 = hdf2arr(item, objprop = 'night_flag')
                arr3 = hdf2arr(item, objprop = 'urban_flag')
                arr4 = hdf2arr(item, objprop = 'segment_landcover') 
                arr5 = hdf2arr2(item, objprop = 'h_mean_canopy')
                arr6 = hdf2arr2(item, objprop = 'h_max_canopy')              
                arr7 = np.asarray([str(hdf2arr3(item, objprop = 'sc_orient')[0])]*(lon.shape[0]))
                arr8 = np.asarray([item]*(lon.shape[0]))
                namelist=[str(os.path.basename(file))[6:14]]*(lon.shape[0])
                f = open (r'.\ATL08_V5.txt','a')
                for i in range(lon.shape[0]):
                    if (lon[i]>=71)&(lon[i]<=137)&(lat[i]<=56)&(lat[i]>=1)&(arr6[i]<999.):#初步筛选范围并剔除无效值
                        outlist = [str(lon[i]),str(lat[i]),str(arr1[i]),str(arr2[i]),str(arr3[i]),str(arr4[i])
                        ,str(arr5[i]),str(arr6[i]),str(arr7[i]),str(arr8[i]),str(namelist[i])]
                        f.write(','.join(outlist)+'\n')
            except KeyError:
                print(file)
                pass
            continue
f.close()
t2 = time.time()
print('共计用时{:.2f}s'.format(t2-t1))

03 结果如下：

04 可视化：

# -*- encoding: utf-8 -*-
'''
@File    :   plot_atl08.py
@Time    :   2022/05/07 14:46:23
@Author  :   HMX 
@Version :   1.0
@Contact :   kzdhb8023@163.com
'''
​
# here put the import lib
import pandas as pd
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.ticker import LatitudeFormatter, LongitudeFormatter
from cnmaps import get_adm_maps
​
def cm2inch(x,y):
    return(x/2.54,y/2.54)
​
size = 10.5
fontdict = {'weight':'bold','size':size,'family':'SimHei'}
mpl.rcParams.update({
    'text.usetex': False,
    'font.family': 'stixgeneral',
    'mathtext.fontset': 'stix',
    "font.family":'serif',
    "font.size": size,
    "mathtext.fontset":'stix',
    "font.serif": ['Times New Roman'],
})
​
df = pd.read_csv(r'.\ATL08_V5.txt',header=0)
​
proj = ccrs.PlateCarree()
fig,ax = plt.subplots(figsize = cm2inch(8,6),subplot_kw={'projection':proj})
region = [71, 137, 1, 56]
ax.set_extent(region, crs=proj)
​
china, sourth_sea = get_adm_maps(level='国', only_polygon=True)
ax.add_geometries(china, crs=ccrs.PlateCarree(), edgecolor='k', facecolor='none')
ax.add_geometries(sourth_sea, crs=ccrs.PlateCarree(), edgecolor='k')
ax.add_feature(cfeature.COASTLINE.with_scale('10m'))
ax.add_feature(cfeature.LAND.with_scale('10m'))
ax.add_feature(cfeature.OCEAN.with_scale('10m'))
# 由于数据字段较多仅展示h_mean_canopy
img= ax.scatter(df.lon,df.lat,c = df.h_max_canopy, s=1, cmap=plt.cm.jet, transform=ccrs.PlateCarree())
​
ax.set_xticks(np.arange(region[0],region[1]+1,22), crs = proj)
ax.set_yticks(np.arange(region[2],region[3]+1,11), crs = proj)
ax.xaxis.set_major_formatter(LongitudeFormatter())
ax.yaxis.set_major_formatter(LatitudeFormatter())
ax.minorticks_on()
​
cb = plt.colorbar(img, pad=0.05,shrink =.8,aspect = 30)
cb.ax.set_yticks(np.arange(10,60,10))
cb.ax.set_ylabel('height/m')
​
plt.tight_layout()
plt.savefig(r'height.png',dpi = 600)
plt.show()

05 可视化结果如下：

​如果对你有帮助的话，请‘点赞’、‘收藏’，‘关注’，你们的支持是我更新的动力。
欢迎关注公众号【森气笔记】。