GalSim icon indicating copy to clipboard operation
GalSim copied to clipboard

define `origin` and `world_origin` for CelestialWCS classes?

Open beckermr opened this issue 2 years ago • 8 comments

@sidneymau and I were working on some DES stuff and ran into an issue where we need the point on the sphere about which a given a WCS uses the projection. We found that the GSFitsWCS class defines a .center property (but that its value was (0, 0) which did not make much sense to us).

However, ignoring that .center was a bit funny in our case, as far as I can tell

  • .center is the point we want but is largely undocumented?
  • the origin and world_origin attributes are not defined

Is there a reason we couldn't set origin to CRPIX or similar and set world_origin to the tangent point on the sphere?

beckermr avatar Nov 07 '23 12:11 beckermr

Yes, center is that point, and yes it's largely undocumented. It's not a general WCS property for all CelestialWCS classes -- it's specific to the GSFitsWCS class. So you can't rely on it if you are using galsim.FitsWCS to generically read in a WCS from a FITS file (the return value isn't necessarily a GSFitsWCS object in that case). But if you are directly enforcing that you have a galsim.GSFitsWCS, then feel free to use it.

Is there a reason we couldn't set origin to CRPIX or similar and set world_origin to the tangent point on the sphere?

You certainly can do that. But neither of those has any requisite connection to the center of the tangent plane projection, so that's just a choice you can make.

rmjarvis avatar Nov 07 '23 14:11 rmjarvis

I must have made a typo when checking this last night! Here is what I see:

>>> coadd_wcs
galsim.GSFitsWCS(_data = ['TAN', array([5000.5, 5000.5]), array([[-7.305555555556e-05, 0.0], [0.0, 7.305555555556e-05]]), coord.CelestialCoord(coord.Angle(1.0591757458530349, coord.radians), coord.Angle(-0.43579901019234296, coord.radians)), None, None, None])
>>> coadd_wcs.center
coord.CelestialCoord(coord.Angle(1.0591757458530349, coord.radians), coord.Angle(-0.43579901019234296, coord.radians))
>>> coadd_wcs.origin
galsim.PositionD(x=0.0, y=0.0)
>>> coadd_wcs.world_origin
Traceback (most recent call last):
  File "<console>", line 1, in <module>
AttributeError: 'GSFitsWCS' object has no attribute 'world_origin'

so we do have .center and it is as expected. .origin is defined but not clear what it's use is. .world_origin is not defined.

sidneymau avatar Nov 07 '23 14:11 sidneymau

.cprix is also available, but maybe it's enough to know that it exists?

>>> coadd_wcs.crpix
array([5000.5, 5000.5])

sidneymau avatar Nov 07 '23 14:11 sidneymau

Oh, right. Yeah, world_origin isn't a thing for CelestialWCS objects. Just Euclidean ones. I forgot about that. I guess we could define world_origin to be wcs.toWorld(wcs.origin). But I'm not sure how useful that is, and it's not necessarily the center of the projection.

If you make it via TanWCS, then the world_origin parameter there does become the wcs.center attribute. But we don't persist that name as an attribute in the constructed wcs.

rmjarvis avatar Nov 07 '23 14:11 rmjarvis

Yeah origin is set to a dummy value in the code. It should probably raise since using the dummy value is almost always wrong and will give unexpected results.

beckermr avatar Nov 07 '23 14:11 beckermr

All the origin/world_origin stuff is really just to make sure we can recenter images (e.g.im = im.setCenter(0,0)) and have the WCS remain consistent in terms of what point on the sky each actual pixel refers to. GSFitsWCS manages this by updating crpix.

rmjarvis avatar Nov 07 '23 14:11 rmjarvis

Right yes. crval and crpix are the fits conventions for these things. I think I expected

  • to not have dummy properties that don't have real meaning attached
  • for there to be origin and world_origin attached to all WCS classes since these concepts exist for each kind of WCS no matter what

I agree that we cannot use the convention that the projection center for the tagent plane is world_origin for all celestial WCSs.

So maybe a small change in API would be to

  • either define origin / world_origin (could be crpix and center for GSFitsWCS) or have them raise useful errors
  • document center or similar as the attribute that holds the center of the projection if a Celestial WCS has a projection in use, otherwise it raises

beckermr avatar Nov 07 '23 14:11 beckermr

either define origin / world_origin (could be crpix and center for GSFitsWCS) or have them raise useful errors

I'm fine with that. I can't remember the reasoning behind the current implementation, but if someone want to track through the unit tests to see what breaks if origin = crpix, I don't have a problem with that. I think world_origin=center isn't something we can guarantee always though after possible calls to newOrigin. Could be wrong though.

document center or similar as the attribute that holds the center of the projection if a Celestial WCS has a projection in use, otherwise it raises

That's not particularly feasible for the other possible CelestialWCS classes. But you could do it for GSFitsWCS and RaDecFunction I think. Or if you want to try to figure out how to make a best-effort attempt for AstropyWCS and PyAstWCS, you could do that and raise an exception if it's too obfuscated to figure out. (Which I think will be the case for some wcs types.)

rmjarvis avatar Nov 07 '23 14:11 rmjarvis