Himawari-8 HSD數據解析及輻射定標

Himawari-8 HSD數據解析及輻射定標

1. HSD數據結構介紹

## Full-disk
 HS_H08_YYYYMMDD_hhmm_Bbb_FLDK_Rjj_Skkll.DAT

 where YYYY: 4-digit year of observation start time (timeline);
       MM: 2-digit month of timeline;
       DD: 2-digit day of timeline;
       hh: 2-digit hour of timeline;
       mm: 2-gidit minutes of timeline;
       bb: 2-digit band number (varies from "01" to "16");
       jj: spatial resolution ("05": 0.5km, "10": 1.0km, "20": 2.0km);
       kk: segment number (varies from "01" to "10"); and
       ll: total number of segments (fixed to "10").

 example: HS_H08_20150728_2200_B01_FLDK_R10_S0110.DAT

​ 具體命名格式如上，其中需要注意的是kk代表segment number，一開始筆者不明其意，下載了單景分辨率1km的全圓盤影像，發現行列號為(1100*11000)，於是很納悶為啥圓盤是這么個長方形。思索了很長時間沒有結果，於是休息了一會，今天突然發現了原因：segement number編號為01-10，其實是將(11000 * 11000)的圓盤圖像分割成了10個小圖像!!!

​ 於是立刻編寫代碼驗證想法，一經測試果真如此，下圖為數據合成后的結果。

2. DAT文件的讀取

​ 了解了Himawari-8 hsd數據的格式后，就可以開始讀取每個DAT文件了。下面這個鏈接給出了DAT文件的數據格式及其python讀取方法，其中將DAT文件按一定format讀取，轉換成了可以打開的txt。

​ DAT數據由12個block組成，每個block有自己的組織形式，如下圖。

​ https://blog.csdn.net/qq_33339770/article/details/102708243?utm_medium=distribute.pc_relevant.none-task-blog-searchFromBaidu-4.control&depth_1-utm_source=distribute.pc_relevant.none-task-blog-searchFromBaidu-4.control

讀取方法如下：1. 對於數據量小的數據可以用f.read()逐字讀取

					   2. 對於數據量大的數據用`np.fromfile()`配合**format**讀取
					   3. 其實不用`f.seek()`也可以用`re`匹配
					   4. `pandas.read_table()`也可以讀取.dat
					   5. **輻射定標**：可在DAT文件中讀取定標系數，進行輻射定標；此外**Albedo/Tbb的反射率定標方法有所不同**。

def read_hsd(inputfile):
    '''
    @desc: 打開單個segment圖幅，讀取.DAT文件中的信息(data/gains/bias...)
    @method: np.fromfile()
    '''
    resolution = int(inputfile[-12:-10])
    if resolution == 10:
        rows = 1100
        cols = 11000
    elif resolution == 20:
        rows = 550
        cols = 5500
    else:
        rows = 2200
        cols = 22000
    # 打開文件
    f = open(inputfile, 'rb')
    # 獲取Observation start time
    f.seek(46)
    imgtime = struct.unpack('d', f.read(8))[0]
    # 獲取Total header length
    f.seek(70)
    header_length = struct.unpack('i', f.read(4))[0]
    # 獲取影像
    formation = [('header', 'S1', header_length), ('pixel', 'i2', rows*cols)]
    imgdata = np.fromfile(inputfile, dtype=formation)['pixel'].reshape(rows, cols)

    '''
    @desc: 下面是輻射定標部分，如果不需要可省略
    '''
    # if resolution != 20:
    #     # 重采樣至550行，5500列，可省略
    #     imgdata = imgdata.reshape(550, int(20/resolution), 5500, int(20/resolution)).mean(axis=(1, 3))
    
    # 獲取Sun's position
    f.seek(510)
    sun_pos1 = struct.unpack('d', f.read(8))[0]
    # print(sun_pos1)
    f.seek(510+8)
    sun_pos2 = struct.unpack('d', f.read(8))[0]
    # print(sun_pos2)
    f.seek(510+8+8)
    sun_pos3 = struct.unpack('d', f.read(8))[0]
    # print(sun_pos3)
    # 獲取Band number
    f.seek(601)
    band_num = int.from_bytes(f.read(2), byteorder='little', signed=False)
    # 獲取Gain
    f.seek(617)
    Gain = struct.unpack('d', f.read(8))[0]
    # print(Gain)
    # 獲取Offset
    f.seek(625)
    Offset = struct.unpack('d', f.read(8))[0]
    # print(Offset)
    # 計算radiance
    radiance = imgdata*Gain + Offset

    # print(radiance[0][1580])
    # 對前6波段定標成反射率
    if band_num <= 6:
        # 獲取radiance to albedo
        f.seek(633)
        cc = struct.unpack('d', f.read(8))[0]
        f.close
        # 計算反射率
        albedo = radiance * cc
        outdata = albedo
    # 后面的波段定標成計算亮溫
    else:
        # 獲取Central wave length
        f.seek(603)
        wv = struct.unpack('d', f.read(8))[0]
        # 獲取radiance to brightness temperature(c0)
        f.seek(633)
        c0 = struct.unpack('d', f.read(8))[0]
        # 獲取radiance to brightness temperature(c1)
        f.seek(641)
        c1 = struct.unpack('d', f.read(8))[0]
        # 獲取radiance to brightness temperature(c2)
        f.seek(649)
        c2 = struct.unpack('d', f.read(8))[0]
        # 獲取Speed of light(c)
        f.seek(681)
        c = struct.unpack('d', f.read(8))[0]
        # 獲取Planck constant(h)
        f.seek(689)
        h = struct.unpack('d', f.read(8))[0]
        # 獲取Boltzmann constant(k)
        f.seek(697)
        k = struct.unpack('d', f.read(8))[0]
        f.close
        # 計算亮溫
        wv = wv * 1e-6
        rad = radiance * 1e6
        Te = h*c/k/wv/(np.log(2*h*c*c/((wv**5)*rad)+1))
        BT = c0 + c1 * Te + c2 * Te * Te
        # 返回值
        outdata = BT

    # 返回：albedo / BT, 時間, 太陽坐標
    sunpos = [sun_pos1, sun_pos2, sun_pos3]
    out = {'pixels': list(outdata.flatten()), 'time': imgtime, 'sun_pos': sunpos}
    return out

3. 經緯度和行列號互轉

​ 與FY-4A類似，Full-disk數據在幾何校正的時候都需要進行經緯度與行列號的轉換，從而把NOM投影轉換成等經緯度投影。兩者計算的唯一差別在於衛星星下點經緯度、行偏移、列偏移的不同(COFF/LFAC...)，其他公式都一樣。

下面給出經緯度與行列號的轉換代碼：

def calParameter(resolution):
    '''
    @desc: 計算Himawari-8 COFF/CFAC/LOFF/LFAC
    @resolution: 10:1km/20:2km
    '''
    if resolution == 20:
        row = 550
        col = 5500
        COFF = LOFF = 2750.5
        CFAC = LFAC = 20466275
    elif resolution == 10:
        row = 1100
        col = 11000
        COFF = LOFF = 5500.5
        CFAC = LFAC = 40932549
    return COFF, LOFF, CFAC, LFAC

def latlon2linecolumn(lat, lon, resolution):
    """
    經緯度轉行列
    (lat, lon) → (line, column)
    resolution：文件名中的分辨率
    line, column
    """
    COFF, LOFF, CFAC, LFAC = calParameter(resolution)
    ea = 6378.137  # 地球的半長軸[km]
    eb = 6356.7523  # 地球的短半軸[km]
    h = 42164  # 地心到衛星質心的距離[km]
    λD = np.deg2rad(140.7)  # 衛星星下點所在經度
    # Step1.檢查地理經緯度
    # Step2.將地理經緯度的角度表示轉化為弧度表示
    lat = np.deg2rad(lat)
    lon = np.deg2rad(lon)
    # Step3.將地理經緯度轉化成地心經緯度
    eb2_ea2 = eb ** 2 / ea ** 2
    λe = lon
    φe = np.arctan(eb2_ea2 * np.tan(lat))
    # Step4.求Re
    cosφe = np.cos(φe)
    re = eb / np.sqrt(1 - (1 - eb2_ea2) * cosφe ** 2)
    # Step5.求r1,r2,r3
    λe_λD = λe - λD
    r1 = h - re * cosφe * np.cos(λe_λD)
    r2 = -re * cosφe * np.sin(λe_λD)
    r3 = re * np.sin(φe)
    # Step6.求rn,x,y
    rn = np.sqrt(r1 ** 2 + r2 ** 2 + r3 ** 2)
    x = np.rad2deg(np.arctan(-r2 / r1))
    y = np.rad2deg(np.arcsin(-r3 / rn))
    # Step7.求c,l
    column = COFF + x * 2 ** -16 * CFAC
    line = LOFF + y * 2 ** -16 * LFAC
    return np.rint(line).astype(np.uint16), np.rint(column).astype(np.uint16)

4. 幾何校正

​ 最后就是進行幾何校正了，需要我們設定目標區的經緯度范圍（如Lon: 90-120 Lat: 30-50），然后根據pixelSize計算出geotrans、col/row等信息。最后根據經緯度反算出行列號，進而求得行列號對應的像元的值，寫入數組，這樣就完成了幾何校正。

1. 切記行列號順序不可弄錯，numpy先行后列，gdal先列后行。

   

def hsd2tif(arr, lonMin, lonMax, latMin, latMax, pixelSize, file_out):
    '''
    @desc: 幾何校正，輸出tif
    '''
    col = np.ceil(lonMax-lonMin)/pixelSize
    row = np.ceil(latMax-latMin)/pixelSize
    col = int(col)
    row = int(row)
    xnew = np.linspace(latMin, latMax, row)
    ynew = np.linspace(lonMin, lonMax, col)
    xnew, ynew = np.meshgrid(xnew, ynew)
    dataGrid = np.zeros((row, col))  # numpy數組先行后列!!!
    
    index = {}
    for i in tqdm(range(row)):
        for j in tqdm(range(col)):
            lat = xnew[i][j]
            lon = ynew[i][j]
            h8Row = 0
            h8Col = 0
            if index.get((lat, lon)) == None:
                h8Row, h8Col = latlon2linecolumn(lat, lon, 10)
                index[(lat, lon)] = h8Row, h8Col
            else:
                h8Row, h8Col = index.get((lat, lon))
                
            if (0<=h8Row<11000) and (0<=h8Col<11000):
                dataGrid[i][j] = arr[h8Row, h8Col]
                print("Worked!" + str(dataGrid[i][j]))
            else:
                print("該坐標(%s, %s)不在影像中"%(lon, lat))

    geotrans = (lonMin, pixelSize, 0, latMax, 0, -pixelSize)
    srs = osr.SpatialReference()
    srs.ImportFromEPSG(4326)
    proj = srs.ExportToWkt()
    raster2tif(dataGrid, geotrans, proj, file_out)

5. 后記

​ 到這里我們的HSD數據解析和輻射定標就完成了。

1. 后續還可以調用6S模型進行大氣校正，這也需要更多的地理信息如(SALA/SAZA/SOLA/SOZA)，https://blog.csdn.net/qq_33339770/article/details/103023725 這里給出了一些參數的計算方法。

2. 根據經緯度反算行列號，計算速度偏慢。

   如果需要快速出圖的話，建議使用已經計算好的經緯度柵格圖進行GLT校正，速度會快很多。（ENVI有GLT校正模塊，gdal應該也有，后續可以建立好經緯度柵格圖然后調用gdal來批量幾何校正）

​