New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
DM-8828: Support proper motions in reference catalogs #111
Conversation
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
General comment while I'm in the middle of the review: should we be still writing doxygen for python docstrings, or should we just replace said doxygen with numpydoc whenever we add to/modify a docstring, since we're going to have to do that anyway.
dtype=float, | ||
doc="Coefficients of ordinary polynomial to convert epoch values to MJD TAI, " | ||
"from zeroth to highest order. To convert from Unix time, use [40587.0, 1.0/86400].", | ||
default=[0.0], |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Would it be better to default to Unix time? I'd expect that is the most common format.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I personally prefer MJD TAI. It is clear, unambiguous and more rigorously defined. But if the code uses astropy time objects then the issue goes away.
epoch_poly = pexConfig.ListField( | ||
dtype=float, | ||
doc="Coefficients of ordinary polynomial to convert epoch values to MJD TAI, " | ||
"from zeroth to highest order. To convert from Unix time, use [40587.0, 1.0/86400].", |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
It's not clear what this conversion is, or how one would figure out the correct polynomial values. Is it ever more than a linear equation? If not, we don't need a polynomial.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
As an alternative, please consider fields for "scale" and "format" as strings and then using them to construct astropy.time objects. This avoids the user having to work out polynomial coefficients and adds support for UTC, which is not linear.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
For the record: I am changing it as per my suggestion.
if len(self.mag_err_column_map) > 0 and not len(self.mag_column_list) == len(self.mag_err_column_map): | ||
raise ValueError("If magnitude errors are provided, all magnitudes must have an error column") | ||
if (sum(1 for col in ("pm_ra_name", "pm_dec_name", "epoch_name") if getattr(self, col) is not None) | ||
not in (0, 3)): |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
The expression in this if
does not parse well. Do we have a way to tie Fields together in pex_config? We probably should.
What about this, which makes the intent much clearer to me:
triplet = [x is None for x in (self.pm_ra_name, self.pm_dec_name, self.epoch_name)]
if not(all(triplet) or (not any(triplet))): # not (all True or all False)
raise
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I want to use this same functionality to validate the coord error and PM error fields as well, so I made a small function to handle it. In order to print a nice error message I found it best to accumulate a list of names of fields that are set, so the logic in that function is basically what's there now, rather than using all
and any
(although I admire the clarity of that approach). Have a look and see what you think.
not in (0, 3)): | ||
raise ValueError("Only all or none of pm_ra_name, pm_dec_name and epoch_name may be specified") | ||
if (self.pm_ra_err_name is not None) != (self.pm_dec_err_name is not None): # XOR | ||
raise ValueError("Only one of pm_ra_err_name and pm_dec_err_name has been specified") | ||
|
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Do you have with asssertRaises
tests for the three above cases? I think it would be worthwhile, and they'd be trivial to write.
""" | ||
if self.config.ra_err_name is None: # IngestIndexedReferenceConfig.validate ensures all or none | ||
record.set(key_map["coord_ra_err"], 0.0) | ||
record.set(key_map["coord_dec_err"], 0.0) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I think it would be better to set the errors to NaN
or None
, to make it very clear that there were no errors in the reference catalog.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I doubt None
is acceptable; this is a float entry. I mildly agree that nan
would be clearer than 0 for "unknown error", but 0 is not unreasonable. Does 0 simplify the math (i.e. avoid the need for special casing unknown error) or does it lead to division by zero?
|
||
|
||
class LoadPanstarrsObjectsTask(LoadIndexedReferenceObjectsTask): | ||
"""Reference catalog loader for ps1_pv3_3pi_20170110 catalog |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Period at the end.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Agreed, but this task went away
eventually adopted for the default of the LSST format: | ||
|
||
- The proper motions are in units of arcsec/yr instead of mas/yr. | ||
- The epoch is in seconds since the Unix epoch instead of MJD. |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I'm concerned about the creation of a ReferenceObjects class designed for exactly one catalog. I worry that we'll end up with similar LoadGaiaDR1FirstCutTask
, etc. We should probably at least call it LoadPanstarrs20170110ObjectsTask
or something, to make it clear that it's for one very specific PS catalog.
Could we instead write something to modify those catalogs to match the above format instead?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I'm fine with having a LoadPanstarrsObjectsTask
, but it seems like this is a good opportunity to write the accepted LSST one as well (although that can be a separate ticket).
In this future ticket, one might imagine making LoadPanstarrsObjectsTask
and LoadLSSTObjectsTask
(naming?) inherit from something slightly more general.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Please, everyone read the docstring, which explains the motivation behind this class: it only exists to provide backward compatibility.
LoadLsstObjectsTask
won't be necessary because we will create catalogs with the correct schema.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
@PaulPrice Ah, yes. I understand. Thank you for the clarification.
- ``units`` : units of proper motion (`lsst.afw.geom.Angle`) | ||
""" | ||
epochUnix = DateTime(float(epoch), DateTime.TAI).nsecs()/1.0e9 # Sec since Unix epoch | ||
timeDiff = (epochUnix - catalog["epoch"])/(3600*24*365.25) # Time difference, years |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This feels like an good use case for astropy.time
. Or do we have an appropriate time model in DateTime? 365.25
isn't technically precise, is it?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
3600 sec/hour *24 hours/day * 365.25 days is indeed defined to be what astronomers mean when they say "year" (where the sec is an SI second). And thus if an astronomer quotes to you a proper motion/year, this is the year they were supposed to have used.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I ended up using astropy.time
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Is that going to be as fast?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I think it will be fast enough. The time difference per object is an array operation (which is no guarantee, but it is bound to help).
@@ -105,6 +106,11 @@ class LoadReferenceObjectsConfig(pexConfig.Config): | |||
itemtype=str, | |||
default={}, | |||
) | |||
requireProperMotion = pexConfig.Field( | |||
doc="Require that the proper motion correction be successful?", |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
What does "be successful" mean in this context? The effect of this config parameter is not clear to me.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I agree. It means that the data needed to apply proper motion corrections is present. I'll try to improve the help string. How about:
"Require that the fields needed to correct proper motion (epoch, pm_ra and pm_dec) are present?"
raErr = "pmRaErr", | ||
decErr = "pmDecErr", | ||
epoch = "epoch", | ||
) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Why define this Struct if we've got a standard schema name for these?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
In addition to some cleanup comments, I have some general concerns about schema contents (allowing fields to be missing instead of filling them with NaN) and names (camelCase pmRa
vs. snake_case coord_ra
), and the use of things that really should have units but don't (i.e. all the various # arcseconds
comments and the 365.25
incorrect hard constants). It would be good to see whether we could fix the units here with astropy.Quantity
or better use of afw.Angles.
I'm also not sure of the need for a specific LoadPanstarrsObjectsTask
(see comment there).
I think I'd need to see this used in situ to better understand its behavior. I'll stare at the other tickets on this branch some and see if it becomes more clear.
refCat.sort() | ||
sourceCat.sort() | ||
return afwTable.unpackMatches(matchCat, refCat, sourceCat) | ||
|
||
def calculateProperMotionScales(self, catalog, epoch): | ||
"""Calculate the scaling factor(s) for proper motion |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Period at end of sentence.
The scaling factor is the number (or numbers, if the value | ||
varies per source) by which to multiply the proper motion | ||
rates in the catalog in order to get the proper motion. | ||
It means calculating the amount of time between the |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I'm not sure what "It means" is supposed to mean here.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I eliminated this -- all reference catalogs will use standard units
``epoch``, which is the mean epoch of the source as a | ||
Modified Julian Date (MJD) in TAI. | ||
epoch : `float` | ||
Epoch to which to move objects (TAI MJD). |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I don't believe TAI MJD
is an actual unit. @r-owen ?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
JD and in turn MJD are defined in times of TAI. One is supposed to explicitly note if one is using an JD or MJD calculated based on a different system.
https://www.iers.org/IERS/EN/Science/Recommendations/resolutionB1.html
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Exactly. But for the record, the epoch argument is now an astropy.time.Time
either a scalar `float` or an `numpy.ndarray`) | ||
- ``units`` : units of proper motion (`lsst.afw.geom.Angle`) | ||
""" | ||
timeDiff = (epoch - catalog[self._properMotionColumns.epoch])/365.25 # Time difference, years |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Here's another 365.25 that worries me (proliferation of random unnamed constants). It's not correct, and we should either use encode it correctly somewhere (daf.base.DateTime
?), or use astropy or equivalent.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
They are all gone
return pipeBase.Struct(values=timeDiff, units=1.0e-3*afwGeom.arcseconds) | ||
|
||
def applyProperMotions(self, catalog, epoch): | ||
"""Apply proper motion to the catalog |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
period at end of sentence.
tests/test_htmIndex.py
Outdated
default_config.epoch_name = "unixtime" | ||
default_config.pm_ra_err_name = "pm_ra_err" | ||
default_config.pm_dec_err_name = "pm_dec_err" | ||
default_config.pm_scale = 1.0e-3 # arcsec/yr --> mas/yr |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
There's a config
that gets created and thrown away in setUpClass
. Why not make it an instance variable and use it here?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I personally think this particular test is more readable by constructing the config each time, though I will change it to not re-use the same variable name for different configs.
tests/test_htmIndex.py
Outdated
self.assertFloatsEqual(references["coord_dec_err"], original["coord_dec_err"]) | ||
|
||
# One year difference --> coordinate difference should be self.proper_motion in RA only | ||
loader.applyProperMotions(references, self.mjd + 365.25) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Another 365.25 hard coded here...
tests/test_htmIndex.py
Outdated
self.assertFloatsAlmostEqual(references["coord_ra"], original["coord_ra"], rtol=1.0e-14) | ||
self.assertFloatsAlmostEqual(references["coord_dec"], original["coord_dec"], rtol=1.0e-14) | ||
self.assertFloatsEqual(references["coord_ra_err"], original["coord_ra_err"]) | ||
self.assertFloatsEqual(references["coord_dec_err"], original["coord_dec_err"]) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Good set of tests here (no correction means no change)!
tests/test_htmIndex.py
Outdated
@@ -70,26 +72,35 @@ def make_sky_catalog(out_path, size=1000): | |||
is_variable = np.random.randint(2, size=size) | |||
extra_col1 = np.random.normal(size=size) | |||
extra_col2 = np.random.normal(1000., 100., size=size) | |||
pm_ra = -np.ones(size)*cls.proper_motion # arcsec/year | |||
pm_dec = np.zeros(size) # arcsec/year |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Worth a comment here to make it clear that the intention is for the proper motion to be explicitly in ra only.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Good suggestion. Better yet I'll put the PM at an angle with cls.proper_motion_amt
and cls.proper_motion_dir
and make them instances of lsst.geom.Angle
tests/test_htmIndex.py
Outdated
self.assertFloatsAlmostEqual(orig.getCoord().bearingTo(ref.getCoord()).asDegrees(), | ||
180.0, rtol=1.0e-9) | ||
self.assertFloatsAlmostEqual(orig.getCoord().separation(ref.getCoord()).asArcseconds(), | ||
self.proper_motion, rtol=1.0e-12) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I'm somewhat surprised the tolerance has to be that high for both direction and distance.
I was about to say you should use assertSpherePointListAlmostEqual
, until I thought about the nature of the test (move X along RA only). So, maybe a comment about why the above test method wouldn't be ideal?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I believe this tolerance was based on writing out rather few digits to to the "original catalog" that is ingested. I have overhauled the code to compare angles and use absolute tolerance. Here is the new code:
self.assertAnglesAlmostEqual(orig.getCoord().separation(ref.getCoord()),
self.proper_motion_amt, maxDiff=1.0e-6*lsst.geom.arcseconds)
self.assertAnglesAlmostEqual(orig.getCoord().bearingTo(ref.getCoord()),
self.proper_motion_dir, maxDiff=1.0e-4*lsst.geom.arcseconds)
the separation is quite a bit better than this, but I didn't want to overdo it. Direction is harder to compute accurately, but the allowed error is quite acceptable.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Shouldn't be using under_scores
in variable names.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Good point. I'll fix it (though I don't think I made up those names; didn't you come up with them?)
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I don't think so. Probably they came in with the origin of the "LSST format catalogs", which I think was back when the policy was unclear and thought to be in transition.
dtype=float, | ||
doc="Coefficients of ordinary polynomial to convert epoch values to MJD TAI, " | ||
"from zeroth to highest order. To convert from Unix time, use [40587.0, 1.0/86400].", | ||
default=[0.0], |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I personally prefer MJD TAI. It is clear, unambiguous and more rigorously defined. But if the code uses astropy time objects then the issue goes away.
epoch_poly = pexConfig.ListField( | ||
dtype=float, | ||
doc="Coefficients of ordinary polynomial to convert epoch values to MJD TAI, " | ||
"from zeroth to highest order. To convert from Unix time, use [40587.0, 1.0/86400].", |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
As an alternative, please consider fields for "scale" and "format" as strings and then using them to construct astropy.time objects. This avoids the user having to work out polynomial coefficients and adds support for UTC, which is not linear.
""" | ||
if self.config.ra_err_name is None: # IngestIndexedReferenceConfig.validate ensures all or none | ||
record.set(key_map["coord_ra_err"], 0.0) | ||
record.set(key_map["coord_dec_err"], 0.0) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I doubt None
is acceptable; this is a float entry. I mildly agree that nan
would be clearer than 0 for "unknown error", but 0 is not unreasonable. Does 0 simplify the math (i.e. avoid the need for special casing unknown error) or does it lead to division by zero?
if self.config.pm_ra_err_name is None or self.config.pm_dec_err_name is None: | ||
return | ||
record.set(key_map["pmRaErr"], row[self.config.pm_ra_err_name]*self.config.pm_scale) | ||
record.set(key_map["pmDecErr"], row[self.config.pm_dec_err_name]*self.config.pm_scale) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Agreed; consistency here would be a big help. Reference catalogs use camelCase for field names so coordRaErr
or just raErr
would be better than coord_ra_err
.
@@ -326,6 +433,13 @@ def add_field(name): | |||
attr_name = 'is_{}_name'.format(flag) | |||
if getattr(self.config, attr_name): | |||
key_map[flag] = schema.addField(flag, 'Flag') | |||
if self.config.pm_ra_name is not None: # pm_dec_name and epoch_name also; by validation |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I agree. Clearly we need to support missing fields for ingest, but once ingested it's easier of all standard fields are always present, including proper motion, proper motion error and coordinate error. And I guess 0 makes more sense than nan
for for unknown proper motion because the math just works.
@@ -195,6 +209,7 @@ def loadPixelBox(self, bbox, wcs, filterName=None, calib=None): | |||
@param[in] wcs WCS (an lsst.afw.geom.SkyWcs) | |||
@param[in] filterName name of camera filter, or None or blank for the default filter | |||
@param[in] calib calibration, or None if unknown | |||
@param[in] epoch Epoch for proper motion correction (MJD TAI), or None |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Please accept epoch as a DateTime (here and in all other methods). This avoids avoid system or scale errors and is now the date is made available in VisitInfo
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I had not yet made this change. I have now done so.
9170ae2
to
0fa73d3
Compare
de7cd25
to
03742d1
Compare
Add support for proper motion and parallax (including error) to `LoadReferenceObjectsTask.makeMinimalSchema` and add coord error.
03742d1
to
06c1f1a
Compare
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Some (hopefully helpful) comment.
@@ -32,6 +35,9 @@ | |||
from lsst.afw.image import fluxFromABMag, fluxErrFromABMagErr | |||
from .indexerRegistry import IndexerRegistry | |||
from .readTextCatalogTask import ReadTextCatalogTask | |||
from .loadReferenceObjects import LoadReferenceObjectsTask | |||
|
|||
_RAD_PER_DEG = math.pi / 180 |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Why not just use np.radians
? Then you could get rid of your math
import.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Two reasons:
- I prefer to avoid numpy for scalars, for fear of inefficiencies and because it tends to confuse readers (who may think it is operating on array-like data)
- I personally dislike
math.radians
andnp.radians
because I feel their meaning is opaque. It implies that there are only two reasonable angular scales: degrees and radians, thus ignoring arcseconds, arcminutes, etc.
- Epoch is available in a column, in a format supported by astropy.time.Time. | ||
- Off-diagonal covariance terms are not available, e.g. | ||
covariance between RA and Dec, or between PM RA and PM Dec. | ||
Thus this is not an appropriate task to ingest Gaia catalogs. |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
It seems a little strange to call out Gaia as the one catalog we cannot ingest when, in fact, there must be many we cannot handle with this task. I get your point that it is the one we have on disk so someone might try.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I called out Gaia because it is the only astrometry catalog I know of that has covariance terms (certainly the other two we use -- PanSTARRS and SDSS -- do not). How about the following:
- There are no off-diagonal covariance terms, such as covariance
between RA and Dec, or between PM RA and PM Dec. Gaia is a well
known example of a catalog that has such terms, and thus should not
be ingested with this task.
@@ -169,6 +271,7 @@ def __init__(self, *args, **kwargs): | |||
self.indexer = IndexerRegistry[self.config.dataset_config.indexer.name]( | |||
self.config.dataset_config.indexer.active) | |||
self.makeSubtask('file_reader') | |||
self.radPerMilliarcsec = math.pi/(180*3600*1000) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
You have already defined _RAD_PER_DEG
, yet you are still calling math.pi
. I think I would still prefer math.radians
and math.degrees
.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
That's a great point. I will define a new constant next to _RAD_PER_DEG
, which isn't exactly what you wanted but is more self-consistent than the current code.
coordKey = catalog.table.getCoordKey() | ||
# Compute the offset of each object due to proper motion | ||
# as components of the arc of a great circle along RA and Dec | ||
offsetsRaRad = catalog["pm_ra"]*timeDiffsYears |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I'd slightly prefer that this be done in a utility library so others can use it if needed. Ideally it would delegate to a third party package to do the calculations, but if we don't have that, so be it.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I will file a ticket to overhaul this code to support covariance terms (which I suspect will not be found in any 3rd party package) and to support parallax.
tests/test_htmIndex.py
Outdated
@@ -232,16 +356,18 @@ def testIngest(self): | |||
config = LoadIndexedReferenceObjectsConfig() | |||
config.ref_dataset_name = "myrefcat" | |||
loader = LoadIndexedReferenceObjectsTask(butler=butler, config=config) | |||
cat = loader.loadSkyCircle(cent, self.search_radius, filterName='a') | |||
cat = loader.loadSkyCircle(cent, self.search_radius, filterName='a').refCat |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Good catch.
e1c8695
to
e40e369
Compare
@parejkoj in response to your main comments:
|
to reference catalog loaders. Known limitations: - It does not apply parallax - It only handles error in RA, Dec and proper motion in RA and Dec; it ignores other errors, including off-diagonal terms
791e2cd
to
1886e8b
Compare
1886e8b
to
e75e696
Compare
LoadIndexedReferenceObjectsTask.loadSkyCircle could return a catalog that was not contigous (despite code that attempted to guarantee otherwise). Add tests for this to test_htmIndex.py. Also document that the catalogs returned by LoadIndexedReferenceObjectsTask.loadSkyCircle and loadSkyBox are always contiguous.
54dd57b
to
af53c6f
Compare
by making the error different in RA and Dec
namely IngestIndexedReferenceTask
IngestIndexedReferenceTask. I also took the liberty of changing the data returned from its `run` method from a `Struct` containing 'result=None' to an empty `Struct`, since `run` did not, and still does not, return anything useful.
26bb0ca
to
dac5aba
Compare
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Couple of docstring comments.
Otherwise, I rescind my previous change request.
@param[in] row dict like object containing magnitude values | ||
@param[in] key_map Map of catalog keys to use in filling the record | ||
def _setMags(self, record, row, key_map): | ||
"""Set flux fields in a record of an indexed catalog. |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Is this the correct docstring for _setMags
?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Yes. I will rename the method to _setFlux
@param[in] radius radius of search region (an lsst.geom.Angle) | ||
@param[in] filterName name of filter, or None for the default filter; | ||
used for flux values in case we have flux limits (which are not yet implemented) | ||
@param[in] epoch Epoch for proper motion and parallax correction | ||
(an astropy.time.Time), or None |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Note something to the effect that if epoch is None, the corrections will not be applied.
"""!Load reference objects that overlap a circular sky region | ||
|
||
@param[in] ctrCoord ICRS center of search region (an lsst.geom.SpherePoint) | ||
@param[in] radius radius of search region (an lsst.geom.Angle) | ||
@param[in] filterName name of filter, or None for the default filter; | ||
used for flux values in case we have flux limits (which are not yet implemented) | ||
@param[in] epoch Epoch for proper motion and parallax correction | ||
(an astropy.time.Time), or None |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
same as my note about about loadSkyCircle()
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Scratch that: you can use docstring inheritance now, so you can delete the docstrings from the derived class's methods (see discussion on #dm-docs).
|
||
@staticmethod | ||
def make_sky_catalog(out_path, size=1000): | ||
@classmethod |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Good case for a mock
? But not on this ticket, probably.
leave the main task doc string alone for now, since I don't know what it is supposed to look like in Sphinx format.
because the units are different
No description provided.