EDIT: solved, thanks to Eric Jensen's suggestion (in a comment on his answer) that I include the correction to my right ascension value as a function of my declination. All other comments were helpful but included information I'd already referenced, it was the the correction to RA that I was missing.
As the title states, I'm trying to use the WCS data added to a FITS header by a plate solve to calculate RA and Dec at any given pixel for an image in decimal form. I am aware that programs like AstroImageJ can display this sort of info, but only in HHMSS/DDMMSS format and only by using the cursor. I'd like to automate this with my own code.
Problem summary:
I have a number of plate solved images of star fields with satellite streaks, and I'd like to automate the extraction of satellite RA/Dec from the images based on their positions within the FITS pixel space with a function that allows me to input FITS header data and pixel location to generate RA/Dec. I know that this is a solved problem because AstroImageJ can give me truth data from a plate solved image when I mouse over certain pixels, but I have so many satellite streaks to measure that using the cursor for RA/Dec readout for each point along them would be prohibitively time consuming.
What I've tried so far:
I have tried to apply the transformations from the FITS header data based on this Caltech guide, but the algorithm here does not match my data. In particular, the RA/Dec I calculate for a pixel's position using this algorithm is offset by up to an arcminute (in a ~1* FOV frame) from the RA/Dec I get using my cursor in AIJ to find the RA/Dec of the same spot, and I haven't managed to figure out why these transformations aren't working.
I know that these transformations assume a linear model and don't directly include distortion, so I included a correction for the SIP distortion but found that the distortions were sub-pixel in my images and couldn't account for the offset.
I then found this StackExchange post of a user with a similar problem and attempted to apply a correction for the gnomonic projection based on this algorithm, and while that seemed to work, I'm still off by a few arcminutes.
So, after using the FITS algorithm, SIP distortion correction, and correcting for the gnomonic projection, I've accounted for everything I can think of but I still can't match AstroImageJ's numbers for a plate solved image. Can anyone help me find the missing piece here?
I would also be open to suggestions on alternate methods to get satellite streak RA/Dec from plate solved images, in case there's a simpler approach I'm missing.