cut from https://forum.cosmoquest.org/showthread.php?137475-How-to-cut-out-a-region-from-an-SDSS-FITS-image-file

The coordinate information in FITS files is encoded in sets of keywords describing the world coordinate system (WCS), which may take se several forms. For direct images, these typically specify a reference point in RA and dec, the pixel location of this reference point, the pixel scale (which may in general differ for each axis), and the orientation of the pixel grid on the sky. In SDSS files, they look like this:
CTYPE1  = 'RA---TAN'
CTYPE2  = 'DEC--TAN'
CUNIT1  = 'deg     '
CUNIT2  = 'deg     '
CRPIX1  = 1.02450000000000E+03 / Column Pixel Coordinate of Ref. Pixel
CRPIX2  = 7.44500000000000E+02 / Row Pixel Coordinate of Ref. Pixel
CRVAL1  = 2.47486643250000E+02 / RA at Reference Pixel
CRVAL2  = 2.43901742700000E+01 / DEC at Reference Pixel
CD1_1   = 6.89565533678154E-05 / RA  degrees per column pixel
CD1_2   = 8.57010978890200E-05 / RA  degrees per row pixel
CD2_1   = 8.56370751953130E-05 / DEC degrees per column pixel
CD2_2   = -6.9012795698925E-05 / DEC degrees per row pixel

As is common, the coordinates assume the so-called tangent-plane projection from the celestial sphere to the flat detector (which is fine for small fields of view unless the optics have very strong radial distortion). CDn_m specifies the change in celestial coordinate n due to a pixel change in image coordinate m; CD1_2 and CD2_1 are 0 when the pixel grid is aligned with the coordinate axes.

The traditional aproach in astrometry is to define standard coordinates xi, eta:
Xi = cd1_1*(x-crpix1) + cd1_2*(y-crpix2)
Eta = cd2_1*(x-crpix1) + cd2_2*(y-crpix2)

and if there is no further distortion (or it has been corrected by resampling the image) transform those coordinates into angular ones:
cot δ sin (RA - CRVAL1) = (ξ) / (sin CRVAL2 + η cos CRVAL2) and cot δ cos (RA - CRVAL1) = (cot CRVAL2 - η sin CRVAL2) / (sin CRVAL2 + η cos CRVAL2)
(I see the Greek letters in those; if they're lost, these are in the "Narrow-field astrometry" section here). I think I got the translation from RA, dec symbols to CRVALx properly... For very small coordinate differences, one can often simply pretend everything is a linear transformation from a notional rectangular RA/dec grid to pixel space.

The defining document for representing celestial coordinates in FITS is by Greisen and Calabretta; most of it deals with assorted all-sky projections rather than narrow-field dorect imaging, where the transformation is rather simpler.

ETA: You may have seen that the various SDSS filter images on a given field have registration differences of several pixels, which is why this exercise needs doing to get matched subimages. For such small offsets, you can use the simplest Cartesian approximation to much better than a pixel accuracy. The important information is, in this case with the same pixel scale and orientation for each image, carried in the CRVAL and CRPIX sets of keywords.