Calculating Satellite and Solar Angles in Python Wrapper

[SWADHIN SATAPATHY]

Hi All,

I recently came across the updated version of RTTOV, which can calculate the satellite and solar angles given the latitude, longitude, and datetimes (including some other optional parameters). I use the Python wrapper PyRTTOV to generate the profiles and run RTTOV as well. Can anyone help me with how to set the data as required by by those specific subroutines and how to call it using the wrapper?

[James Hocking]

Hi,

I’m afraid this isn’t possible as we have not implemented interfaces for these angle calculation routines in the wrapper.

I will look at implementing this in a future release.

Best wishes,
James

[SWADHIN SATAPATHY]

Hi James,

Thanks for the response. For now I have created a couple of functions which basically are supposed to do the same things as the subroutines. It’ll be great if you can look into them and let me know if there’s any issue.

def calculate_geo_sat_angles(lat, lon, geo_sat_lon, geo_sat_lat=None, geo_sat_height=None):
earth_radius = 6371.0 # km
deg2rad = np.pi / 180.0
rad2deg = 180.0 / np.pi
zenmaxv9 = 86.0
zenmax = 88.0

if geo_sat_lat is None:
sat_lat = 0.0
else:
sat_lat = geo_sat_lat

if geo_sat_height is None:
sat_height = 35800.0
else:
sat_height = geo_sat_height

r_ratio = earth_radius / (sat_height + earth_radius)

# max_sat_zen = max_sat_zen * (1.0 – 1e-5)

cos_dlon = cos(radians(lon – geo_sat_lon))

cos_alpha = sin(radians(lat)) * sin(radians(sat_lat)) + \
cos(radians(lat)) * cos(radians(sat_lat)) * cos_dlon

tan_satzen_denom = cos_alpha – r_ratio

if tan_satzen_denom > 0:
sin_alpha = np.sqrt(1.0 – cos_alpha * cos_alpha)
tan_satzenithang = sin_alpha / tan_satzen_denom
sat_zen = degrees(atan(tan_satzenithang))
else:
sat_zen = 90.0

# if sat_zen > max_sat_zen:
# raise ValueError(“Satellite zenith angle exceeds RTTOV maximum”)

sin_lat = sin(radians(lat))
sin_dlon = sin(radians(lon – geo_sat_lon))

if sin_dlon == 0.0 and -sin_lat * cos_dlon == 0.0:
sat_azim = 360.0 # nadir, azimuth is undefined
else:
sat_azim = degrees(atan2(sin_dlon, -sin_lat * cos_dlon))

sat_azim = 360.0 – sat_azim if sat_azim < 0 else sat_azim

return sat_zen, sat_azim

def equation_of_time(date):
# Calculate Julian day
year, month, day = date
if month < 3:
year -= 1
month += 12
A = year // 100
B = 2 – A + (A // 4)
julian_day = math.floor(365.25 * (year + 4716)) + math.floor(30.6001 * (month + 1)) + day + B – 1524.5

# Calculate mean solar anomaly
mean_solar_anomaly = (0.9856 * (julian_day – 2451545)) % 360

# Calculate equation of center
equation_of_center = (1.9148 * math.sin(math.radians(mean_solar_anomaly))) + \
(0.02 * math.sin(math.radians(2 * mean_solar_anomaly))) + \
(0.0003 * math.sin(math.radians(3 * mean_solar_anomaly)))

# Calculate ecliptic longitude
ecliptic_longitude = (mean_solar_anomaly + equation_of_center + 180 + 102.9372) % 360

# Calculate equation of time
obliquity_of_ecliptic = 23.4397
y = math.tan(math.radians(obliquity_of_ecliptic / 2)) ** 2
eccentricity_earth_orbit = 0.0167
equation_of_time = 4 * (math.degrees(y * math.sin(2 * math.radians(ecliptic_longitude))) –
(eccentricity_earth_orbit * math.sin(math.radians(mean_solar_anomaly))) +
(0.020 * math.sin(math.radians(2 * mean_solar_anomaly))))

# Calculate sun declination
sin_declination = math.sin(math.radians(math.asin(math.sin(math.radians(obliquity_of_ecliptic)) *
math.sin(math.radians(ecliptic_longitude)))))

return equation_of_time, sin_declination

def calculate_solar_angles(latitude, longitude, year, month, day, hour, minute, second):
# Calculate equation of time and sun declination
eqoftime, sundeclination = equation_of_time((year, month, day))

# Ensure valid latitude and longitude
if latitude < -90 or latitude > 90:
raise ValueError(“Latitude must be in the range [-90, 90] degrees”)
longitude = longitude % 360.0

# Calculate true solar time and solar hour angle
timeoffset = eqoftime + (4 * longitude)
truesolartime = (hour * 60) + minute + (second / 60) + timeoffset
hourangle = (truesolartime / 4 – 180) * math.pi / 180

# Calculate trigonometric values
sin_lat = math.sin(latitude * math.pi / 180)
cos_lat = math.cos(latitude * math.pi / 180)
sin_dec = sundeclination
cos_dec = math.cos(math.asin(sundeclination))
sin_hour = math.sin(hourangle)
cos_hour = math.cos(hourangle)

# Calculate solar zenith angle
cos_theta = (sin_lat * sin_dec) + (cos_lat * cos_dec * cos_hour)
if abs(cos_theta) > 1:
cos_theta = math.copysign(1, cos_theta)
sunzenangle = math.acos(cos_theta) * 180 / math.pi

# Calculate solar azimuth angle
tan_azim_num = cos_dec * sin_hour
tan_azim_den = (sin_lat * cos_dec * cos_hour) – (cos_lat * sin_dec)
sunazangle = math.atan2(tan_azim_num, tan_azim_den) * 180 / math.pi + 180

return sunzenangle, sunazangle

[James Hocking]

Hi,

If you have ported the calculations correctly then they should be fine.

If you want to validate your code you could look at the unit tests for these two routines:
src/test/rttov_test_calc_solar_angles.F90
src/test/rttov_test_calc_geo_sat_angles.F90

In both cases you can see the test inputs and reference output values in the DATA statements in these source files. Your routines should replicate the outputs for the given inputs.

Best wishes,
James

[Author]

Posts