Sentinel-2A/B数据预处理之观测角度数据导出

最近分析了下哨兵二号数据，其角度数据让我烦不胜烦，写一个小程序将角度数据导出来为tif影像。
一、角度数据来源
 在类似下面S2B_MSIL1C_20180112T025049_N0206_R132_T50RPV_20180112T062213.SAFE\GRANULE\L1C_T50RPV_A004448_20180112T025254\MTD_TL.xml地址的文件中保存有哨兵数据的观测角度信息，但打开来是这样的：

看的人头皮发麻
二、tif角度数据生成
我用的是python程序
导入库：

from osgeo import gdal
import numpy as np
import xmltodict

提取经纬度坐标信息，使用任何一个样例文件位置：

path=r'S2B_MSIL1C_20180112T025049_N0206_R132_T50RPV_20180112T062213.SAFE\GRANULE\L1C_T50RPV_A004448_20180112T025254\IMG_DATA\T50RPV_20180112T025049_B02.jp2'#影像位置
root = gdal.Open(path)
Proj = root.GetProjection()
geotrans = root.GetGeoTransform()
im_geotrns = (geotrans[0],5000.0,geotrans[2],geotrans[3],geotrans[4],-5000.0)#调整分辨率为5000m
#坐标投影信息获取完毕

#导出数据为字典
xml_file=r'S2B_MSIL1C_20180112T025049_N0206_R132_T50RPV_20180112T062213.SAFE\GRANULE\L1C_T50RPV_A004448_20180112T025254\MTD_TL.xml'
    with open(xml_file, encoding="UTF-8") as xml_file1:
            parsed_data = xmltodict.parse(xml_file1.read())
            xml_file1.close()
vzn_m=int(float(parsed_data["n1:Level-1C_Tile_ID"]["n1:Geometric_Info"]["Tile_Angles"]["Mean_Viewing_Incidence_Angle_List"]["Mean_Viewing_Incidence_Angle"][1]["ZENITH_ANGLE"]["#text"]))
vaz_m=int(float(parsed_data["n1:Level-1C_Tile_ID"]["n1:Geometric_Info"]["Tile_Angles"]["Mean_Viewing_Incidence_Angle_List"]["Mean_Viewing_Incidence_Angle"][1]["AZIMUTH_ANGLE"]["#text"]))
szn_m=int(float(parsed_data["n1:Level-1C_Tile_ID"]["n1:Geometric_Info"]["Tile_Angles"]["Mean_Sun_Angle"]["ZENITH_ANGLE"]["#text"]))
saz_m=int(float(parsed_data["n1:Level-1C_Tile_ID"]["n1:Geometric_Info"]["Tile_Angles"]["Mean_Sun_Angle"]["AZIMUTH_ANGLE"]["#text"]))
#导出4个角度数据的均值(其路径位置都一样)
#导出4个单独角数据
"""
注意
除太阳天顶角与方位角，其他各波段观测角都是分别保存且分为7个部分
"""
    szn1=parsed_data["n1:Level-1C_Tile_ID"]["n1:Geometric_Info"]["Tile_Angles"]["Sun_Angles_Grid"]["Zenith"]["Values_List"]["VALUES"]
    saz1=parsed_data["n1:Level-1C_Tile_ID"]["n1:Geometric_Info"]["Tile_Angles"]["Sun_Angles_Grid"]["Azimuth"]["Values_List"]["VALUES"]
    height=len(szn1)
    width=len(szn1[0].split(" "))
    szn=np.zeros((height,width),dtype=np.uint8)
    saz=np.zeros((height,width),dtype=np.uint8)
    
    for i in range(height):
        for j in range(width):
            szn[i][j]=int(float(szn1[i].split(" ")[j]))
            saz[i][j]=int(float(saz1[i].split(" ")[j]))#太阳角提取完毕
    vzn1=parsed_data["n1:Level-1C_Tile_ID"]["n1:Geometric_Info"]["Tile_Angles"]["Viewing_Incidence_Angles_Grids"][28]["Zenith"]["Values_List"]["VALUES"]
    vaz1=parsed_data["n1:Level-1C_Tile_ID"]["n1:Geometric_Info"]["Tile_Angles"]["Viewing_Incidence_Angles_Grids"][28]["Azimuth"]["Values_List"]["VALUES"]

    vzn=np.zeros((height,width),dtype=np.uint8)
    vaz=np.zeros((height,width),dtype=np.int16)
    
    for d in range(28,35):
        vzn1=parsed_data["n1:Level-1C_Tile_ID"]["n1:Geometric_Info"]["Tile_Angles"]["Viewing_Incidence_Angles_Grids"][d]["Zenith"]["Values_List"]["VALUES"]
        vaz1=parsed_data["n1:Level-1C_Tile_ID"]["n1:Geometric_Info"]["Tile_Angles"]["Viewing_Incidence_Angles_Grids"][d]["Azimuth"]["Values_List"]["VALUES"]
        for i in range(height):
            for j in range(width):
                if vzn[i][j]==0 or vaz[i][j]==0:#注意7个部分直接相加存在部分区域重复相加的问题
                    vzn[i][j]=vzn[i][j]+int(float(vzn1[i].replace("NaN","0").split(" ")[j]))
                    vaz[i][j]=vaz[i][j]+int(float(vaz1[i].replace("NaN","0").split(" ")[j]))#观测角提取完毕
                    #zn天顶角 az方位角
    write_img(angle_name, Proj, im_geotrans, szn, saz, vzn,vaz)#保存文件到tif
def write_img(filename, proj, im_geotrans, da1, da2, da3,da4):
    driver = gdal.GetDriverByName("GTiff")  # 数据类型必须有，因为要计算需要多大内存空间
    dataset = driver.Create(filename, im_width, im_height, 4, gdal.GDT_Byte)
    dataset.SetGeoTransform(im_geotrans)  # 写入仿射变换参数
    dataset.SetProjection(im_proj)  # 写入投影
    dataset.GetRasterBand(1).WriteArray(da1)  # szn
    dataset.GetRasterBand(2).WriteArray(da2)  # saz
    dataset.GetRasterBand(3).WriteArray(da3)  # vzn
	dataset.GetRasterBand(4).WriteArray(da4)  # vaz
    del dataset,driver