<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Regrid from healpix
We show how to combine the healpix scheme with cygrid, using the healpy package. We will generate a simulated, full-sky data set and grid it onto a smaller fits image.
We start by adjusting the notebook settings.
Step1:  We attempt to limit our dependencies as much as possible, but astropy,  healpy, and wcsaxes needs to be available on your machine if you want to re-run the calculations. We can highly recommend anaconda as a scientific python platform.
Step2:  Let's define some plotting kwargs for the images.
Step3:  We start off by setting the basic healpix parameters which will define our coordinate system. It's given entirely by the nside. For more details, check the paper by Gorski et al. (2005). We stick to ring-ordering throughout this notebook.
Step4:  To generate a simulated fullsky map, we set up some power law and use healpy to create an image based on this.
Step5:  The random_field is easily generated by using hp.synfast.
Step6:  We manually create the header onto which we will grid the data. The target field is centered about $(l, b)=(180^\circ, 10^\circ)$. You need to make sure that all the celestial information are available.
Step7:  We start the gridding by initating the gridder with the header.
Step8:  Similar to the examples in previous notebooks, we need to define the coordinates for the input image. For healpix data, this is really straightforward and can easily be done with the healpy package.
Step9:  The gridding kernel is of key importance for the entire gridding process. cygrid allows you to specify the shape of the kernel (e.g. elliptical Gaussian or tapered sinc) and its size.
In our example, we use a symmetrical Gaussian (i.e. the major and minor axis of the kernel are identical). In that case, we need to furthermore specify kernelsize_sigma, the sphere_radius up to which the kernel will be computed, and the maximum acceptable healpix resolution for which we recommend kernelsize_sigma/2.
We refer to section 3.5 of the paper ('a minimal example') for a short discussion of the kernel parameters.
Step10:  After the kernel has been set, we perform the actual gridding by calling grid() with the coordinates and the data.
Step11:  To get the gridded data, we simply call get_datacube().`.
Step12:  And here is how our gridded field looks like.
Step13:  Finally, we can write the fits-image to disk. To make it two-dimensional, we have to create a celestial header. The WCS module can take care of this.
<ASSISTANT_TASK:>
Python Code:
%load_ext autoreload
%autoreload 2
%matplotlib inline
%config InlineBackend.figure_format = 'retina'
Explanation: Regrid from healpix
We show how to combine the healpix scheme with cygrid, using the healpy package. We will generate a simulated, full-sky data set and grid it onto a smaller fits image.
We start by adjusting the notebook settings.
End of explanation
import numpy as np
import matplotlib.pyplot as plt
import healpy as hp
from astropy.io import fits
from astropy.coordinates import SkyCoord
from astropy.wcs import WCS
import cygrid
Explanation: We attempt to limit our dependencies as much as possible, but astropy,  healpy, and wcsaxes needs to be available on your machine if you want to re-run the calculations. We can highly recommend anaconda as a scientific python platform.
End of explanation
imkw = dict(origin='lower', interpolation='nearest')
Explanation: Let's define some plotting kwargs for the images.
End of explanation
NSIDE = 256
NPIX = hp.nside2npix(NSIDE)
LMAX = 256
Explanation: We start off by setting the basic healpix parameters which will define our coordinate system. It's given entirely by the nside. For more details, check the paper by Gorski et al. (2005). We stick to ring-ordering throughout this notebook.
End of explanation
power = (lambda k: 12 / (k + 1) ** 2.5)
cl = power(np.arange(LMAX))
plt.loglog(cl)
plt.xlabel(r'$l$')
plt.ylabel(r'$C_l$')
Explanation: To generate a simulated fullsky map, we set up some power law and use healpy to create an image based on this.
End of explanation
random_field = hp.synfast(cl, NSIDE, verbose=False)
hp.mollview(random_field, xsize=300)
Explanation: The random_field is easily generated by using hp.synfast.
End of explanation
header = fits.Header()
pixsize = 14 / 60.
header['SIMPLE'] = 'T'
header['BITPIX'] = -32
header['NAXIS'] = 2
header['NAXIS1'] = 101
header['NAXIS2'] = 101
header['CDELT1'] = -pixsize
header['CDELT2'] = pixsize
header['CRPIX1'] = 51
header['CRPIX2'] = 51
header['CRVAL1'] = 180.
header['CRVAL2'] = 10.
header['LATPOLE'] = 90.
header['CTYPE1'] = 'GLON-SIN'
header['CTYPE2'] = 'GLAT-SIN'
wcs = WCS(header)
Explanation: We manually create the header onto which we will grid the data. The target field is centered about $(l, b)=(180^\circ, 10^\circ)$. You need to make sure that all the celestial information are available.
End of explanation
gridder = cygrid.WcsGrid(header)
Explanation: We start the gridding by initating the gridder with the header.
End of explanation
theta, phi = hp.pix2ang(NSIDE, np.arange(NPIX))
lons = np.rad2deg(phi).astype(np.float64)
lats = (90. - np.rad2deg(theta)).astype(np.float64)
Explanation: Similar to the examples in previous notebooks, we need to define the coordinates for the input image. For healpix data, this is really straightforward and can easily be done with the healpy package.
End of explanation
kernelsize_fwhm = 12. / 60.  # 12 arcminutes
# see https://en.wikipedia.org/wiki/Full_width_at_half_maximum
kernelsize_sigma = kernelsize_fwhm / np.sqrt(8 * np.log(2))
sphere_radius = 3. * kernelsize_sigma
gridder.set_kernel(
    'gauss1d',
    (kernelsize_sigma,),
    sphere_radius,
    kernelsize_sigma / 2.
    )
Explanation: The gridding kernel is of key importance for the entire gridding process. cygrid allows you to specify the shape of the kernel (e.g. elliptical Gaussian or tapered sinc) and its size.
In our example, we use a symmetrical Gaussian (i.e. the major and minor axis of the kernel are identical). In that case, we need to furthermore specify kernelsize_sigma, the sphere_radius up to which the kernel will be computed, and the maximum acceptable healpix resolution for which we recommend kernelsize_sigma/2.
We refer to section 3.5 of the paper ('a minimal example') for a short discussion of the kernel parameters.
End of explanation
gridder.grid(lons, lats, random_field)
Explanation: After the kernel has been set, we perform the actual gridding by calling grid() with the coordinates and the data.
End of explanation
datacube = gridder.get_datacube()
Explanation: To get the gridded data, we simply call get_datacube().`.
End of explanation
fig = plt.figure()
ax = fig.add_axes([0.1, 0.1, 0.8, 0.8], projection=wcs.celestial)
lon, lat = ax.coords
lon.set_axislabel('Galactic Longitude')
lat.set_axislabel('Galactic Latitude')
ax.imshow(datacube, **imkw)
ax.coords.grid(color='white', alpha=0.5, linestyle='solid')
Explanation: And here is how our gridded field looks like.
End of explanation
cel_header = wcs.celestial.to_header()
# fits.writeto('fullsky_regridded.fits', datacube, cel_header)
Explanation: Finally, we can write the fits-image to disk. To make it two-dimensional, we have to create a celestial header. The WCS module can take care of this.
End of explanation
<END_TASK>