Calculate when a satellite horizon is overlapping target area on ground

Question

I’m trying to write some python code to work out when the Horizon of the International Space Station is overlapping with a target area on ground. In my case I’ve chosen China. I’ve got some code that works but I am not sure if I am mathematically approaching the problem correctly and could do with some advice.

The code I have uses shapely to define a rough polygon using polar coordinates for the national boundary of China. I then use ephem to compute the ground coordinates and elevation of the International Space Station. The elevation is in meters so I use that to calculate the distance to the horizon from an eyeball at that elevation (the radius of the horizon circle).

The idea next would be to create a circular polygon (.buffer(1)) with the International Space Station in the centre and simply see if this intersects with the polygon of China using shapely built in functions.

But it gets a bit tricky because the China polygon is in polar coordinates and the horizon distance is in meters. So, to account for the fact that lines of longitude are wide at the equator and narrow at the poles, I have to do a calculation to first work out the distance between lines of longitude at the current latitude and then use that to work out how many polar degrees the circle needs to be in x and y directions. I then make an ellipse of those dimensions and see if the intersection is there.

import time
import ephem
import math
from shapely.geometry import Polygon, Point
import shapely.affinity
import datetime

R = 6371 * 1000.0

dist_1_deg_long_equator = 111321
dist_1_deg_lat_equator = 111000

def get_dist_to_horizon(elevation_meters):
    eye = R + elevation_meters
    return R * math.acos(R / eye)

china = Polygon([
    (48.74, 87.17),
    (39.09, 74.05),
    (33.25, 79.13),
    (28.24, 86.44),
    (29.58, 95.74),
    (26.33, 98.82),
    (24.29, 97.89),
    (21.85, 100.82),
    (23.57, 105.32),
    (21.65, 108.14),
    (23.04, 116.23),
    (27.10, 120.31),
    (30.63, 122.08),
    (39.81, 124.09),
    (46.87, 133.90),
    (53.37, 121.73),
    (46.57, 119.57),
    (41.64, 105.26),
    (42.75, 96.28),
    (45.30, 90.76)])

name = "ISS (ZARYA)             "
line1 = "1 25544U 98067A   19094.21920345  .00002412  00000-0  46183-4 0  9994"
line2 = "2 25544  51.6444   8.9214 0002426 147.8175  11.8704 15.52483300163762"

iss = ephem.readtle(name, line1, line2)

while True:
    iss.compute()
    iss_horizon_radius = get_dist_to_horizon(iss.elevation)    

    dist_1_deg_long_current = math.cos(iss.sublat) * dist_1_deg_long_equator

    x_fact = iss_horizon_radius / dist_1_deg_long_current  # include more lines of longitude further from equator
    y_fact = iss_horizon_radius / dist_1_deg_lat_equator  # latitude lines are approx equal distant

    iss_horizon_circle = Point(math.degrees(iss.sublat), math.degrees(iss.sublong)).buffer(1)
    iss_horizon_ellipse = shapely.affinity.scale(iss_horizon_circle, x_fact, y_fact, origin='center')

    if iss_horizon_ellipse.intersects(china):
        print("ISS horizon is over China! %s" % datetime.datetime.utcnow())

    time.sleep(30)

Is there a better way to do this? For instance convert everything into Cartesian and then just have the horizon as a perfect circle? Any advice would be appreciated! Thanks in advance