Skip to content

Commit

Permalink
Chance localWcs unittest to Gnomonic wcs.
Browse files Browse the repository at this point in the history 
Several tests/changes were performed during the corse
of this ticket including different wcs approximations.
To see the full set of attempts see:
https://jira.lsstcorp.org/browse/DM-23946

Debug test fuctors

Test different declinations

Change dec range

Debug setting declination.

Add random xy locations.

Change LocalWcs functors and tests.

Debug and change distance to 3 vectors.

addapt functors and unittests.

Debug unittests.
  • Loading branch information
@morriscb
morriscb committed Apr 15, 2020
1 parent 13db8fb commit 3b9ae51
Show file tree
Hide file tree
Showing 2 changed files with 227 additions and 75 deletions.
153 changes: 131 additions & 22 deletions python/lsst/pipe/tasks/functors.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -901,10 +901,10 @@ def _func(self, df):
return (df[self.colXX]*df[self.colYY] - df[self.colXY]**2)**0.25


class ComputePixelScale(Functor):
"""Compute the local pixel scale from the stored CDMatrix.
class LocalWcs(Functor):
"""Computations using the stored localWcs.
"""
name = "Pixel Scale"
name = "LocalWcsOperations"

def __init__(self,
colCD_1_1,
Expand All @@ -918,50 +918,155 @@ def __init__(self,
self.colCD_2_2 = colCD_2_2
super().__init__(**kwargs)

def computeDeltaRaDec(self, x, y, cd11, cd12, cd21, cd22):
"""Compute the distance on the sphere from x2, y1 to x1, y1.
Parameters
----------
x : `pandas.Series`
X pixel coordinate.
y : `pandas.Series`
Y pixel coordinate.
cd11 : `pandas.Series`
[1, 1] element of the local Wcs affine transform.
cd11 : `pandas.Series`
[1, 1] element of the local Wcs affine transform.
cd12 : `pandas.Series`
[1, 2] element of the local Wcs affine transform.
cd21 : `pandas.Series`
[2, 1] element of the local Wcs affine transform.
cd22 : `pandas.Series`
[2, 2] element of the local Wcs affine transform.
Returns
-------
raDecTuple : tuple
RA and dec conversion of x and y given the local Wcs. Returned
units are in radians.
"""
return (x * cd11 + y * cd12, x * cd21 + y * cd22)

def computeSkySeperation(self, ra1, dec1, ra2, dec2):
"""Compute the local pixel scale conversion.
Parameters
----------
ra1 : `pandas.Series`
Ra of the first coordinate in radians.
dec1 : `pandas.Series`
Dec of the first coordinate in radians.
ra2 : `pandas.Series`
Ra of the second coordinate in radians.
dec2 : `pandas.Series`
Dec of the second coordinate in radians.
Returns
-------
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))

def getSkySeperationFromPixel(self, x1, y1, x2, y2, cd11, cd12, cd21, cd22):
"""Compute the distance on the sphere from x2, y1 to x1, y1.
Parameters
----------
x1 : `pandas.Series`
X pixel coordinate.
y1 : `pandas.Series`
Y pixel coordinate.
x2 : `pandas.Series`
X pixel coordinate.
y2 : `pandas.Series`
Y pixel coordinate.
cd11 : `pandas.Series`
[1, 1] element of the local Wcs affine transform.
cd11 : `pandas.Series`
[1, 1] element of the local Wcs affine transform.
cd12 : `pandas.Series`
[1, 2] element of the local Wcs affine transform.
cd21 : `pandas.Series`
[2, 1] element of the local Wcs affine transform.
cd22 : `pandas.Series`
[2, 2] element of the local Wcs affine transform.
Returns
-------
Distance : `pandas.Series`
Arcseconds per pixel at the location of the local WC
"""
ra1, dec1 = self.computeDeltaRaDec(x1, y1, cd11, cd12, cd21, cd22)
ra2, dec2 = self.computeDeltaRaDec(x2, y2, cd11, cd12, cd21, cd22)
# Great circle distance for small separations.
return self.computeSkySeperation(ra1, dec1, ra2, dec2)


class ComputePixelScale(LocalWcs):
"""Compute the local pixel scale from the stored CDMatrix.
"""
name = "PixelScale"

@property
def columns(self):
return [self.colCD_1_1, self.colCD_1_2,
self.colCD_2_1, self.colCD_2_2]
return [self.colCD_1_1,
self.colCD_1_2,
self.colCD_2_1,
self.colCD_2_2]

def pixelScale(self, cd11, cd12, cd21, cd22):
"""Compute the local pixel scale conversion.
def pixelScaleArcseconds(self, cd11, cd12, cd21, cd22):
"""Compute the local pixel to scale conversion in arcseconds.
Parameters
----------
cd11 : `pandas.Series`
[1, 1] element of the local CDMatricies.
[1, 1] element of the local Wcs affine transform in radians.
cd11 : `pandas.Series`
[1, 1] element of the local Wcs affine transform in radians.
cd12 : `pandas.Series`
[1, 2] element of the local CDMatricies.
[1, 2] element of the local Wcs affine transform in radians.
cd21 : `pandas.Series`
[2, 1] element of the local CDMatricies.
cd2 : `pandas.Series`
[2, 2] element of the local CDMatricies.
[2, 1] element of the local Wcs affine transform in radians.
cd22 : `pandas.Series`
[2, 2] element of the local Wcs affine transform in radians.
Returns
-------
pixScale : `pandas.Series`
Arcseconds per pixel at the location of the local WC
"""
return 3600 * np.sqrt(np.fabs(cd11 * cd22 - cd12 * cd21))
return 3600 * np.degrees(np.sqrt(np.fabs(cd11 * cd22 - cd12 * cd21)))

def _func(self, df):
return self.pixelScale(df[self.colCD_1_1], df[self.colCD_1_2],
df[self.colCD_2_1], df[self.colCD_2_2])
return self.pixelScaleArcseconds(df[self.colCD_1_1],
df[self.colCD_1_2],
df[self.colCD_2_1],
df[self.colCD_2_2])


class ConvertPixelToArcseconds(ComputePixelScale):
"""Convert a value in units pixels to units arcseconds.
"""
name = "Pixel scale converter"

def __init__(self,
col,
colCD_1_1,
colCD_1_2,
colCD_2_1,
colCD_2_2, **kwargs):
colCD_2_2,
**kwargs):
self.col = col
super().__init__(colCD_1_1, colCD_1_2, colCD_2_1, colCD_2_2, **kwargs)
super().__init__(colCD_1_1,
colCD_1_2,
colCD_2_1,
colCD_2_2,
**kwargs)

@property
def name(self):
Expand All @@ -970,12 +1075,16 @@ def name(self):
@property
def columns(self):
return [self.col,
self.colCD_1_1, self.colCD_1_2,
self.colCD_2_1, self.colCD_2_2]
self.colCD_1_1,
self.colCD_1_2,
self.colCD_2_1,
self.colCD_2_2]

def _func(self, df):
return df[self.col] * self.pixelScale(df[self.colCD_1_1], df[self.colCD_1_2],
df[self.colCD_2_1], df[self.colCD_2_2])
return df[self.col] * self.pixelScaleArcseconds(df[self.colCD_1_1],
df[self.colCD_1_2],
df[self.colCD_2_1],
df[self.colCD_2_2])


class ReferenceBand(Functor):
Expand Down
149 changes: 96 additions & 53 deletions tests/test_functors.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,6 +29,8 @@

import lsst.daf.base as dafBase
import lsst.afw.geom as afwGeom
import lsst.geom as geom
import lsst.meas.base as measBase
import lsst.utils.tests

# TODO: Remove skipUnless and this try block DM-22256
Expand All @@ -40,7 +42,7 @@
HsmTraceSize, PsfHsmTraceSizeDiff, HsmFwhm,
LocalPhotometry, LocalNanojansky, LocalNanojanskyErr,
LocalMagnitude, LocalMagnitudeErr,
ComputePixelScale, ConvertPixelToArcseconds)
LocalWcs, ComputePixelScale, ConvertPixelToArcseconds)
havePyArrow = True
except ImportError:
havePyArrow = False
Expand Down Expand Up @@ -368,57 +370,98 @@ def testConvertPixelToArcseconds(self):
"""Test calculations of the pixel scale and conversions of pixel to
arcseconds.
"""
wcs = self._makeWcs()
cdMatrix = wcs.getCdMatrix(wcs.getPixelOrigin())
self.dataDict["dipoleSep"] = np.full(self.nRecords, 10)
self.dataDict["base_LocalWcs_CDMatrix_1_1"] = np.full(self.nRecords,
cdMatrix[0, 0])
self.dataDict["base_LocalWcs_CDMatrix_1_2"] = np.full(self.nRecords,
cdMatrix[0, 1])
self.dataDict["base_LocalWcs_CDMatrix_2_1"] = np.full(self.nRecords,
cdMatrix[1, 0])
self.dataDict["base_LocalWcs_CDMatrix_2_2"] = np.full(self.nRecords,
cdMatrix[1, 1])
parq = self.simulateMultiParquet(self.dataDict)
func = ComputePixelScale("base_LocalWcs_CDMatrix_1_1",
"base_LocalWcs_CDMatrix_1_2",
"base_LocalWcs_CDMatrix_2_1",
"base_LocalWcs_CDMatrix_2_2")
df = parq.toDataFrame(columns={"dataset": "meas",
"filter": "HSC-G",
"columns": ["dipoleSep",
"base_LocalWcs_CDMatrix_1_1",
"base_LocalWcs_CDMatrix_1_2",
"base_LocalWcs_CDMatrix_2_1",
"base_LocalWcs_CDMatrix_2_2"]})
pixelScale = func.pixelScale(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")])
self.assertTrue(np.allclose(pixelScale.values,
wcs.getPixelScale().asArcseconds(),
rtol=0,
atol=1e-10))

# Test functors
val = self._funcVal(func, parq)
self.assertTrue(np.allclose(pixelScale.values,
val.values,
atol=1e-13,
rtol=0))

func = ConvertPixelToArcseconds("dipoleSep",
"base_LocalWcs_CDMatrix_1_1",
"base_LocalWcs_CDMatrix_1_2",
"base_LocalWcs_CDMatrix_2_1",
"base_LocalWcs_CDMatrix_2_2")
val = self._funcVal(func, parq)
self.assertTrue(np.allclose(pixelScale.values * 10,
val.values,
atol=1e-13,
rtol=0))

def _makeWcs(self):
dipoleSep = 10
np.random.seed(1234)
testPixelDeltas = np.random.uniform(-100, 100, size=(self.nRecords, 2))
import lsst.afw.table as afwTable
localWcsPlugin = measBase.EvaluateLocalWcsPlugin(
None,
"base_LocalWcs",
afwTable.SourceTable.makeMinimalSchema(),
None)
for dec in np.linspace(-90, 90, 10):
for x, y in zip(np.random.uniform(2 * 1109.99981456774, size=10),
np.random.uniform(2 * 560.018167811613, size=10)):

center = geom.Point2D(x, y)
wcs = self._makeWcs(dec)
skyOrigin = wcs.pixelToSky(center)

linAffMatrix = localWcsPlugin.makeLocalTransformMatrix(wcs,
center)
self.dataDict["dipoleSep"] = np.full(self.nRecords, dipoleSep)
self.dataDict["slot_Centroid_x"] = np.full(self.nRecords, x)
self.dataDict["slot_Centroid_y"] = np.full(self.nRecords, y)
self.dataDict["someCentroid_x"] = x + testPixelDeltas[:, 0]
self.dataDict["someCentroid_y"] = y + testPixelDeltas[:, 1]
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,
linAffMatrix[0, 1])
self.dataDict["base_LocalWcs_CDMatrix_2_1"] = np.full(self.nRecords,
linAffMatrix[1, 0])
self.dataDict["base_LocalWcs_CDMatrix_2_2"] = np.full(self.nRecords,
linAffMatrix[1, 1])
parq = self.simulateMultiParquet(self.dataDict)
func = LocalWcs("base_LocalWcs_CDMatrix_1_1",
"base_LocalWcs_CDMatrix_1_2",
"base_LocalWcs_CDMatrix_2_1",
"base_LocalWcs_CDMatrix_2_2")
df = parq.toDataFrame(columns={"dataset": "meas",
"filter": "HSC-G",
"columns": ["dipoleSep",
"slot_Centroid_x",
"slot_Centroid_y",
"someCentroid_x",
"someCentroid_y",
"base_LocalWcs_CDMatrix_1_1",
"base_LocalWcs_CDMatrix_1_2",
"base_LocalWcs_CDMatrix_2_1",
"base_LocalWcs_CDMatrix_2_2"]})

# Exercise the full set of functions in LocalWcs.
sepRadians = func.getSkySeperationFromPixel(
df[("meas", "HSC-G", "someCentroid_x")] - df[("meas", "HSC-G", "slot_Centroid_x")],
df[("meas", "HSC-G", "someCentroid_y")] - df[("meas", "HSC-G", "slot_Centroid_y")],
0.0,
0.0,
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")])

# Test functor values against afw SkyWcs computations.
for centX, centY, sep in zip(testPixelDeltas[:, 0],
testPixelDeltas[:, 1],
sepRadians.values):
afwSepRadians = skyOrigin.separation(
wcs.pixelToSky(x + centX, y + centY)).asRadians()
self.assertAlmostEqual(1 - sep / afwSepRadians, 0, places=6)

# Test the pixel scale computation.
func = ComputePixelScale("base_LocalWcs_CDMatrix_1_1",
"base_LocalWcs_CDMatrix_1_2",
"base_LocalWcs_CDMatrix_2_1",
"base_LocalWcs_CDMatrix_2_2")
pixelScale = self._funcVal(func, parq)
self.assertTrue(np.allclose(
wcs.getPixelScale(center).asArcseconds(),
pixelScale.values,
rtol=1e-8,
atol=0))

func = ConvertPixelToArcseconds("dipoleSep",
"base_LocalWcs_CDMatrix_1_1",
"base_LocalWcs_CDMatrix_1_2",
"base_LocalWcs_CDMatrix_2_1",
"base_LocalWcs_CDMatrix_2_2")
val = self._funcVal(func, parq)
self.assertTrue(np.allclose(pixelScale.values * dipoleSep,
val.values,
atol=1e-16,
rtol=1e-16))

def _makeWcs(self, dec=53.1595451514076):
"""Create a wcs from real CFHT values.
Returns
Expand All @@ -437,7 +480,7 @@ def _makeWcs(self):
metadata.set("EQUINOX", 2000.)

metadata.setDouble("CRVAL1", 215.604025685476)
metadata.setDouble("CRVAL2", 53.1595451514076)
metadata.setDouble("CRVAL2", dec)
metadata.setDouble("CRPIX1", 1109.99981456774)
metadata.setDouble("CRPIX2", 560.018167811613)
metadata.set("CTYPE1", 'RA---SIN')
Expand Down

0 comments on commit 3b9ae51

Please sign in to comment.