...
 
Commits (6)
......@@ -87,7 +87,7 @@ def Moffat2D(fwhm, beta, shape, center=None, normalize=True):
alpha = fwhm / (2 * np.sqrt(2**(1 / beta) - 1))
amplitude = (beta - 1) * (np.pi * alpha**2)
if center is None:
x0, y0 = np.array(shape) / 2
x0, y0 = np.array(shape) / 2 - np.array([0.5,0.5])
else:
x0, y0 = center
xx, yy = np.mgrid[:shape[0], :shape[1]]
......
......@@ -184,5 +184,5 @@ def test_combine_fsf():
assert_allclose(fsf.get_beta(7000), fsf1.get_beta(7000), rtol=1.e-6)
fsf, cube = combine_fsf([fsf1, fsf2])
assert_allclose(fsf.get_fwhm(7000), 0.397959, rtol=1.e-3)
assert_allclose(fsf.get_beta(7000), 1.843269, rtol=1.e-3)
assert_allclose(fsf.get_fwhm(7000), 0.397959, rtol=1.e-2)
assert_allclose(fsf.get_beta(7000), 1.843269, rtol=1.e-2)
......@@ -725,10 +725,10 @@ class WCS:
np.clip(res, (0, 0), (self.naxis2 - 1, self.naxis1 - 1),
out=res)
return res
def coord(self, spaxel=False, relative=False, center=None,
def coord(self, spaxel=False, relative=False, center=None,
polar=False, unit=None, reshape=False, mask=None):
"""Return the full coordinate array of the image
"""Return the full coordinate array of the image
Parameters
----------
......@@ -757,38 +757,43 @@ class WCS:
array of [dec,ra], [p,q], [delta_dec, delta_ra], [delta_p, delta_q],
[r, theta]
"""
x,y = np.indices((self.naxis2,self.naxis1))
if not spaxel:
coord = np.mgrid[:self.naxis2, :self.naxis1].reshape(2, -1).T
coord = self.pix2sky(coord)
x,y = coord.T
x,y = np.meshgrid(x,y)
"""
pixcrd = np.indices((self.naxis2, self.naxis1))
if spaxel:
x, y = pixcrd
else:
coord = self.pix2sky(pixcrd.reshape(2, -1).T).T
x, y = coord.reshape(pixcrd.shape)
if mask is not None:
x = x[~mask]
y = y[~mask]
if relative or polar:
if center is None:
if spaxel:
center = 0.5*np.array([self.naxis2-1,self.naxis1-1])
center = 0.5 * np.array([self.naxis2 - 1, self.naxis1 - 1])
else:
center = self.get_center()
x0,y0 = center
center = self.get_center()
x0, y0 = center
x = x - x0
y = y - y0
if (not spaxel) and (unit is not None):
if not spaxel and unit is not None:
x = UnitArray(x, self.unit, unit)
y = UnitArray(y, self.unit, unit)
if polar:
theta = np.arctan2(y, x)
rho = np.hypot(x, y)
x = rho
y = theta
y = theta
if reshape:
x,y = np.vstack([x.ravel(), y.ravel()])
x = np.sort(np.unique(x))
y = np.sort(np.unique(y))
return x,y
x = np.unique(x)
y = np.unique(y)
return x, y
def pix2sky(self, x, unit=None):
"""Convert image pixel indexes (y,x) to world coordinates (dec,ra).
......@@ -1120,7 +1125,7 @@ class WCS:
[self.naxis2 - 1, self.naxis1 - 1]]
pixsky = self.pix2sky(pixcrd, unit=unit)
return np.hstack([pixsky.min(axis=0), pixsky.max(axis=0)])
def get_center(self, unit=None):
"""Return the center (dec,ra) of the image array.
......@@ -1138,8 +1143,8 @@ class WCS:
"""
crange = self.get_range(unit=unit)
center = 0.5*np.array([crange[0]+crange[2], crange[1]+crange[3]])
return center
center = 0.5 * np.array([crange[0] + crange[2], crange[1] + crange[3]])
return center
def get_start(self, unit=None):
"""Return the [dec,ra] coordinates of pixel (0,0).
......
......@@ -696,7 +696,7 @@ class CubeMosaic(CubeList):
def pycombine(self, nmax=2, nclip=5.0, var='propagate', nstop=2, nl=None,
header=None, mad=False):
crpix_out = self.wcs.wcs.wcs.crpix[::-1]
pos = np.array([crpix_out - cube.wcs.wcs.wcs.crpix[::-1]
pos = np.array([np.rint(crpix_out - cube.wcs.wcs.wcs.crpix[::-1])
for cube in self.cubes], dtype=int)
shapes = np.array([cube.shape[1:] for cube in self.cubes])
......
......@@ -142,42 +142,53 @@ class TestWCS:
assert wcs2.get_crval1() == 0.0
assert wcs2.get_crval2() == 0.0
wcs2.set_crval1(-2, unit=2*u.pix)
assert wcs2.get_crval1(unit=2*u.pix) == -2.0
wcs2.set_crval2(-2, unit=2*u.pix)
assert wcs2.get_crval2(unit=2*u.pix) == -2.0
wcs2.set_crval1(-2, unit=2 * u.pix)
assert wcs2.get_crval1(unit=2 * u.pix) == -2.0
wcs2.set_crval2(-2, unit=2 * u.pix)
assert wcs2.get_crval2(unit=2 * u.pix) == -2.0
def test_coord(self):
"""WCS class: testing coord"""
wcs = WCS(shape=(4,5), crval=(-23.0,5.0), cdelt=(0.2/3600.,-0.2/3600.),
wcs = WCS(shape=(4, 5),
crval=(-23.0, 5.0),
cdelt=(0.2 / 3600., -0.2 / 3600.),
deg=True)
p,q = wcs.coord(spaxel=True, relative=True)
assert (p.mean() == 0) & (q.mean() == 0)
p,q = wcs.coord(spaxel=True, relative=True, reshape=True)
assert (p.shape == (4,)) & (q.shape == (5,))
p,q = wcs.coord(spaxel=True, relative=True, center=(5.2,6.3))
assert (p.min() == -5.2) & (q.min() == -6.3)
dec,ra = wcs.coord()
p, q = wcs.coord(spaxel=True, relative=True)
assert p.mean() == 0
assert q.mean() == 0
assert p.shape == (4, 5)
assert q.shape == (4, 5)
p, q = wcs.coord(spaxel=True, relative=True, reshape=True)
assert p.shape == (4,)
assert q.shape == (5,)
p, q = wcs.coord(spaxel=True, relative=True, center=(5.2, 6.3))
assert p.min() == -5.2
assert q.min() == -6.3
assert p.shape == (4, 5)
assert q.shape == (4, 5)
dec, ra = wcs.coord()
assert_allclose(ra.min(), 4.99987929, rtol=1e-6)
assert_allclose(dec.min(), -23.000083333, rtol=1e-6)
dec,ra = wcs.coord(relative=True, unit='arcsec')
dec, ra = wcs.coord(relative=True, unit='arcsec')
assert_allclose(ra.min(), -0.4345444, rtol=1e-6)
assert_allclose(dec.min(),-0.3000001, rtol=1e-6)
r,theta = wcs.coord(polar=True, unit='arcsec')
assert_allclose(dec.min(), -0.3000001, rtol=1e-6)
r, theta = wcs.coord(polar=True, unit='arcsec')
assert_allclose(r.max(), 0.5280425, rtol=1e-6)
assert_allclose(theta.mean(),0.314159265, rtol=1e-6)
mask = np.empty(shape=(4,5), dtype=bool)
mask[:,:] = False
mask[0,0] = True
p,q = wcs.coord(spaxel=True, mask=mask)
assert (p.shape == (19,)) & (q.shape == (19,))
assert_allclose(theta.mean(), 0.314159265, rtol=1e-6)
mask = np.empty(shape=(4, 5), dtype=bool)
mask[:, :] = False
mask[0, 0] = True
p, q = wcs.coord(spaxel=True, mask=mask)
assert p.shape == (19,)
assert q.shape == (19,)
class TestWaveCoord:
......@@ -239,13 +250,13 @@ class TestWaveCoord:
assert wave.get_start(unit=u.nm) == 0.0
assert wave.get_end(unit=u.nm) == 9.0
wave.set_crval(2, unit=2*u.nm)
wave.set_crval(2, unit=2 * u.nm)
assert wave.get_crval() == 4.0
wave.set_crpix(2)
assert wave.get_crpix() == 2.0
wave.set_step(2)
assert wave.get_step() == 2.0
wave.set_step(2, unit=2*u.nm)
wave.set_step(2, unit=2 * u.nm)
assert wave.get_step() == 4.0
def test_rebin(self):
......