Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
def fullsky_geometry(res=None, shape=None, dims=(), proj="car"):
"""Build an enmap covering the full sky, with the outermost pixel centers
at the poles and wrap-around points. Assumes a CAR (clenshaw curtis variant)
projection for now."""
assert proj == "car", "Only CAR fullsky geometry implemented"
if shape is None:
res = np.zeros(2)+res
shape = utils.nint(([1*np.pi,2*np.pi]/res) + (1,0))
ny, nx = shape
assert abs(res[0] * (ny-1) - np.pi) < 1e-8, "Vertical resolution does not evenly divide the sky; this is required for SHTs."
assert abs(res[1] * nx - 2*np.pi) < 1e-8, "Horizontal resolution does not evenly divide the sky; this is required for SHTs."
wcs = wcsutils.WCS(naxis=2)
# Note the reference point is shifted by half a pixel to keep
# the grid in bounds, from ra=180+cdelt/2 to ra=-180+cdelt/2.
wcs.wcs.crval = [res[0]/2/utils.degree,0]
wcs.wcs.cdelt = [-360./nx,180./(ny-1)]
wcs.wcs.crpix = [nx//2+0.5,ny//2+1]
wcs.wcs.ctype = ["RA---CAR","DEC--CAR"]
return dims+(ny,nx), wcs
"""Read an enmap from the specified hdf file. Two formats
are supported. The old enmap format, which simply used
a bounding box to specify the coordinates, and the new
format, which uses WCS properties. The latter is used if
available. With the old format, plate carree projection
is assumed. Note: some of the old files have a slightly
buggy wcs, which can result in 1-pixel errors."""
import h5py
with h5py.File(fname,"r") as hfile:
data = hfile["data"]
hwcs = hfile["wcs"]
header = astropy.io.fits.Header()
for key in hwcs:
header[key] = fix_python3(hwcs[key].value)
if wcs is None:
wcs = wcsutils.WCS(header).sub(2)
proxy = ndmap_proxy_hdf(data, wcs, fname=fname, threshold=sel_threshold)
return read_helper(proxy, sel=sel, box=box, pixbox=pixbox, geometry=geometry, wrap=wrap, mode=mode, delayed=delayed)
decrange = (decrange-np.mean(decrange))*1.1+np.mean(decrange)
decrange = np.array([max(-np.pi/2,decrange[0]),min(np.pi/2,decrange[1])])
decrange /= utils.degree
wdec = np.abs(decrange[1]-decrange[0])
# The shortest wavelength in the alm is about 2pi/lmax. We need at least
# two samples per mode.
res = 180./lmax/oversample
# Set up an intermediate coordinate system for the SHT. We will use
# CAR coordinates conformal on the equator, with a pixel on each pole.
# This will give it clenshaw curtis pixelization.
nx = utils.nint(360/res)
nytot = utils.nint(180/res)
# First set up the pixelization for the whole sky. Negative cdelt to
# make sharp extra happy. Not really necessary, but makes some things
# simpler later.
wcs = wcsutils.WCS(naxis=2)
wcs.wcs.ctype = ["RA---CAR","DEC--CAR"]
wcs.wcs.crval = [ra_ref,0]
wcs.wcs.cdelt = [360./nx,-180./nytot]
wcs.wcs.crpix = [nx/2.0+1,nytot/2.0+1]
# Then find the subset that includes the dec range we want
y1= utils.nint(wcs.wcs_world2pix(0,decrange[0],0)[1])
y2= utils.nint(wcs.wcs_world2pix(0,decrange[1],0)[1])
y1, y2 = min(y1,y2), max(y1,y2)
# Offset wcs to start at our target range
ny = y2-y1
wcs.wcs.crpix[1] -= y1
# Construct the map. +1 to put extra pixel at pole when we are fullsky
if verbose: print("Allocating shape %s dtype %s intermediate map" % (dims+(ny+1,nx),np.dtype(dtype).char))
tmap = enmap.zeros(dims+(ny+1,nx),wcs,dtype=dtype)
return tmap
def read_fits(fname, hdu=None, sel=None, box=None, pixbox=None, geometry=None, wrap="auto", mode=None, sel_threshold=10e6, wcs=None, delayed=False, verbose=False):
"""Read an enmap from the specified fits file. By default,
the map and coordinate system will be read from HDU 0. Use
the hdu argument to change this. The map must be stored as
a fits image. If sel is specified, it should be a slice
that will be applied to the image before reading. This avoids
reading more of the image than necessary. Instead of sel,
a coordinate box [[yfrom,xfrom],[yto,xto]] can be specified."""
if hdu is None: hdu = 0
hdu = astropy.io.fits.open(fname)[hdu]
ndim = len(hdu.shape)
if hdu.header["NAXIS"] < 2:
raise ValueError("%s is not an enmap (only %d axes)" % (fname, hdu.header["NAXIS"]))
if wcs is None:
with warnings.catch_warnings():
wcs = wcsutils.WCS(hdu.header).sub(2)
proxy = ndmap_proxy_fits(hdu, wcs, fname=fname, threshold=sel_threshold, verbose=verbose)
return read_helper(proxy, sel=sel, box=box, pixbox=pixbox, geometry=geometry, wrap=wrap, mode=mode, delayed=delayed)
def read_fits_geometry(fname, hdu=None):
"""Read an enmap wcs from the specified fits file. By default,
the map and coordinate system will be read from HDU 0. Use
the hdu argument to change this. The map must be stored as
a fits image."""
if hdu is None: hdu = 0
with utils.nowarn():
hdu = astropy.io.fits.open(fname)[hdu]
if hdu.header["NAXIS"] < 2:
raise ValueError("%s is not an enmap (only %d axes)" % (fname, hdu.header["NAXIS"]))
with warnings.catch_warnings():
wcs = wcsutils.WCS(hdu.header).sub(2)
shape = tuple([hdu.header["NAXIS%d"%(i+1)] for i in range(hdu.header["NAXIS"])[::-1]])
return shape, wcs