Compute off-nadir angle of vector from satellite to ground

I'd like to write Python code to specify the orbit of a satellite with Keplerian elements, specify a point on the Earth with latitude, longitude, and altitude, specify a time, and compute the angle between two vectors:
- the vector from the satellite to the specified point on the Earth
- the vector from the satellite to the center of the Earth.

I know I can use poliastro to define the orbit and propagate it to the specified time. The hard part is representing the satellite and the Earth point in the same coordinate system.

poliastro currently doesn't specify a coordinate system. Someone in their chat room told me that Earth orbits are in GCRS. astropy can convert GCRS to ITRS, which is an Earth-centered Earth-fixed frame:

import math
import numpy as np
from astropy import units as u
from astropy.time import Time
from poliastro.bodies import Earth
from poliastro.twobody import Orbit
from astropy.coordinates import SkyCoord

def lla2ecef(lat, lon, alt):
    """ Convert lat/lon/alt to cartesian position in ECEF frame.

    Origin is center of Earth. +x axis goes through lat/lon (0, 0).
    +y axis goes through lat/lon (0, 90). +z axis goes through North Pole.

    lat: number, geodetic latitude in degrees
    lon: number, longitude in degrees
    alt: number, altitude above WGS84 ellipsoid, in km

    Returns: tuple (x, y, z) coordinates, in km.

    Source: "Department of Defense World Geodetic System 1984"
    Page 4-4
    National Imagery and Mapping Agency
    Last updated June, 2004
    NIMA TR8350.2
    """

    lon = lon * math.pi/180.0  # Convert to radians
    lat = lat * math.pi/180.0  # Convert to radians

    # WGS84 ellipsoid constants:
    a = 6378.137 #equatorial radius, in km
    e = 8.1819190842622e-2

    # intermediate calculation: prime vertical radius of curvature
    N =  a/math.sqrt(1 - e**2*math.sin(lat)**2)

    #results
    x = (N + alt)*math.cos(lat)*math.cos(lon)
    y = (N + alt)*math.cos(lat)*math.sin(lon)
    z = ((1 - e**2)*N + alt)*math.sin(lat)

    return (x, y, z)

epoch = Time(2018, format='decimalyear', scale='tai') #01-01-2018 00:00:00, in TAI
propagation_time = 9000 #seconds
semi_major_axis = 10000 #km
eccentricity = 0.1
inclination = 50 #deg
raan = 70 #deg
arg_perigee = 60 #deg
true_anomaly = 80 #deg
orbit = Orbit.from_classical(
    Earth, 
    semi_major_axis*u.km, 
    eccentricity*u.one,  
    inclination*u.deg, raan*u.deg, 
    arg_perigee*u.deg, 
    true_anomaly*u.deg, 
    epoch)
propagated_orbit = orbit.propagate(propagation_time*u.s)
pos_gcrs = propagated_orbit.state.r
sky_gcrs = SkyCoord(
    representation_type='cartesian', 
    x=pos_gcrs[0], y=pos_gcrs[1], z=pos_gcrs[2],
    frame='gcrs',
    obstime=(epoch + propagation_time*u.s))
pos_ecef = sky_gcrs.transform_to('itrs')
pos_s = np.array((pos_ecef.x.to(u.km).value, 
                  pos_ecef.y.to(u.km).value, 
                  pos_ecef.z.to(u.km).value))

lat = 40 #deg
lon = 50 #deg
alt = 0.06 #km
pos_t = np.array(lla2ecef(lat, lon, alt))

#Compute angle at satellite between target and center of earth
v1 = pos_t - pos_s
v2 = -pos_s
angle = math.acos(np.dot(v1, v2)/(np.linalg.norm(v1)*np.linalg.norm(v2)))

#convert to degrees
angle = angle*180.0/math.pi