太阳天顶角SZA的计算

#encoding=utf-8
import xlutils
import math
import datetime
import pandas as pd

hourlist = ['00','01','02','03','04','05','06','07','08','09','10','11','12','13','14','15','16','17','18','19','20','21','22','23']
begin = datetime.date(2019, 8, 6)
end = datetime.date(2022, 12, 31)
for i in range((end - begin).days + 1):
    date = begin + datetime.timedelta(days=i)
    year = str(date).split('-')[0]
    month = str(date).split('-')[1]
    day = str(date).split('-')[2]
    for hour in hourlist:
            filename = r"C:\Users\point.csv"
            data = pd.read_csv(filename, header=0)
            df = pd.DataFrame(data)
            nrows = df.shape[0]  # 行数
            ncols = df.shape[1]  # 列数
            df.insert(loc=5, column='sza',value=-1)

            for row in range(nrows):
                station = df.iloc[row]['pointid']
                lon = df.iloc[row]['lon']
                lat = df.iloc[row]['lat']
                def tt(currentLon):
                    shangValue = int(currentLon / 15)
                    yushuValue = math.fabs(currentLon % 15)
                    if (yushuValue <= 7.5):
                        timeZone = shangValue
                    else:
                        if (currentLon > 0):
                            timeZone = shangValue + 1
                        else:
                            timeZone = shangValue - 1
                    return timeZone
                TimeZone = tt(lon)
                # 年积日DOY的计算。用到儒略日 Julian day(由通用时转换到儒略日)
                JD0 = int(365.25 * (float(year) - 1)) + int(30.6001 * (1 + 13)) + 1 + float(hour) / 24 + 1720981.5
                if float(month) <= 2:
                    JD2 = int(365.25 * (float(year) - 1)) + int(30.6001 * (float(month) + 13)) + float(day) + float(
                        hour) / 24 + 1720981.5
                else:
                    JD2 = int(365.25 * float(year)) + int(30.6001 * (float(month) + 1)) + float(day) + float(
                        hour) / 24 + 1720981.5
                DOY = JD2 - JD0 + 1

                # N0 sitar=θ
                # 太阳赤纬角ED1转为黄赤交角ED
                N0 = 79.6764 + 0.2422 * (float(year) - 1985) - int((float(year) - 1985) / 4.0)
                sitar = 2 * math.pi * (DOY - N0) / 365.2422
                ED1 = 0.3723 + 23.2567 * math.sin(sitar) + 0.1149 * math.sin(2 * sitar) - 0.1712 * math.sin(
                    3 * sitar) - 0.758 * math.cos(sitar) + 0.3656 * math.cos(2 * sitar) + 0.0201 * math.cos(
                    3 * sitar)
                ED = ED1 * math.pi / 180  # ED本身有符号

                # # (经度纠正）小时尺度，不需要
                # if lon >= 0:
                #     if TimeZone == -13:
                #         dLon = lon - (math.floor((lon*10-75)/150)+1)*15.0
                #     else:
                #         dLon = lon - TimeZone*15.0 #地球上某一点与其所在时区中心的经度差
                # else:
                #     if TimeZone == -13:
                #         dLon = (math.floor((lon*10-75)/150)+1)*15.0- lon
                #     else:
                #         dLon = TimeZone*15.0- lon

                Eq = 0.0028 - 1.9857 * math.sin(sitar) + 9.9059 * math.sin(2 * sitar) - 7.0924 * math.cos(
                    sitar) - 0.6882 * math.cos(2 * sitar)  # 时差
                # gtdt1 = hour[m] + min[m]/60.0 + sec[m]/3600.0 + dLon/15 #地方时
                gtdt1 = float(hour) + TimeZone  # 地方时
                gtdt = gtdt1 + Eq / 60.0  # 真太阳时
                dTimeAngle1 = 15.0 * (gtdt - 12)  # 太阳时角
                dTimeAngle = dTimeAngle1 * math.pi / 180  # 太阳时角
                latitudeArc = lat * math.pi / 180  # 纬度 -90~90
                # 天顶角计算公式
                CosZenithAngle = (
                        math.sin(ED) * math.sin(latitudeArc) + math.cos(ED) * math.cos(latitudeArc) * math.cos(
                    dTimeAngle))
                ZenithAngleARC = math.acos(CosZenithAngle)
                ZenithAngle = ZenithAngleARC * 180 / math.pi

                print('id：' + str(station) + ' 太阳天顶角(deg)：%f' % (ZenithAngle))

                df.loc[row, 'sza'] = ZenithAngle
            new_name ="E:\\test\\"+year+month+day+hour+'.csv'
            df.to_csv(new_name, index=False)
            print(new_name)