关于生成Sentinel-2 L1C逐像素的太阳天顶角、太阳方位角、传感器天顶角、传感器方位角tif

Sentinel-2 L1C的太阳角度是以XML文件给的，在I:\S2B_MSIL1C_20201106T134209_N0500_R124_T22KCA_20230410T094912\S2B_MSIL1C_20201106T134209_N0500_R124_T22KCA_20230410T094912.SAFE\GRANULE\L1C_T22KCA_A019169_20201106T134212\MTD_TL.xml这个路径下面，但是是给的5km一个角度值，但是实际项目需要逐像素的角度文件，所以就需要对这个进行处理
主要利用python解析xml，用arcpy进行插值计算

import xml.etree.ElementTree as ET
import numpy as np
import math
from osgeo import ogr,osr
from osgeo import gdal
import os
import arcpy

from math import sqrt, cos, sin, tan, pi, asin, acos, atan, atan2
class OutOfRangeError(Exception):
	pass
a = 6378137.0                   # WGS 84 semi-major axis in meters
b = 6356752.314                 # WGS 84 semi-minor axis in meters
ecc = 1.0 - b / a * b / a       # WGS 84 ellipsoid eccentricity
todeg = 180.0 / pi		# Converts radians to degrees
K0 = 0.9996

E = 0.00669438
E2 = E * E
E3 = E2 * E
E_P2 = E / (1.0 - E)

SQRT_E = math.sqrt(1 - E)
_E = (1 - SQRT_E) / (1 + SQRT_E)
_E2 = _E * _E
_E3 = _E2 * _E
_E4 = _E3 * _E
_E5 = _E4 * _E

M1 = (1 - E / 4 - 3 * E2 / 64 - 5 * E3 / 256)
M2 = (3 * E / 8 + 3 * E2 / 32 + 45 * E3 / 1024)
M3 = (15 * E2 / 256 + 45 * E3 / 1024)
M4 = (35 * E3 / 3072)

P2 = (3. / 2 * _E - 27. / 32 * _E3 + 269. / 512 * _E5)
P3 = (21. / 16 * _E2 - 55. / 32 * _E4)
P4 = (151. / 96 * _E3 - 417. / 128 * _E5)
P5 = (1097. / 512 * _E4)

R = 6378137

ZONE_LETTERS = "CDEFGHJKLMNPQRSTUVWXX"
def 创建矢量点(AngleObs,shpname,raster_file):
    raster_ds = gdal.Open(raster_file)
    raster_srs = osr.SpatialReference()
    raster_srs.ImportFromWkt(raster_ds.GetProjection())
    driver = ogr.GetDriverByName("ESRI Shapefile")
    #output_layer = output_ds.CreateLayer("new_layer_name", srs=raster_srs, geom_type=ogr.wkbPolygon)
    data_source = driver.CreateDataSource(shpname)
    layer = data_source.CreateLayer('new_points', srs=raster_srs, geom_type=ogr.wkbPoint)
    f_ = ogr.FieldDefn('sunZenith', ogr.OFTReal)
    f_.SetWidth(32)
    f_.SetPrecision(16)
    layer.CreateField(f_)
    f_ = ogr.FieldDefn('sunAzimuth', ogr.OFTReal)
    f_.SetWidth(32)
    f_.SetPrecision(16)
    layer.CreateField(f_)
    for angle in AngleObs['obs']:
        point = ogr.Geometry(ogr.wkbPoint)
        point.AddPoint(angle[0], angle[1])
        feature = ogr.Feature(layer.GetLayerDefn())
        feature.SetField('sunZenith', angle[2])
        feature.SetField('sunAzimuth', angle[3])
        feature.SetGeometry(point)
        layer.CreateFeature(feature)

def from_latlon(latitude, longitude, force_zone_number=None):
    if not -80.0 <= latitude <= 84.0:
        raise OutOfRangeError('latitude out of range (must be between 80 deg S and 84 deg N)')
    if not -180.0 <= longitude <= 180.0:
        raise OutOfRangeError('longitude out of range (must be between 180 deg W and 180 deg E)')
    lat_rad = math.radians(latitude)
    lat_sin = math.sin(lat_rad)
    lat_cos = math.cos(lat_rad)

    lat_tan = lat_sin / lat_cos
    lat_tan2 = lat_tan * lat_tan
    lat_tan4 = lat_tan2 * lat_tan2

    if force_zone_number is None:
        zone_number = latlon_to_zone_number(latitude, longitude)
    else:
        zone_number = force_zone_number

    zone_letter = latitude_to_zone_letter(latitude)

    lon_rad = math.radians(longitude)
    central_lon = zone_number_to_central_longitude(zone_number)
    central_lon_rad = math.radians(central_lon)

    n = R / math.sqrt(1 - E * lat_sin**2)
    c = E_P2 * lat_cos**2

    a = lat_cos * (lon_rad - central_lon_rad)
    a2 = a * a
    a3 = a2 * a
    a4 = a3 * a
    a5 = a4 * a
    a6 = a5 * a

    m = R * (M1 * lat_rad -
             M2 * math.sin(2 * lat_rad) +
             M3 * math.sin(4 * lat_rad) -
             M4 * math.sin(6 * lat_rad))

    easting = K0 * n * (a +
                        a3 / 6 * (1 - lat_tan2 + c) +
                        a5 / 120 * (5 - 18 * lat_tan2 + lat_tan4 + 72 * c - 58 * E_P2)) + 500000

    northing = K0 * (m + n * lat_tan * (a2 / 2 +
                                        a4 / 24 * (5 - lat_tan2 + 9 * c + 4 * c**2) +
                                        a6 / 720 * (61 - 58 * lat_tan2 + lat_tan4 + 600 * c - 330 * E_P2)))

    if latitude < 0:
        northing += 10000000

    return easting, northing, zone_number, zone_letter


def latitude_to_zone_letter(latitude):
    if -80 <= latitude <= 84:
        return ZONE_LETTERS[int(latitude + 80) >> 3]
    else:
        return None


def latlon_to_zone_number(latitude, longitude):
    if 56 <= latitude < 64 and 3 <= longitude < 12:
        return 32

    if 72 <= latitude <= 84 and longitude >= 0:
        if longitude <= 9:
            return 31
        elif longitude <= 21:
            return 33
        elif longitude <= 33:
            return 35
        elif longitude <= 42:
            return 37

    return int((longitude + 180) / 6) + 1


def zone_number_to_central_longitude(zone_number):
    return (zone_number - 1) * 6 - 180 + 3
# Inverse (X/Y to lat/long) UTM projection
def utm_inv( Zone, X, Y, a=6378137.0, b=6356752.31414 ):
        if Zone < 0 :
                FNorth = 10000000.0     # Southern hemisphere False Northing
        else:
                FNorth = 0.0	        # Northern hemisphere False Northing
        FEast = 500000.0                # UTM False Easting
        Scale = 0.9996                  # Scale at CM (UTM parameter)
        LatOrigin = 0.0                 # Latitude origin (UTM parameter)
        CMDeg = -177 + (abs(int(Zone))-1)*6
        CM = float(CMDeg)*pi/180.0      # Central meridian (based on zone)
        ecc = 1.0 - b/a*b/a
        ep = ecc/(1.0-ecc)
        M0 = a*((1.0-ecc*(0.25+ecc*(3.0/64.0+ecc*5.0/256.0)))*LatOrigin
               -ecc*(0.375+ecc*(3.0/32.0+ecc*45.0/1024.0))*sin(2.0*LatOrigin)
               +ecc*ecc*(15.0/256.0+ecc*45.0/1024.0)*sin(4.0*LatOrigin)
               -ecc*ecc*ecc*35.0/3072.0*sin(6.0*LatOrigin))
        M = M0+(Y-FNorth)/Scale
        Mu = M/(a*(1.0-ecc*(0.25+ecc*(3.0/64.0+ecc*5.0/256.0))))
        e1 = (1.0-sqrt(1-ecc))/(1.0+sqrt(1.0-ecc))
        Phi1 = Mu+(e1*(1.5-27.0/32.0*e1*e1)*sin(2.0*Mu)
                   +e1*e1*(21.0/16.0-55.0/32.0*e1*e1)*sin(4.0*Mu)
                   +151.0/96.0*e1*e1*e1*sin(6.0*Mu)
                   +1097.0/512.0*e1*e1*e1*e1*sin(8.0*Mu))  
        slat = sin(Phi1)
        clat = cos(Phi1)
        Rn1 = a/sqrt(1.0-ecc*slat*slat)
        T1 = slat*slat/clat/clat
        C1 = ep*clat*clat
        R1 = Rn1*(1.0-ecc)/(1.0-ecc*slat*slat)
        D = (X-FEast)/Rn1/Scale
        # Calculate Lat/Lon
        Lat = Phi1 - (Rn1*slat/clat/R1*(D*D/2.0
                        -(5.0+3.0*T1+10.0*C1-4.0*C1*C1-9.0*ep)*D*D*D*D/24.0
                        +(61.0+90.0*T1+298.0*C1+45.0*T1*T1-252.0*ep-3.0*C1*C1)*D*D*D*D*D*D/720.0))
        Lon = CM + (D-(1.0+2.0*T1+C1)*D*D*D/6.0+(5.0-2.0*C1+28.0*T1-3.0*C1*C1+8.0*ep+24.0*T1*T1)
                    *D*D*D*D*D/120.0)/clat
        
        return (Lat, Lon)

def LOSVec( Lat, Lon, Zen, Az ):
	LSRx = ( -sin(Lon), cos(Lon), 0.0 )
	LSRy = ( -sin(Lat)*cos(Lon), -sin(Lat)*sin(Lon), cos(Lat) )
	LSRz = ( cos(Lat)*cos(Lon), cos(Lat)*sin(Lon), sin(Lat) )
	LOS = ( sin(Zen)*sin(Az), sin(Zen)*cos(Az), cos(Zen) )
	Sat = ( LOS[0]*LSRx[0] + LOS[1]*LSRy[0] + LOS[2]*LSRz[0],
	        LOS[0]*LSRx[1] + LOS[1]*LSRy[1] + LOS[2]*LSRz[1],
	        LOS[0]*LSRx[2] + LOS[1]*LSRy[2] + LOS[2]*LSRz[2] )
	Rn = a / sqrt( 1.0 - ecc *sin(Lat)*sin(Lat))
	Gx = ( Rn*cos(Lat)*cos(Lon),
		Rn*cos(Lat)*sin(Lon),
		Rn*(1-ecc)*sin(Lat) )
	return ( Sat, Gx )
def get_angleobssun( XML_File ):
    tree = ET.parse(XML_File)
    root = tree.getroot()
    for child in root:
            if child.tag[-12:] == 'General_Info':
                geninfo = child
            if child.tag[-14:] == 'Geometric_Info':
                geoinfo = child

    for segment in geninfo:
        if segment.tag == 'TILE_ID':
            tile_id = segment.text.strip()

    for segment in geoinfo:
        if segment.tag == 'Tile_Geocoding':
            frame = segment
        if segment.tag == 'Tile_Angles':
                angles = segment

    for box in frame:
        if box.tag == 'HORIZONTAL_CS_NAME':
            czone = box.text.strip()[-3:]
            hemis = czone[-1:]
            zone = int(czone[:-1])
        if box.tag == 'Size' and box.attrib['resolution'] == '60':
            for field in box:
                if field.tag == 'NROWS':
                    nrows = int(field.text)
                if field.tag == 'NCOLS':
                    ncols = int(field.text)
        if box.tag == 'Geoposition' and box.attrib['resolution'] == '60':
            for field in box:
                if field.tag == 'ULX':
                    ulx = float(field.text)
                if field.tag == 'ULY':
                    uly = float(field.text)
    if hemis == 'S':
        lzone = -zone
    else:
        lzone = zone
    Angle = { 'zone' : zone, 'hemis' : hemis, 'nrows' : nrows, 'ncols' : ncols, 'ul_x' : ulx, 'ul_y' : uly, 'obs' : [] }
    sunangle = { 'zone' : zone, 'hemis' : hemis, 'nrows' : nrows, 'ncols' : ncols, 'ul_x' : ulx, 'ul_y' : uly, 'obs' : [] }
    for angle in angles:
        #print(angle.tag)
        if angle.tag == 'Sun_Angles_Grid':	
            for bset in angle:
                if bset.tag == 'Zenith':
                    zenith = bset
                if bset.tag == 'Azimuth':
                    azimuth = bset	
            for field in zenith:
                if field.tag == 'COL_STEP':
                    col_step = int(field.text)
                if field.tag == 'ROW_STEP':
                    row_step = int(field.text)
                if field.tag == 'Values_List':
                    zvallist = field
            for field in azimuth:
                if field.tag == 'Values_List':
                    avallist = field
            for rindex in range(len(zvallist)):
                zvalrow = zvallist[rindex]
                avalrow = avallist[rindex]
                zvalues = zvalrow.text.split(' ')
                avalues = avalrow.text.split(' ')
                values = list(zip( zvalues, avalues ))
                ycoord = uly - rindex * row_step
                #for cindex in values:
                for cindex in range(len(values)):
                    xcoord = ulx + cindex * col_step
                    (lat, lon) = utm_inv( lzone, xcoord, ycoord )
                    if ( values[cindex][0] != 'NaN' and values[cindex][1] != 'NaN' ):
                        zen = float( values[cindex][0] ) / todeg    #转换为弧度制
                        az = float( values[cindex][1] ) / todeg
                        zen1 = float( values[cindex][0] ) 
                        az1 = float( values[cindex][1] ) 
                        ( Sat, Gx ) = LOSVec( lat, lon, zen, az )
                        observe = [xcoord, ycoord, Sat, Gx ]
                        Angle['obs'].append( observe )
                        sunobserve = [xcoord, ycoord, zen1, az1 ]
                        sunangle['obs'].append( sunobserve )
    return (tile_id, Angle, sunangle)

XML_File='I:\\HLSwGF6\\bx11_6\S2B_MSIL2A_20201106T134209_N0214_R124_T22KCA_20201106T160551\\S2B_MSIL2A_20201106T134209_N0214_R124_T22KCA_20201106T160551.SAFE\\GRANULE\\L2A_T22KCA_A019169_20201106T134212\\MTD_TL.xml'

# Sudipta spatial subset setting
sul_lat = sul_lon = slr_lat = slr_lon = None
(Tile_ID, AngleObs,sunangles) = get_angleobssun( XML_File )
Tile_Base = Tile_ID.split('.')
# 应该是使gdal库支持中文编码、中文字段、中文路径
gdal.SetConfigOption("GDAL_ARRAY_OPEN_BY_FILENAME","TRUE")
创建矢量点(sunangles,'D:\\文献\\1-s2.0-S0034425717303991-mmc1\\point','I:\\HLSwGF6\\bx11_6\\S2B_MSIL1C_20201106T134209_N0500_R124_T22KCA_20230410T094912\\S2B_MSIL1C_20201106T134209_N0500_R124_T22KCA_20230410T094912.SAFE\\GRANULE\\L1C_T22KCA_A019169_20201106T134212\\IMG_DATA\\T22KCA_20201106T134209_B03.jp2')
#插值为30m的tif
arcpy.gp.Idw_sa("D:\\文献\\1-s2.0-S0034425717303991-mmc1\\point\\new_points.shp", "sunAzimuth", "D:\\文献\\1-s2.0-S0034425717303991-mmc1\\point\\output4.tif", 30, 2)

传感器的两个角度是从文献中下载的代码，贴出文章名可以去官网下载，亲测可以直接使用
Sentinel-2 MultiSpectral Instrument (MSI) data processing for aquatic science applications: Demonstrations and validations
Supplementary data to this article can be found online at http://dx.doi.org/10.1016/j.rse.2017.08.033.