Skip to content

Commit

Permalink
Debug and change distance to 3 vectors.
Browse files Browse the repository at this point in the history
  • Loading branch information
@morriscb
morriscb committed Apr 14, 2020
1 parent c4a1638 commit b50fac0
Show file tree
Hide file tree
Showing 2 changed files with 71 additions and 45 deletions.
32 changes: 21 additions & 11 deletions python/lsst/pipe/tasks/functors.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -901,7 +901,7 @@ def _func(self, df):
return (df[self.colXX]*df[self.colYY] - df[self.colXY]**2)**0.25


class LocalWcs(Fuctor):
class LocalWcs(Functor):
"""Computations using the stored localWcs.
"""
name = "LocalWcsOperations"
Expand Down Expand Up @@ -956,6 +956,17 @@ def computeRaDec(self, x, y, cd11, cd12, cd21, cd22, ra, dec):
"""
return (x * cd11 + y * cd12 + ra, x * cd21 + y * cd22 + dec)

def computeVector(self, ra, dec):
"""
"""
vector = np.empty((len(ra), 3))
vector[:, 2] = np.sin(dec)

cosDec = np.cos(dec)
vector[:, 0] = np.cos(ra) * cosDec
vector[:, 1] = np.sin(ra) * cosDec
return vector

def computeSkySeperation(self, ra1, dec1, ra2, dec2):
"""Compute the local pixel scale conversion.
Expand All @@ -966,7 +977,7 @@ def computeSkySeperation(self, ra1, dec1, ra2, dec2):
dec1 : `pandas.Series`
Dec of the first coordinate in radians.
ra2 : `pandas.Series`
Ra of the second coordinate in radians.
RA of the second coordinate in radians.
dec2 : `pandas.Series`
Dec of the second coordinate in radians.
Expand All @@ -975,12 +986,9 @@ def computeSkySeperation(self, ra1, dec1, ra2, dec2):
dist : `pandas.Series`
Distance on the sphere in radians.
"""
deltaDec = dec2 - dec1
deltaRa = ra2 - ra1
return 2 * np.arcsin(
np.sqrt(
np.sin(deltaDec / 2) ** 2 +
np.cos(dec2) * np.cos(dec1) * np.sin(deltaRa / 2) ** 2))
vectors1 = self.computeVector(ra1, dec1)
vectors2 = self.computeVector(ra2, dec2)
return np.arccos(np.dot(vectors1, vectors2.transpose()))

def getSkySeperationFromPixel(self, x1, y1, x2, y2, cd11, cd12, cd21, cd22, ra, dec):
"""Compute the distance on the sphere from x2, y1 to x1, y1.
Expand Down Expand Up @@ -1062,9 +1070,11 @@ def pixelScaleArcseconds(self, x, y, cd11, cd12, cd21, cd22, ra, dec):
pixScale : `pandas.Series`
Arcseconds per pixel at the location of the local WC
"""
return 3600 * np.degrees(
self.getSkySeperationFromPixel(
x - 0.5, y, x + 0.5, y, cd11, cd12, cd21, cd22, ra, dec))
xdiff = self.getSkySeperationFromPixel(
x - 0.5, y, x + 0.5, y, cd11, cd12, cd21, cd22, ra, dec)
ydiff = self.getSkySeperationFromPixel(
x, y - 0.5, x, y + 0.5, cd11, cd12, cd21, cd22, ra, dec)
return 3600 * np.degrees((xdiff + ydiff) / 2)

def _func(self, df):
return self.pixelScaleArcseconds(df[self.colX],
Expand Down
84 changes: 50 additions & 34 deletions tests/test_functors.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -370,15 +370,17 @@ def testConvertPixelToArcseconds(self):
arcseconds.
"""
np.random.seed(1234)
refPoint = geom.Point(100, 100)
refPoint = geom.Point2D(100, 100)
for dec in np.linspace(-90, 90, 10):
for x, y in zip(np.random(1, 1024, 100), np.random(1, 1153, 10)):
for x, y in zip(np.random.uniform(1, 2 * 1109.99981456774, 10),
np.random.uniform(1, 2 * 560.018167811613, 10)):
point = geom.Point2D(x, y)
refPoint = geom.Point2D(x + 1, y + 1)
wcs = self._makeWcs(dec)
linAffMatrix = wcs.linearizePixelToSky(point, geom.degrees).getMatrix()
linAffMatrix = wcs.linearizePixelToSky(point, geom.radians).getMatrix()
self.dataDict["dipoleSep"] = np.full(self.nRecords, 10)
self.dataDict["slot_Centroid_x"] = np.full(self.nRecords, x)
self.dataDict["slot_Centroid_y"] = np.full(self.nRecords, x)
self.dataDict["slot_Centroid_y"] = np.full(self.nRecords, y)
self.dataDict["base_LocalWcs_CDMatrix_1_1"] = np.full(self.nRecords,
linAffMatrix[0, 0])
self.dataDict["base_LocalWcs_CDMatrix_1_2"] = np.full(self.nRecords,
Expand Down Expand Up @@ -416,51 +418,63 @@ def testConvertPixelToArcseconds(self):
df[("meas", "HSC-G", "slot_Centroid_y")],
df[("meas", "HSC-G", "base_LocalWcs_CDMatrix_1_1")],
df[("meas", "HSC-G", "base_LocalWcs_CDMatrix_1_2")],
df[("meas", "HSC-G", "base_LocalWcs_CDMatrix_1_1")],
df[("meas", "HSC-G", "base_LocalWcs_CDMatrix_1_2")],
df[("meas", "HSC-G", "base_LocalWcs_CDMatrix_2_1")],
df[("meas", "HSC-G", "base_LocalWcs_CDMatrix_2_2")],
df[("meas", "HSC-G", "base_LocalWcs_SkyOrigin_RA")],
df[("meas", "HSC-G", "base_LocalWcs_SkyOrigin_DEC")])
coord = wcs.pixelToSky(point)
refCoord = wcs.pixelToSky(refPoint)

# Test functor values against afw SkyWcs computations.
self.assertTrue(np.allclose(ra.values,
coord.getRa().asRadians(),
rtol=1e-8,
atol=1e-8))
rtol=0,
atol=1e-16))
self.assertTrue(np.allclose(dec.values,
coord.getDec().asRadians(),
rtol=1e-8,
atol=1e-8))

sep = func.computeSkySeperation(ra,
dec,
refCoord.getRa().asRadians(),
refCoord.getDec().asRadians())
self.assertTrue(np.allclose(sep.values,
coord.getDec().asRadians(),
rtol=1e-8,
atol=1e-8))
self.assertTrue(np.allclose(
coord.separation(refCoord).asRadians(),
sep.values,
atol=1e-8,
rtol=1e-8))
rtol=0,
atol=1e-16))
raRef, decRef = func.computeRaDec(
x + 1 / np.sqrt(2),
y + 1 / np.sqrt(2),
df[("meas", "HSC-G", "base_LocalWcs_CDMatrix_1_1")],
df[("meas", "HSC-G", "base_LocalWcs_CDMatrix_1_2")],
df[("meas", "HSC-G", "base_LocalWcs_CDMatrix_2_1")],
df[("meas", "HSC-G", "base_LocalWcs_CDMatrix_2_2")],
df[("meas", "HSC-G", "base_LocalWcs_SkyOrigin_RA")],
df[("meas", "HSC-G", "base_LocalWcs_SkyOrigin_DEC")])
self.assertTrue(np.allclose(raRef.values,
refCoord.getRa().asRadians(),
rtol=0,
atol=1e-2))
self.assertTrue(np.allclose(decRef.values,
refCoord.getDec().asRadians(),
rtol=0,
atol=1e-6))
sepRaDec = func.computeSkySeperation(ra,
dec,
raRef,
decRef)
"""
import pdb; pdb.set_trace()
self.assertTrue(np.allclose(sepRaDec,
coord.separation(refCoord).asRadians(),
rtol=0,
atol=1e-16))
sep = func.getSkySeperationFromPixel(df[("meas", "HSC-G", "slot_Centroid_x")],
df[("meas", "HSC-G", "slot_Centroid_y")],
refPoint.X(),
refPoint.Y(),
df[("meas", "HSC-G", "base_LocalWcs_CDMatrix_1_1")],
df[("meas", "HSC-G", "base_LocalWcs_CDMatrix_1_2")],
x + 1,
y + 1,
df[("meas", "HSC-G", "base_LocalWcs_CDMatrix_1_1")],
df[("meas", "HSC-G", "base_LocalWcs_CDMatrix_1_2")],
df[("meas", "HSC-G", "base_LocalWcs_CDMatrix_2_1")],
df[("meas", "HSC-G", "base_LocalWcs_CDMatrix_2_2")],
df[("meas", "HSC-G", "base_LocalWcs_SkyOrigin_RA")],
df[("meas", "HSC-G", "base_LocalWcs_SkyOrigin_DEC")])
self.assertTrue(np.allclose(
coord.separation(refCoord).asRadians(),
sep.values,
atol=1e-8,
rtol=1e-8))
rtol=0,
atol=1e-7))
func = ComputePixelScale("slot_Centroid_x",
"slot_Centroid_y",
Expand All @@ -471,11 +485,12 @@ def testConvertPixelToArcseconds(self):
"base_LocalWcs_SkyOrigin_RA",
"base_LocalWcs_SkyOrigin_DEC")
pixelScale = self._funcVal(func, parq)
import pdb; pdb.set_trace()
self.assertTrue(np.allclose(
wcs.pixelScale(point).asArcseconds(),
wcs.getPixelScale(point).asArcseconds(),
pixelScale.values,
atol=1e-8,
rtol=1e-8))
atol=1e-6,
rtol=0))
func = ConvertPixelToArcseconds("dipoleSep",
"slot_Centroid_x",
Expand All @@ -491,6 +506,7 @@ def testConvertPixelToArcseconds(self):
val.values,
atol=1e-13,
rtol=1e-13))
"""

def _makeWcs(self, dec=53.1595451514076):
"""Create a wcs from real CFHT values.
Expand Down

0 comments on commit b50fac0

Please sign in to comment.