今天也是要介紹衛星資料，但是用hdf格式儲存的。至於為什麼氣象局開放資料沒有統一用GeoTiff格式提供衛星檔案，這我就不得而知了。

hdf全名是Hierarchical Data Format，中文稱為層級資料格式，這是一種專門在組織及儲存大量資料的一種格式。
詳細的介紹可以點這。

我們依舊照著步驟，先下載xml檔案，

利用瀏覽器將xml檔案開啟，如下圖所示

畫黃色螢光筆之處，是真正的資料下載連結。

所以第一步還是要寫個下載程式檔案。

import xmltodict as xdict
import  requests

with open("O-A0043-001.xml", "r", encoding="utf-8") as fn:
    satdict = xdict.parse(fn.read())
url = satdict["cwbopendata"]["dataset"]["resource"]["uri"]

fnout = url.split("/")[-1]

r = requests.get(url, stream=True)
with open(fnout, 'wb') as fd:
    for chunk in r.iter_content(chunk_size=128):
        fd.write(chunk)

執行完成上述指令後，應該就會有O-A0043-001.hdf檔案下載完成。

這一份hdf檔案其實是hdf4格式儲存的，這邊介紹使用pyhdf套件處理hdf4格式資料。
先下載套件，指令參考如下
conda install -c conda-forge pyhdf
或
pip 3 install pyhdf

開始讀取吧，程式碼如下

from pyhdf.SD import SD, SDC 
sathdf  = SD("O-A0043-001.hdf",SDC.READ)
print(dir(sathdf))

print(sathdf.datasets()) #看資料檔頭

讀取資料之後，先來看這筆資料有什麼內容，使用print(dir(sathdf))應該會得到如下的訊息，
['class', 'del', 'delattr', 'dict', 'dir', 'doc', 'eq', 'format', 'ge', 'getattr', 'getattribute', 'gt', 'hash', 'init', 'init_subclass', 'le', 'lt', 'module', 'ne', 'new', 'reduce', 'reduce_ex', 'repr', 'setattr', 'sizeof', 'str', 'subclasshook', 'weakref', '_id', 'attr', 'attributes', 'create', 'datasets', 'end', 'info', 'nametoindex', 'reftoindex', 'select', 'setfillmode']

這邊常用的方法為attributes, select 跟datasets。
attributes功用：是看一些基本的設定檔，譬如這份文件的資料格、創建時間、觀測時間及投影資訊等。
datasets功用：是看資料的檔頭名稱，可以想像成python字典物件的key，主要是給select方法使用，譬如使用datasets方法後，有一份資料為lat，那這個lat就可以成為select方法的參數。
select功用：在datasets功用以說明。

print(sathdf.datasets())會得到如下資訊

可以得知這是一份提供3種資料的層級資料(好像有點怪，哈)，知道檔頭之後，就要把資料抓出來囉。

以下無抓出資料的程式碼，以抓出頻道03反射率為例

"""
由sathdf.datasets()的資訊得知頻道03反射率為refl_0_65um_nom_atmos_corr。
但hdf4儲存資料有的規則為每一個資料有自己對應的偏差(add_offset)及刻度因子(scale_factor)，必須要還原。
使用sathdf.select("變數").attributes()，就可以得到相關資訊。
會這樣做的主要目的是降低資料儲存空間。
由於每一個資料都要這樣做，故建議寫個function。這邊就先徒法煉鋼一下。
"""

print(sathdf.select("refl_0_65um_nom_atmos_corr").attributes())
scale_factor = sathdf.select("refl_0_65um_nom_atmos_corr").attributes()["scale_factor"]
add_offset = sathdf.select("refl_0_65um_nom_atmos_corr").attributes()["add_offset"]
#還原成正確的值
hdf_true = sathdf.select("refl_0_65um_nom_atmos_corr").get()*scale_factor + add_offset

print(sathdf.select("refl_0_65um_nom_atmos_corr").attributes())會得到以下資訊，
{'SCALED': 1, 'add_offset': 59.0, 'scale_factor': 0.0018616290763020515, 'units': '%', 'standard_name': 'toa_bidirectional_pseudo_reflectance_ 0_65um_atmos_corr', 'long_name': 'observed pseudo-reflectance, atmospherically corrected, at the nominal wavelength of 0_65um', 'coordinates': 'longitude latitude', 'valid_range': [-32767, 32767], 'actual_range': [-2.0, 120.0], 'actual_missing': -999.0}
可以看到其實有蠻多訊息的，有add_offset, scale_factor，尤其有個actual_range，這就是在提醒經過還原後正確的值會介於的範圍。

解完資料之後，一樣視覺化囉，由於頻道03反射率通常做為可見功的紅色波段，所以我就用紅色表示。

import matplotlib.pyplot as plt
import cartopy.crs as crs
import cartopy.io.shapereader as shpreader
from cartopy.feature import ShapelyFeature
import matplotlib.colors as mcolors

fig = plt.figure(figsize=(12,8))
axs = plt.axes(projection=crs.PlateCarree())
tw_shp = ShapelyFeature(shpreader.Reader("D:\ithome\/tw_shp\COUNTY_MOI_1090820.shp").geometries(), crs.PlateCarree(),facecolor="none",edgecolor='k')
axs.gridlines(draw_labels=True)
ims  = axs.imshow(hdf_true/255,cmap="Reds",extent=(100, 140,9,40),transform=crs.PlateCarree())
axs.add_feature(tw_shp)
plt.colorbar(ims)
filename = sathdf.attributes()["FILENAME"].split("_")
obstime = filename[2] + " " +filename[3] + " UTC"
plt.title(obstime)