今天與第9天一樣，處理GeoTiff檔案
只不過是換成高解析可見光反照率資料-東亞

相關資料的下載如第9天，只不過xml檔名變更為O-B0055-001.xml

以下就是程式範例囉

import zipfile, requests
import xmltodict as xdict

with open("O-B0055-001.xml", "r", encoding="utf-8") as fn:
    satdict = xdict.parse(fn.read())
    
url = satdict["cwbopendata"]["dataset"]["resource"]["uri"] #取得下載連結
print(url)
r = requests.get(url, stream=True) #下載檔案為連續檔案於記憶體

#將下載的檔案變為壓縮檔
with open("test.zip", 'wb') as fd:
    for chunk in r.iter_content(chunk_size=128):
        fd.write(chunk)

#解壓縮到該資料夾下的test2資料夾
extractdir = "test2"
with zipfile.ZipFile("test.zip", 'r') as zip_ref:
    listOfFileNames = zip_ref.namelist()
    for fileName in listOfFileNames:
        if fileName.endswith('.gtif'):
            zip_ref.extract(fileName, extractdir)
print(extractdir +"/" +fileName)

上面的程式執行完成後，應該會得到
https://cwbopendata.s3.ap-northeast-1.amazonaws.com/MSC/O-B0055-001.gtif.zip
test2/data4/hsd/PG261/2022-09-24_1000.B03.LCC.gtif

那現在就是要讀取2022-09-24_1000.B03.LCC.gtif
這邊介紹使用gdal來解，安裝方式如下
conda install -c conda-forge gdal
或
pip3 install gdal

安裝完成後就開始解吧

gtifpath = extractdir +"/" +fileName
sattif = gdal.Open(gtifpath,gdal.GA_ReadOnly)
data  = sattif.GetRasterBand(1) #取值
data2arr = data.ReadAsArray() #僵值轉換為array
albdo = data2arr/1000. #說明文件有提到要將數值除以1000

上面的程式就將亮度溫度讀取完成

下面提供簡易的視覺化(沒有加投影資訊)，這邊因違反照率愈強，即是雲物體，故將cmap變成反轉

fig , axs = plt.subplots(1,2,figsize=(8,8))
axs[0].imshow(plt.imread("LCC_VIS_TRGB_2750-2022-09-24-10-10.jpg"))
ishow = axs[1].imshow(albdo,cmap="Blues_r",extent=(lonst,lonend,latst,latend))
plt.colorbar(ishow,fraction=0.05)
plt.title(obsTime)

左邊的圖為氣象局觀測的真實色雲圖，拿來驗證用的
右邊的圖是透過開放資料提供的反照率繪製而成，看起來我們是有將資料讀對的