I'm happy to see OP has used Brandon Rhodes' sgp4 python library to obtain periapsis and apoapsis from a TLE. I've always used the skyfield python library for this which is from the same developer and builds on the sgp4 library but has way more features.
From https://celestrak.org/NORAD/elements/table.php?tleFile=stations&title=Space%20Stations&orbits=1&pointsPerRev=90&frame=1 as of right now:
Int. Des. Cat No. name T (min) i (deg) Apo x Peri (km) ecc TLE age (day)
1998-067A 25544 ISS (ZARYA) 92.96 51.64 423 x 416 0.0004540 0.33
which cites the most recent TLE
ISS (ZARYA)
1 25544U 98067A 21356.70730882 .00006423 00000-0 12443-3 0 9994
2 25544 51.6431 130.5342 0004540 343.5826 107.2903 15.49048054317816
What I'll do is propagate the satellite for +/- 2 days around it's epoch (time of best accuracy) and ask for the apo and peri from sgp4 and also plot the 'altitude' and the 'height' above Earth's surface at 30 second intervals to sample at a fine enough step size so that we can just use a simple numerical extrema detector to get pretty close to the perigees and apogees.
What's interesting is that these do not at all coincide with the maxima and minima of elevation above the Earth's surface! Since the orbit is very very close to circular, the Earth's equatorial bulge tends to force a second minimum altitude near the ascending and descending nodes, separate from the perigee.
What is even more interesting is that I can not reproduce Celestrak's 425 x 416!
For orbital altitudes that go into things like 423 x 416 km, we take the distances from Earth's geocenter and subtract 6378.137 kilometers; Earth's geoid's equatorial radius.
This is simply a convention.
We can also see the eccentricity is slowly increasing over a four day period centered on the TLE's epoch.

import numpy as np
import matplotlib.pyplot as plt
from skyfield.api import load, EarthSatellite, Loader, wgs84
import datetime
loaddata = Loader('~/Documents/fishing/SkyData') # avoids multiple copies of large files
ts = loaddata.timescale() # include builtin=True if you want to use older files (you may miss some leap-seconds)
eph = loaddata('de421.bsp') # a small one, fine for this
satname, L1, L2 = """ISS (ZARYA)
1 25544U 98067A 21356.70730882 .00006423 00000-0 12443-3 0 9994
2 25544 51.6431 130.5342 0004540 343.5826 107.2903 15.49048054317816""".splitlines()
sat = EarthSatellite(L1, L2, satname)
hw_days, dt_sec = 2.0, 10.
daysec = 24 * 3600 # seconds
earth_eq_radius = 6378.137 # km
days = np.arange(-hw_days * daysec, hw_days * daysec + 1, dt_sec) / daysec
times = ts.tt_jd(sat.epoch.tt + days)
g = sat.at(times) # geocentric position object
p = g.position.km # geocentric positions (km)
d = g.distance().km
h = d - earth_eq_radius
el = wgs84.subpoint(g).elevation.km # height above the wgs84 geoid
def get_maxes_and_mins(a):
"""simple imperfect way to find points above or below both of those on either side
and misses those where two are exactly equal as well as possible endpoints"""
maxes = np.where((a[1:-1] > a[2:]) * (a[1:-1] > a[:-2]))[0]
mins = np.where((a[1:-1] < a[2:]) * (a[1:-1] < a[:-2]))[0]
return maxes, mins
h_maxes, h_mins = get_maxes_and_mins(h)
el_maxes, el_mins = get_maxes_and_mins(el)
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=[12, 6])
ax1.plot(days, h, label='height')
ax1.plot(days, el, label='elevation')
ax1.plot(days[h_maxes], h[h_maxes], '.k')
ax1.plot(days[h_mins], h[h_mins], '.k')
ax1.plot(days[el_maxes], el[el_maxes], '.k')
ax1.plot(days[el_mins], el[el_mins], '.k')
ax1.set_xlim(-0.5, 0.5)
ax1.set_ylim(410, 440)
ax1.legend()
ax2.plot(days[h_maxes], h[h_maxes], '-')
ax2.plot(days[h_mins], h[h_mins], '-')
ax2.plot(days[el_maxes], el[el_maxes], '-')
ax2.plot(days[el_mins], el[el_mins], '-')
ax2.set_xlim(-hw_days, hw_days)
ax2.set_ylim(412, 437)
epoch_string = datetime.datetime(*sat.epoch.utc[:5]).isoformat().replace('T', ' ') + 'UTC'
ax2.set_xlabel('days relative to TLE epoch')
plt.suptitle('ISS TLE epoch: ' + epoch_string)
plt.show()