Satellite Application Facility for Numerical Weather Prediction › Forums › RTTOV › RTTOV v13 › RTTOV v13 General Discussion › Calculating Satellite and Solar Angles in Python Wrapper
Tagged: RTTOV WRAPPER
- This topic has 3 replies, 2 voices, and was last updated 7 months, 4 weeks ago by James Hocking.
-
AuthorPosts
-
April 6, 2024 at 9:14 am #49417SWADHIN SATAPATHYParticipant
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?
April 6, 2024 at 9:15 am #49419James HockingKeymasterHi,
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,
JamesApril 8, 2024 at 3:52 pm #49431SWADHIN SATAPATHYParticipantHi 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.0if geo_sat_lat is None:
sat_lat = 0.0
else:
sat_lat = geo_sat_latif geo_sat_height is None:
sat_height = 35800.0
else:
sat_height = geo_sat_heightr_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_dlontan_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 + 180return sunzenangle, sunazangle
April 8, 2024 at 3:56 pm #49439James HockingKeymasterHi,
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.F90In 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 -
AuthorPosts
- You must be logged in to reply to this topic.