#?liding?2019.12.30
#+++++++++++++++讀取MAIAC.hdf數(shù)據(jù)++++++++++++++++++++++
#該數(shù)據(jù)比較怪異，一個(gè)是組織結(jié)構(gòu)
#?（MAIAC的反演算法造成的，一個(gè)數(shù)據(jù)里有多軌數(shù)據(jù)，這里使用了均值）
#另一個(gè)是投影，無(wú)經(jīng)緯度信息?（使用了投影包）
#+++++++++++++++讀取MAIAC.hdf數(shù)據(jù)+++++++++++++++++++++++
#from?osgeo?import?gdalosr
import?pyproj
from?pathlib?import?Path
import?xarray?as?xr
import?numpy?as?np
in_dir?=?Path（‘./dat‘）
ot_dir?=?Path（‘./result‘）
k?=?0
for?i?in?in_dir.glob（‘*.hdf‘）:
????inp?=?i
????print（inp）
????with?xr.open_dataset（inpengine?=?‘netcdf4‘）?as?ds:
????????AOD_1k?=?ds[‘Optical_Depth_055‘].values??#讀取數(shù)值
????????attr_?=?ds.attrs??#讀取坐標(biāo)
????#數(shù)據(jù)是多軌的，需要進(jìn)行平均
????aod_550?=?np.full（[AOD_1k.shape[1]AOD_1k.shape[2]]np.nan）
????for?ii?in?range（AOD_1k.shape[1]）:
????????for?jj?in?range（AOD_1k.shape[2]）:
????????????aod_550[iijj]?=?np.nanmean（AOD_1k[:iijj]）?#需要使用避免nan的均值
????#投影是GCTP_SNSOID的，需要轉(zhuǎn)化，只有左上角和右下角坐標(biāo)
????for?key_val?in?attr_.items（）:
????????if?key_?==?‘Structmetadata.0‘:
????????????s_val?=?val.split（）
????for?ii?in?s_val:
????????if?‘pperLeftPointMtrs‘?in?ii:
????????????x_up?=?float（ii[ii.index（‘（‘）+1:ii.index（‘‘）]）
????????????y_up?=?float（ii[ii.index（‘‘）+1:-1]）
????????if?‘LowerRightMtrs‘?in?ii:
????????????x_down?=?float（ii[ii.index（‘（‘）+1:ii.index（‘‘）]）
????????????y_down?=?float（ii[ii.index（‘‘）+1:-1]）
????#坐標(biāo)網(wǎng)格化
????grid_xgrid_y?=?np.meshgrid（np.linspace（x_upx_downAOD_1k.shape[1]）np.linspace（y_upy_downAOD_1k.shape[2]））
????p?=?pyproj.Proj（‘+proj=sinu?+lon_0=0?+x_0=0?+y_0=0?+a=6371007.181?+b=6371007.181?+units=m?+no_defs‘）??#modis坐標(biāo)投影系統(tǒng)（可替換，還有2個(gè)一樣的），查找https://spatialreference.org/ref/
????lat?=?np.full（grid_x.shapenp.nan）?#經(jīng)緯度空矩陣，便于存放數(shù)據(jù)
????lon?=?np.full（grid_x.shapenp.nan）
????#循環(huán)計(jì)算每個(gè)格點(diǎn)的投影，注意檢查是否符合規(guī)律，避免弄反
????for?l?in?range（lat.shape[0]）:
????????for?r?in?range（lat.shape[1]）:
????????????lon[lr]?lat[lr]?=?p（grid_x[lr]?grid_y[lr]?inverse=True）
????#保存（注意名字）
????ot_file?=?ot_dir/（i.name[9:16]+‘_‘+i.name[-8:-5]）
????np.savez（ot_filelon=lonlat=lataod550=aod_550）

?屬性????????????大小?????日期????時(shí)間???名稱
-----------?---------??----------?-----??----
?????文件????????2342??2020-05-22?10:09??read_MAIAC.py