python代码计算太阳高度角和太阳方位角

simone's record

已于 2024-04-12 14:15:16 修改

阅读量1.5k

点赞数 18

文章标签： python 学习方法

于 2024-04-12 12:30:42 首次发布

本文链接：https://blog.csdn.net/mumu__fairy/article/details/137675319

版权

本文详细介绍了如何使用Python的ephem和pysolar库以及手动计算方法来求解给定经纬度在特定日期和时间的太阳高度角和方位角。作者通过比较不同方法的结果，强调了手算方法5的准确性，并提供了数学公式和示例代码。

摘要由CSDN通过智能技术生成

计算(109.2, 35.1) 2021/12/22 09:00:00太阳高度角和太阳方位角

本文使用了5种方式对其进行求解，最终应用了方法5的运算结果，读者可直接跳至方法5查看结果和python代码

计算(109.2, 35.1) 2021/12/22 09:00:00太阳高度角和太阳方位角

方法1：使用python库ephem

结果：Sun altitude: 13:33:46.2, Sun azimuth: 267:57:01.4

import ephem  # 获取特定地点和日期时间的太阳高度角和方位角

def calculate_sun_position(latitude, longitude, date):
    obs = ephem.Observer()
    obs.lat, obs.lon = latitude, longitude
    obs.date = ephem.Date(date)
    sun = ephem.Sun()
    sun.compute(obs)
    return sun.alt,sun.az  # 返回太阳高度角和方位角

latitude = 35.1
longitude = 109.2
date = '2021/12/22 9:00:00'  # 某一天的正午时间
sun_altitude, sun_azimuth = calculate_sun_position(latitude, longitude, date)
print(f"Sun altitude: {sun_altitude}, Sun azimuth: {sun_azimuth}")

方法2：使用python库pysolar

结果：Sun altitude: 10.52866037565445, Sun azimuth: 128.60185059293036

import datetime
import pysolar.solar
from pytz import timezone
from timezonefinder import TimezoneFinder

# 定义观测点经纬度
latitude = 35.1  # 纬度#38.053839
longitude = 109.2  # 经度#106.697611
elevation = 1000 # 高度

# 创建一个TimezoneFinder实例
tf = TimezoneFinder(in_memory=True)
# 使用TimezoneFinder查找对应的时区名称
timezone_name = tf.timezone_at(lat=latitude, lng=longitude)
# 将时区名称转换为pytz对象
timezone_obj = timezone(timezone_name)
print(f"经纬度({latitude}, {longitude})所处的时区是：{timezone_name}")
# 获取当前日期和时间
date_time = datetime.datetime(2021,12,22,9,0,0,0,timezone_obj)

# 计算当前时刻的太阳高度角（单位：度）
sun_altitude = pysolar.solar.get_altitude(latitude, longitude, date_time, elevation)
sun_azimuth = pysolar.solar.get_azimuth(latitude, longitude, date_time, elevation)
print(f"Sun altitude: {sun_altitude}, Sun azimuth: {sun_azimuth}")

方法3：使用网页在线计算器1

结果：Sun altitude: 20.15858896258514, Sun azimuth: 140.893388236611

方法4：使用网页在线计算器2

结果：Sun altitude: 11.39, Sun azimuth: 50.40

方法5：根据公式手算

结果：Sun altitude: 11.13924683296049, Sun azimuth: 50.66025597336174

根据网上找到的例子，使用太阳等高线原理计算(30N, 120E)，夏至日10:00的太阳高度角是62.5°。使用解天文三角形的公式代入进行验证，太阳高度角≈62.5°，该法亦接近方法4所求结果，暂定此法为准确结果，转为python代码如下：

# -*- coding:utf-8 -*-
import math
from datetime import datetime

ASIN = math.asin
SIN = math.sin
COS = math.cos
RADIANS = math.radians
DEGREES = math.degrees

#时差: 
def TimeDifference(lon):
    t = (lon -120) / 15
    return t

#真太阳时(T：北京时间)
def TrueSolarTime(T, t):
    ST = T + t
    return ST

#时角
def HourAngle(ST):
    h = 15 * (ST -12)
    return h

#太阳赤纬(N为积日，一年中的第几天)
def SolarDeclination(N):
    DEC = -23.443 * COS(RADIANS(round(360/365, 8) * (N + 10))) #三角函数需要用弧度制
    return DEC
    
#积日
def DayOfYear(time):
    day_of_year = time.timetuple().tm_yday  # 获取当年的第几天
    return day_of_year

#计算时角
def CalculateHourAngle(lon, T):
    t = TimeDifference(lon)
    ST = TrueSolarTime(T, t)
    H = HourAngle(ST)
    return H
    
#计算太阳赤纬
def CalculateSolarDeclination(time):
    N = DayOfYear(time)
    DEC = SolarDeclination(N)
    return DEC
    
def CalculateSunPosition(latitude, longitude, time):
    T = time.hour + round(1/60, 8) * time.minute + round(1/3600, 8) * time.second
    print("Beijing Time: ", T)
    H = CalculateHourAngle(longitude, T)
    print("Hour Angle: ", H)
    DEC = CalculateSolarDeclination(time)
    print("DEC: ", DEC)
    sun_altitude = DEGREES(ASIN(SIN(RADIANS(latitude))*SIN(RADIANS(DEC))+COS(RADIANS(DEC))*COS(RADIANS(latitude))*COS(RADIANS(H))))
    sun_azimuth = DEGREES(ASIN(round((-COS(RADIANS(DEC))*SIN(RADIANS(H)) / COS(RADIANS(sun_altitude))), 8)))
    print("altitude_angle_deg:", sun_altitude)
    print("azimuth_angle_deg:", sun_azimuth)
    return sun_altitude, sun_azimuth


time_str = '2021-12-22 9:00:00'
time = datetime.strptime(time_str, '%Y-%m-%d %H:%M:%S')
sun_altitude, sun_azimuth = CalculateSunPosition(35.1, 109.2, time)

运算结果：