Himawari-8HSD数据解析及辐射定标
Himawari-8 HSD数据解析及辐射定标
1. HSD数据结构介绍
## Full-disk
HS_H08_YYYYMMDD_hhmm_Bbb_FLDK_Rjj_Skkll.DAT
where YYYY: 4-digit year of obrvation 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: gment number (varies from "01" to "10"); and
ll: total number of gments (fixed to "10").
example: HS_H08_20150728_2200_B01_FLDK_R10_S0110.DAT
具体命名格式如上，其中需要注意的是kk代表gment number，⼀开始笔者不明其意，下载了单景分辨率1km的全圆盘影像，发现⾏列号为(1100*11000)，于是很纳闷为啥圆盘是这么个长⽅形。思索了很长时间没有结果，于是休息了⼀会，今天突然发现了原因：gement number编号为01-10，其实是将(11000 * 11000)的圆盘图像分割成了10个⼩图像!!!
于是⽴刻编写代码验证想法，⼀经测试果真如此，下图为数据合成后的结果。
2. DAT⽂件的读取
了解了Himawari-8 hsd数据的格式后，就可以开始读取每个DAT⽂件了。下⾯这个链接给出了DAT⽂件的数据格式及其python读取⽅法，其中将DAT⽂件按⼀定format读取，转换成了可以打开的txt。
DAT数据由12个block组成，每个block有⾃⼰的组织形式，如下图。
读取⽅法如下：1. 对于数据量⼩的数据可以⽤f.read()逐字读取
2. 对于数据量⼤的数据⽤`np.fromfile()`配合**format**读取
3. 其实不⽤`f.ek()`也可以⽤`re`匹配
4. `ad_table()`也可以读取.dat
5. **辐射定标**：可在DAT⽂件中读取定标系数，进⾏辐射定标；此外**Albedo/Tbb的反射率定标⽅法有所不同**。def read_hsd(inputfile):
'''
@desc: 打开单个gment图幅，读取.DAT⽂件中的信息(data/)
@method: np.fromfile()
'''
resolution = int(inputfile[-12:-10])
if resolution == 10:
rows = 1100
cols = 11000
elif resolution == 20:
俞敏洪的励志故事rows = 550
cols = 5500
el:
rows = 2200
cols = 22000
# 打开⽂件
f = open(inputfile, 'rb')
# 获取Obrvation start time
f.ek(46)
imgtime = struct.unpack('d', f.read(8))[0]
# 获取Total header length
f.ek(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 = shape(550, int(20/resolution), 5500, int(20/resolution)).mean(axis=(1, 3))
# 获取Sun's position
f.ek(510)
sun_pos1 = struct.unpack('d', f.read(8))[0]
# print(sun_pos1)
f.ek(510+8)
sun_pos2 = struct.unpack('d', f.read(8))[0]
# print(sun_pos2)
f.ek(510+8+8)
sun_pos3 = struct.unpack('d', f.read(8))[0]
# print(sun_pos3)
# 获取Band number
f.ek(601)
band_num = int.from_ad(2), byteorder='little', signed=Fal)
# 获取Gain
f.ek(617)
Gain = struct.unpack('d', f.read(8))[0]
# print(Gain)
# 获取Offt
f.ek(625)
Offt = struct.unpack('d', f.read(8))[0]
# print(Offt)
# 计算radiance
radiance = imgdata*Gain + Offt
皮蛋豆腐做法
# print(radiance[0][1580])
# 对前6波段定标成反射率
if band_num <= 6:
# 获取radiance to albedo
f.ek(633)
cc = struct.unpack('d', f.read(8))[0]
f.clo
# 计算反射率
albedo = radiance * cc
outdata = albedo
# 后⾯的波段定标成计算亮温
el:
# 获取Central wave length
f.ek(603)
wv = struct.unpack('d', f.read(8))[0]
# 获取radiance to brightness temperature(c0)
f.ek(633)
c0 = struct.unpack('d', f.read(8))[0]
# 获取radiance to brightness temperature(c1)
f.ek(641)
c1 = struct.unpack('d', f.read(8))[0]
# 获取radiance to brightness temperature(c2)
f.ek(649)
c2 = struct.unpack('d', f.read(8))[0]
# 获取Speed of light(c)
f.ek(681)
c = struct.unpack('d', f.read(8))[0]
# 获取Planck constant(h)
f.ek(689)
h = struct.unpack('d', f.read(8))[0]
# 获取Boltzmann constant(k)
f.ek(697)幼儿绘画入门
k = struct.unpack('d', f.read(8))[0]
f.clo
# 计算亮温
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投影转换成等经纬度投影。两者计算的唯⼀差别在于卫星星下点经纬度、⾏偏移、列偏移的不同()，其他公式都⼀样。
下⾯给出经纬度与⾏列号的转换代码：
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
((lat, lon)) == None:
h8Row, h8Col = latlon2linecolumn(lat, lon, 10)
端午节包粽子index[(lat, lon)] = h8Row, h8Col
el:
h8Row, h8Col = ((lat, lon))
文档背景图
if (0<=h8Row<11000) and (0<=h8Col<11000):
dataGrid[i][j] = arr[h8Row, h8Col]
print("Worked!" + str(dataGrid[i][j]))
小学英语怎么学el:
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)，这⾥给出了⼀些参数的计算⽅法。
2. 根据经纬度反算⾏列号，计算速度偏慢。
如果需要快速出图的话，建议使⽤已经计算好的经纬度栅格图进⾏GLT校正，速度会快很多。（ENVI有GLT校正模块，gdal应该也有，后续可以建⽴好经纬度栅格图然后调⽤gdal来批量⼏何校正）
