With reference to #2244. Specifically addressing @mhvk's point 1 in #2244 (comment). It would be very useful if we could add the above frames to the coordinates package. I don't think this should be too hard - I've included a first stab below.

The transformation framework allows you to calculate the location of the observatory in GCRS coordinates (see below). @eteq's transformation agrees with the code I posted in #2244 (comment) to within 70m, so that's a tiny part of #4268 ticked off.

from astropy.coordinates import sites
import astropy.coordinates as coord
from astropy.time import Time
from astropy import units
site_list = sites.get_builtin_sites()
WHT = site_list['lapalma']
now  = Time.now()
itrs = coord.ITRS(WHT.itrs.cartesian,obstime=now)
gcrs = itrs.transform_to(coord.GCRS(obstime=now))

The position and velocity of the heliocenter and barycenter can be found using erfa.epv00, which is is consistent with the JPL ephemeris within 11.2 km = 4micro-sec and 5 mm/s, so the addition of a Heliocentric and Barycentric frame to the coordinates package doesn't look hard, e.g

def gcrs_to_bcrs(gcrs_loc):
    tdb = gcrs_loc.obstime.tdb
    #returns BCRS position and velocity of Earth in AU and AU/d
    h_pv,b_pv = erfa.epv00(tdb.jd1,tdb.jd2)    
    # assume scalar obstime! and unpack position array
    b_pos = b_pv[0,:] * units.AU
    new_rep = CartesianRepresentation(gcrs.cartesian.xyz + b_pos)
    return BCRS(new_rep,obstime=gcrs_loc.obstime)

One snag with the example above is that it assumes the obstime of the frame is a scalar.