How to cut out a region from an SDSS FITS image file?

[Jean Tate]

FITS files of Fields contain the science data in SDSS images, one per band (or filter).

Using PyFITS, you can select the image data for any region within a 1489 x 2048 pixel FITS image. Five times (one for each band/filter).

But how to select the (five) regions so they are each centered on the same (RA, Dec) point (within half a pixel)? ???

I'm guessing you can work out the position, in pixel coordinates, of a particular (RA, Dec) point using the Header fields "RA" ("1st row - Right ascension of telescope boresigh"), "Dec" ("1st row - Declination of telescope boresight (d"), "SPA" ("1st row - Camera col position angle wrt north ("), etc ... but how? ??? ???

I just posted this in the Galaxy Zoo forum. I'll be interested to see how the responses differ ...

[ngc3314]

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.

[StupendousMan]

It is in many cases easier to do this the hard way. Extract five individual FITS images. Examine them by eye, and pick out one reasonably bright object which appears in all five images. Measure the (x, y) coordinates of this object in each image. Then, if you're lucky -- and for the SDSS, you probably WILL be lucky -- the scale and rotation of all five images will be almost exactly the same. In order to register the images, all you need to do is to shift each image a little bit in (x, y) so that the chosen object has the same position.

Do so.

In general, of course, you would need to use the full-blown WCS information, as ngc3314 has described. In the case of the SDSS, this shouldn't be necessary.

[Jean Tate]

A very big THANK YOU to you both! :-*

ngc3314, I particularly like the breadth of your post; in terms of understanding it was just what I was looking for (I like to understand how the motor works, not just know how to turn the key to start it).

StupendousMan, that is one cool approach! 8) I just love these kinds of "if you don't have a ruler handy, here's how you can measure distance" sorts of things!

You guys make it a great pleasure to learn how to *DO* astronomy, with the fantastic data which is SDSS (and 2MASS and HST and ...), without the benefit of being a full-time student in a formal university program.