完整代碼見(jiàn) Github：https://github.com/ChaoQiezi/read_fy3d_tshs，由于代碼中注釋較為詳細(xì)，因此博客中部分操作一筆帶過(guò)。

01 FY3D的HDF轉(zhuǎn)TIFF

1.1 數(shù)據(jù)集說(shuō)明

FY3D TSHS數(shù)據(jù)集是二級(jí)產(chǎn)品(TSHS即MWTS/MWHS 融合大氣溫濕度廓線/穩(wěn)定度指數(shù)/位勢(shì)高度產(chǎn)品)；
文件格式為HDF5；
空間分辨率為33KM（星下點(diǎn)）；
范圍為全球區(qū)域；(FY3D為極軌衛(wèi)星，因此對(duì)于得到的單個(gè)數(shù)據(jù)集并沒(méi)有完全覆蓋全球區(qū)域)，類似下方（其余地區(qū)均未掃描到）：

提取數(shù)據(jù)涉及的多個(gè)數(shù)據(jù)集：

位于DATA/MWHS_Ch_BT和DATA/MWTS_Ch_BT 路徑下的兩個(gè)數(shù)據(jù)集分別為MWTS 通道亮溫和MWHS 通道亮溫，這是需要提取的數(shù)據(jù)集；
另外位于GEO/Latitude 和GEO/Longitude 路徑下的兩個(gè)數(shù)據(jù)集分別經(jīng)緯度數(shù)據(jù)集，用于確定像元的地理位置以便后續(xù)進(jìn)行重投影；

1.1.1 MWTS 通道亮溫和MWHS 通道亮溫?cái)?shù)據(jù)集的說(shuō)明

這是關(guān)于數(shù)據(jù)集的基本屬性說(shuō)明：

在產(chǎn)品說(shuō)明書中關(guān)于數(shù)據(jù)集的說(shuō)明如下：

可以發(fā)現(xiàn)，其共有三個(gè)維度，譬如MWHS_Ch_BT數(shù)據(jù)集的Shape為(15, 1212, 90)，其表示極軌衛(wèi)星的傳感器15個(gè)通道(即波段數(shù))在每條掃描線總共 90 個(gè)觀測(cè)像元的MWHS亮溫值，此處掃描線共有1212個(gè)。對(duì)于MWTS數(shù)據(jù)集類似格式。（這意味著三個(gè)維度中并沒(méi)有空間上的關(guān)系）

1.1.2 經(jīng)緯度數(shù)據(jù)集的說(shuō)明

這是官方產(chǎn)品說(shuō)明對(duì)于經(jīng)緯度數(shù)據(jù)集的介紹：

可以發(fā)現(xiàn)，經(jīng)緯度數(shù)據(jù)集的Shape均為(1212, 90)，這正好對(duì)應(yīng)前文提及的兩個(gè)數(shù)據(jù)集的所有像元（除去波段數(shù)），其表示每條1212條掃描線上的90個(gè)觀測(cè)像元的經(jīng)緯度。

1.2 讀取HDF5文件數(shù)據(jù)集

def read_h5(hdf_path, ds_path, scale=True):"""讀取指定數(shù)據(jù)集并進(jìn)行預(yù)處理:param hdf_path: 數(shù)據(jù)集所在HDF文件的絕對(duì)路徑:param ds_path: 數(shù)據(jù)集在HDF文件內(nèi)的路徑:return: 返回處理好的數(shù)據(jù)集"""with h5py.File(hdf_path) as f:# 讀取目標(biāo)數(shù)據(jù)集屬矩陣和相關(guān)屬性ds = f[ds_path]ds_values = np.float32(ds[:])  # 獲取數(shù)據(jù)集valid_range = ds.attrs['valid_range']  # 獲取有效范圍slope = ds.attrs['Slope'][0]  # 獲取斜率(類似scale_factor)intercept = ds.attrs['Intercept'][0]  # 獲取截距(類似add_offset)""""原始數(shù)據(jù)集可能存在縮放(可能是為存儲(chǔ)空間全部存為整數(shù)(需要通過(guò)斜率和截距再還原為原來(lái)的值,或者是需要進(jìn)行單位換算甚至物理量的計(jì)算例如最常見(jiàn)的DN值轉(zhuǎn)大氣層表觀反射率等(這多出現(xiàn)于一級(jí)產(chǎn)品的輻射定標(biāo), 但二級(jí)產(chǎn)品可能因?yàn)閱挝粨Q算等也可能出現(xiàn)));如果原始數(shù)據(jù)集不存在縮放, 那么Slope屬性和Intercept屬性分別為1和0;這里為了確保所有迭代的HDF文件的數(shù)據(jù)集均正常得到, 這里依舊進(jìn)行slope和intercept的讀取和計(jì)算(盡管可能冗余)"""# 目標(biāo)數(shù)據(jù)集的預(yù)處理(無(wú)效值, 范圍限定等)ds_values[(ds_values < valid_range[0]) | (ds_values > valid_range[1])] = np.nanif scale:ds_values = ds_values * slope + intercept  # 還原縮放"""Note: 這里之所以選擇是否進(jìn)行縮放, 原因?yàn)榻?jīng)緯度數(shù)據(jù)集中的slope為1, intercept也為1, 但是進(jìn)行縮放后超出地理范圍1°即出現(xiàn)了90.928對(duì)于緯度。其它類似, 因此認(rèn)為這里可能存在問(wèn)題如果進(jìn)行縮放, 所以對(duì)于經(jīng)緯度數(shù)據(jù)集這里不進(jìn)行縮放"""return ds_values

上述代碼用于讀取指定HDF5文件的指定數(shù)據(jù)集的數(shù)組/矩陣，scale參數(shù)用于是否對(duì)數(shù)據(jù)集進(jìn)行slope和intercept線性轉(zhuǎn)換。

1.3 重組和重投影

這一部分是整個(gè)數(shù)據(jù)處理的核心。

1.3.1 重組

def reform_ds(ds, lon, lat, reform_range=None):"""重組數(shù)組:param ds: 目標(biāo)數(shù)據(jù)集(三維):param lon: 對(duì)應(yīng)目標(biāo)數(shù)據(jù)集的經(jīng)度數(shù)據(jù)集():param lat: 對(duì)應(yīng)目標(biāo)數(shù)據(jù)集的緯度數(shù)據(jù)集(二維):param reform_range: 重組范圍, (lon_min, lat_max, lon_max, lat_min), 若無(wú)則使用全部數(shù)據(jù):return: 以元組形式依次返回: 重組好的目標(biāo)數(shù)據(jù)集, 經(jīng)度數(shù)據(jù)集, 緯度數(shù)據(jù)集(均為二維數(shù)組)"""# 裁選指定地理范圍的數(shù)據(jù)集if reform_range:lon_min, lat_max, lon_max, lat_min = reform_rangex, y = np.where((lon > lon_min) & (lon < lon_max) & (lat > lat_min) & (lat < lat_max))ds = ds[:, x, y]lon = lon[x, y]lat = lat[x, y]else:ds = ds.reshape(ds.shape[0], -1)lon = lon.flatten()lat = lat.flatten()# 無(wú)效值去除(去除地理位置為無(wú)效值的元素)valid_pos = ~np.isnan(lat) & ~np.isnan(lon)ds = ds[:, valid_pos]lon = lon[valid_pos]lat = lat[valid_pos]# 重組數(shù)組的初始化bands = []for band in ds:reform_ds_size = np.int32(np.sqrt(band.size))  # int向下取整band = band[:reform_ds_size ** 2].reshape(reform_ds_size, reform_ds_size)bands.append(band)else:lon = lon[:reform_ds_size ** 2].reshape(reform_ds_size, reform_ds_size)lat = lat[:reform_ds_size ** 2].reshape(reform_ds_size, reform_ds_size)bands = np.array(bands)return bands, lon, lat

這部分是對(duì)原始數(shù)據(jù)集(此處理中待重組數(shù)據(jù)集shape為(15, 1212, 90))進(jìn)行重組。
重組原理即基于經(jīng)緯度數(shù)據(jù)集，依據(jù)裁剪范圍將滿足地理范圍內(nèi)的所有有效像元一維化，然后重新reshape為二維數(shù)組，數(shù)組行列數(shù)均為原維度的平方根。

1.3.2 重投影

def data_glt(out_path, src_ds, src_x, src_y, out_res, zoom_scale=6, glt_range=None, windows_size=9):"""基于經(jīng)緯度數(shù)據(jù)集對(duì)目標(biāo)數(shù)據(jù)集進(jìn)行GLT校正/重投影(WGS84), 并輸出為TIFF文件:param out_path: 輸出tiff文件的路徑:param src_ds: 目標(biāo)數(shù)據(jù)集:param src_x: 對(duì)應(yīng)的橫軸坐標(biāo)系(對(duì)應(yīng)地理坐標(biāo)系的經(jīng)度數(shù)據(jù)集):param src_y: 對(duì)應(yīng)的縱軸坐標(biāo)系(對(duì)應(yīng)地理坐標(biāo)系的緯度數(shù)據(jù)集):param out_res: 輸出分辨率(單位: 度/°):param zoom_scale::return: None"""if glt_range:# lon_min, lat_max, lon_max, lat_min = -180.0, 90.0, 180.0, -90.0lon_min, lat_max, lon_max, lat_min = glt_rangeelse:lon_min, lat_max, lon_max, lat_min = np.nanmin(src_x), np.nanmax(src_y), \np.nanmax(src_x), np.nanmin(src_y)zoom_lon = zoom(src_x, (zoom_scale, zoom_scale), order=0)  # 0為最近鄰插值zoom_lat = zoom(src_y, (zoom_scale, zoom_scale), order=0)# # 確保插值結(jié)果正常# zoom_lon[(zoom_lon < -180) | (zoom_lon > 180)] = np.nan# zoom_lat[(zoom_lat < -90) | (zoom_lat > 90)] = np.nanglt_cols = np.ceil((lon_max - lon_min) / out_res).astype(int)glt_rows = np.ceil((lat_max - lat_min) / out_res).astype(int)deal_bands = []for src_ds_band in src_ds:glt_ds = np.full((glt_rows, glt_cols), np.nan)glt_lon = np.full((glt_rows, glt_cols), np.nan)glt_lat = np.full((glt_rows, glt_cols), np.nan)geo_x_ix, geo_y_ix = np.floor((zoom_lon - lon_min) / out_res).astype(int), \np.floor((lat_max - zoom_lat) / out_res).astype(int)glt_lon[geo_y_ix, geo_x_ix] = zoom_longlt_lat[geo_y_ix, geo_x_ix] = zoom_latglt_x_ix, glt_y_ix = np.floor((src_x - lon_min) / out_res).astype(int), \np.floor((lat_max - src_y) / out_res).astype(int)glt_ds[glt_y_ix, glt_x_ix] = src_ds_band# write_tiff('H:\\Datasets\\Objects\\ReadFY3D\\Output\\py_lon.tiff', [glt_lon],#            [lon_min, out_res, 0, lat_max, 0, -out_res])# write_tiff('H:\\Datasets\\Objects\\ReadFY3D\\Output\\py_lat.tiff', [glt_lat],#            [lon_min, out_res, 0, lat_max, 0, -out_res])# # 插值# interpolation_ds = np.full_like(glt_ds, fill_value=np.nan)# jump_size = windows_size // 2# for row_ix in range(jump_size, glt_rows - jump_size):#     for col_ix in range(jump_size, glt_cols - jump_size):#         if ~np.isnan(glt_ds[row_ix, col_ix]):#             interpolation_ds[row_ix, col_ix] = glt_ds[row_ix, col_ix]#             continue#         # 定義當(dāng)前窗口的邊界#         row_start = row_ix - jump_size#         row_end = row_ix + jump_size + 1  # +1 因?yàn)榍衅话Y(jié)束索引#         col_start = col_ix - jump_size#         col_end = col_ix + jump_size + 1#         rows, cols = np.ogrid[row_start:row_end, col_start:col_end]#         distances = np.sqrt((rows - row_ix) ** 2 + (cols - col_ix) ** 2)#         window_ds = glt_ds[(row_ix - jump_size):(row_ix + jump_size + 1),#                     (col_ix - jump_size):(col_ix + jump_size + 1)]#         if np.sum(~np.isnan(window_ds)) == 0:#             continue#         distances_sort_pos = np.argsort(distances.flatten())#         window_ds_sort = window_ds[np.unravel_index(distances_sort_pos, distances.shape)]#         interpolation_ds[row_ix, col_ix] = window_ds_sort[~np.isnan(window_ds_sort)][0]deal_bands.append(glt_ds)# print('處理波段: {}'.format(len(deal_bands)))# if len(deal_bands) == 6:#     breakwrite_tiff(out_path, deal_bands, [lon_min, out_res, 0, lat_max, 0, -out_res])# write_tiff('H:\\Datasets\\Objects\\ReadFY3D\\Output\\py_underlying.tiff', [interpolation_ds], [lon_min, out_res, 0, lat_max, 0, -out_res])# write_tiff('H:\\Datasets\\Objects\\ReadFY3D\\Output\\py_lon.tiff', [glt_lon], [x_min, out_res, 0, y_max, 0, -out_res])# write_tiff('H:\\Datasets\\Objects\\ReadFY3D\\Output\\py_lat.tiff', [glt_lat], [x_min, out_res, 0, y_max, 0, -out_res])

這是對(duì)重組后的數(shù)組進(jìn)行重投影，其基本思路就是對(duì)經(jīng)緯度數(shù)據(jù)集進(jìn)行zoom (congrid)，將其重采樣放大譬如此處為原行列數(shù)的六倍。
接著再將zoom后的經(jīng)緯度數(shù)據(jù)集按照角點(diǎn)信息套入到輸出glt數(shù)組中，而對(duì)重組后的目標(biāo)數(shù)組直接套入，無(wú)需進(jìn)行zoom操作。
接著對(duì)套入的目標(biāo)數(shù)組進(jìn)行最近鄰插值，如果沒(méi)有進(jìn)行插值，情況如下：

進(jìn)行最近鄰插值后（9*9窗口內(nèi)的最近有效像元為填充值），結(jié)果如下：

但是由于數(shù)據(jù)集巨大，關(guān)于最近鄰插值此處進(jìn)行了注釋，將該操作移至后面鑲嵌操作之后，數(shù)據(jù)量大大減少，所花費(fèi)時(shí)間成本也極大降低。

2 鑲嵌和最近鄰插值

2.1 鑲嵌

def img_mosaic(mosaic_paths: list, return_transform: bool = True, mode: str = 'last'):"""該函數(shù)用于對(duì)列表中的所有TIFF文件進(jìn)行鑲嵌:param mosaic_paths: 多個(gè)TIFF文件路徑組成的字符串列表:param return_transform: 是否一同返回仿射變換:param mode: 鑲嵌模式, 默認(rèn)是Last(即如果有存在像元重疊, mosaic_paths中靠后影像的像元將覆蓋其),可選: last, mean, max, min:return:"""# 獲取鑲嵌范圍x_mins, x_maxs, y_mins, y_maxs = [], [], [], []for mosaic_path in mosaic_paths:ds = gdal.Open(mosaic_path)x_min, x_res, _, y_max, _, y_res_negative = ds.GetGeoTransform()x_size, y_size = ds.RasterXSize, ds.RasterYSizex_mins.append(x_min)x_maxs.append(x_min + x_size * x_res)y_mins.append(y_max+ y_size * y_res_negative)y_maxs.append(y_max)else:y_res = -y_res_negativeband_count = ds.RasterCountds = Nonex_min, x_max, y_min, y_max = min(x_mins), max(x_maxs), min(y_mins), max(y_maxs)# 鑲嵌col = ceil((x_max - x_min) / x_res).astype(int)row = ceil((y_max - y_min) / y_res).astype(int)mosaic_imgs = []  # 用于存儲(chǔ)各個(gè)影像for ix, mosaic_path in enumerate(mosaic_paths):mosaic_img = np.full((band_count, row, col), np.nan)  # 初始化ds = gdal.Open(mosaic_path)ds_bands = ds.ReadAsArray()# 計(jì)算當(dāng)前鑲嵌范圍start_row = floor((y_max - (y_maxs[ix] - x_res / 2)) / y_res).astype(int)start_col = floor(((x_mins[ix] + x_res / 2) - x_min) / x_res).astype(int)end_row = (start_row + ds_bands.shape[1]).astype(int)end_col = (start_col + ds_bands.shape[2]).astype(int)mosaic_img[:, start_row:end_row, start_col:end_col] = ds_bandsmosaic_imgs.append(mosaic_img)# 判斷鑲嵌模式if mode == 'last':mosaic_img = mosaic_imgs[0].copy()for img in mosaic_imgs:mask = ~np.isnan(img)mosaic_img[mask] = img[mask]elif mode == 'mean':mosaic_imgs = np.asarray(mosaic_imgs)mask = np.isnan(mosaic_imgs)mosaic_img = np.ma.array(mosaic_imgs, mask=mask).mean(axis=0).filled(np.nan)elif mode == 'max':mosaic_imgs = np.asarray(mosaic_imgs)mask = np.isnan(mosaic_imgs)mosaic_img = np.ma.array(mosaic_imgs, mask=mask).max(axis=0).filled(np.nan)elif mode == 'min':mosaic_imgs = np.asarray(mosaic_imgs)mask = np.isnan(mosaic_imgs)mosaic_img = np.ma.array(mosaic_imgs, mask=mask).min(axis=0).filled(np.nan)else:raise ValueError('不支持的鑲嵌模式: {}'.format(mode))if return_transform:return mosaic_img, [x_min, x_res, 0, y_max, 0, -y_res]return mosaic_img

這里的鑲嵌不僅僅可以解決相同地理位置的鑲嵌，也可以解決不同地理位置的拼接，拼接方式支持最大最小值計(jì)算、均值計(jì)算、Last等模式。
思路非常簡(jiǎn)單，就是將輸入的每個(gè)數(shù)據(jù)集重新裝入統(tǒng)一地理范圍的箱子（數(shù)組）中，使所有數(shù)組的角點(diǎn)信息一致，然后基于數(shù)據(jù)集個(gè)數(shù)這一維度進(jìn)行Mean、Max、Min等計(jì)算，得到鑲嵌數(shù)組。

2.2 最近鄰插值

def window_interp(arr, distances):if np.sum(~np.isnan(arr)) == 0:return np.nan# 距離最近的有效像元arr = arr.flatten()arr_sort = arr[np.argsort(distances)]if np.sum(~np.isnan(arr_sort)) == 0:return np.nanelse:return arr_sort[~np.isnan(arr_sort)][0]

思路與之前一致，不過(guò)重構(gòu)為函數(shù)了。

具體代碼見(jiàn)項(xiàng)目：https://github.com/ChaoQiezi/read_fy3d_tshs

這是原始數(shù)據(jù)集：

這是目標(biāo)結(jié)果：