CODEWITHSUNDEEP
pythonscipyfitsastropy
Wolfgang
Wolfgang

Reputation: 349

Splicing image array (FITS file) using coordinates from header

I am trying to splice a fits array based on the latitudes provided from the Header. However, I cannot seem to do so with my knowledge of Python and the documentation of astropy. The code I have is something like this:

from astropy.io import fits
import numpy as np

Wise1 = fits.open('Image1.fits')
im1 = Wise1[0].data

im1 = np.where(im1 > *latitude1, 0, im1)

newhdu = fits.PrimaryHDU(im1)
newhdulist = fits.HDUList([newhdu])
newhdulist.writeto('1b1_Bg_Removed_2.fits')

Here latitude1 would be a value in degrees, recognized after being called from the header. So there are two things I need to accomplish:

  1. How to call the header to recognize Galactic Latitudes?
  2. Splice the array in such a way that it only contains values for the range of latitudes, with everything else being 0.

Upvotes: 0

Views: 1973

Answers (1)

keflavich
keflavich

Reputation: 19205

I think by "splice" you mean "cut out" or "crop", based on the example you've shown.

astropy.nddata has a routine for world-coordinate-system-based (i.e., lat/lon or ra/dec) cutouts

However, in the simple case you're dealing with, you just need the coordinates of each pixel. Do this by making a WCS:

from astropy import wcs
w = wcs.WCS(Wise1[0].header)
xx,yy = np.indices(im.shape)
lon,lat = w.wcs_pix2world(xx,yy,0)

newim = im[lat > my_lowest_latitude]

But if you want to preserve the header information, you're much better off using the cutout tool, since you then do not have to manually manage this.

from astropy.nddata import Cutout2D
from astropy import coordinates
from astropy import units as u

# example coordinate - you'll have to figure one out that's in your map
center = coordinates.SkyCoord(mylon*u.deg, mylat*u.deg, frame='fk5')

# then make an array cutout
co = nddata.Cutout2D(im, center, size=[0.1,0.2]*u.arcmin, wcs=w)

# create a new FITS HDU
hdu = fits.PrimaryHDU(data=co.data, header=co.wcs.to_header())

# write to disk
hdu.writeto('cropped_file.fits')

An example use case is in the astropy documentation.

Upvotes: 3

Related Questions