GalSim-developers / GalSim

The modular galaxy image simulation toolkit. Documentation:
http://galsim-developers.github.io/GalSim/
Other
223 stars 105 forks source link

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

Open beckermr opened 9 months ago

beckermr commented 9 months ago

@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

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

rmjarvis commented 9 months ago

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.

sidneymau commented 9 months ago

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 commented 9 months ago

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

>>> coadd_wcs.crpix
array([5000.5, 5000.5])
rmjarvis commented 9 months ago

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.

beckermr commented 9 months ago

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.

rmjarvis commented 9 months ago

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.

beckermr commented 9 months ago

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

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

rmjarvis commented 9 months ago

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.)