I'm no SPICE expert but here are some potential solutions (unless of course you want to try Skyfield's Almanac methods!)
Firsrt possibility, but all these may not work
this answer links to this comment links to SPICE Tutorials, Updated December 11, 2019 links to the 69 slide presentation SPICE Geometry Finder (GF) Subsystem; Searching for times when specified
geometric conditions occur dated January 2020.
- Much SPICE software computes a geometry parameter at a given time, t, i.e. x = f(t). Example: on 2011 MAR 30 14:57:08, what is the spacecraft’s altitude
above Mars?
- The Geometry Finder subsystem does the inverse: it finds times when specified geometric events occur
GF provides two primary types of event-finding APIs
- Boolean: a geometric condition (an event) is true or false
– Numeric: a geometric quantity has a given value, is within a given range or has achieved a local or global maximum or minimum
In your case you would specify a position on a body, then look for times when elevation was 0 or 0+/- small.
Another possibility, this seems to address exactly what you are looking for
https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/spkcpo_c.html
Examples
The numerical results shown for these examples may differ across
platforms. The results depend on the SPICE kernels used as
input, the compiler and supporting libraries, and the machine
specific arithmetic implementation.
Compute apparent solar azimuth and elevation as seen from a
specified surface point on the earth.
Task Description
In this example we'll use the location of the DSN station
DSS-14 as our surface point.
We'll perform the solar azimuth and elevation computation two
ways:
Using a station frame kernel to provide the
specification of a topocentric reference frame
centered at DSS-14.
Computing inline the transformation from the earth-fixed,
earth-centered frame ITRF93 to a topocentric frame
centered at DSS-14.
Note that results of the two computations will differ
slightly. There are three sources of the differences:
The station position is time-dependent due to tectonic
plate motion, and epochs of the station positions used
to specify the axes of the topocentric frame are
different in the two cases. This gives rise to different
orientations of the frame's axes relative to the frame
ITRF93.
The two computations use different earth radii; this
results in computation of different geodetic latitudes
of the station. This difference also affects the
topocentric frame orientation relative to ITRF93.
The station movement between ET and the epoch at which
the DSS-14_TOPO frame is specified contributes a very
small offset---on the order of 10 cm---to the station-sun
position vector, expressed in the ITRF93 frame.
and there's a lot more there.
It seems you can do the same thing in Python using SpicyPy
Find the azimuth and elevation of the apparent position of the Moon as seen from the DSN station DSS-13 by the following steps:
- Find the apparent position vector of the Moon relative to the DSN station DSS-13 in the topocentric reference frame DSS-13_TOPO at epoch ET. Use light time and stellar aberration corrections.
For this step, you'll need to have loaded a station SPK file providing geocentric station position vectors, as well as a frame kernel specifying topocentric reference frames centered at the respective DSN stations. (Other kernels will be needed as well; you must choose these.)
Convert the position vector to latitudinal coordinates. Use the routine spiceypy.reclat for this computation.
Compute the Moon's azimuth and elevation as follows: azimuth is the negative of topocentric longitude and lies within the range 0-360 degrees; elevation is equal to the topocentric latitude. Display the results in degrees.