astropy
clarification on RectangleSkyRegion?
Hi!
This is a documentation question. Are RectangleSkyRegion sides (and the cross from the middle) following great circles or small circles?
This would help with the conversion to MOCs, to be sure that we're following what is intended here. For now, we're doing great circles everywhere.
Linked to this MOCPy issue and the IVOA ADQL / STC definition of a Box (which is not totally clear to us).
It depends, and you can perhaps help refine the definitions.
These regions are defined by a center, a width, a height, and an angle. https://github.com/astropy/regions/blob/f807d23ec92e5d2b70af4f4febbfca378f822f72/regions/shapes/rectangle.py#L366
They have no other properties - there is no definition of how the corners are connected. That definition is left to the PixelSkyRegion, which requires a projection (a WCS) to be defined.
Internally, the "Sky" regions get converted to "Pixel" regions (using the WCS) for most operations. The conversion simply converts the central sky coord to pixel coords and the widths/heights are converted from angular units to pixel units using a single pixel scale near the central coord. The "Sky" regions aren't really meant to represent regions on the celestial sphere. That would involve much more complicated pixel <-> sky region conversions, especially if the WCS contains distortions.
The docs state:
Sky regions are regions that are defined using celestial coordinates. Please note they are not defined as regions on the celestial sphere, but rather are meant to represent shapes on an image. They simply use sky coordinates instead of pixel coordinates to define their position. The remaining shape parameters are converted to pixels using the pixel scale of the image.
Thanks for the clarification. Then you would say that this is intended behavior for RectangleSkyRegions ?
from astropy.wcs import WCS
import regions
import astropy.units as u
from astropy.coordinates import SkyCoord
# define a box
center = SkyCoord(42, 43, unit="deg", frame="fk5")
box = regions.RectangleSkyRegion(center, 12 * u.deg, 6 * u.deg)
# define two similar WCS but in two projections
wcs_par = WCS("""WCSAXES = 2 / Number of coordinate axes
CRPIX1 = 250.5 / Pixel coordinate of reference point
CRPIX2 = 250.5 / Pixel coordinate of reference point
CDELT1 = -0.026810473472933 / [deg] Coordinate increment at reference point
CDELT2 = 0.026810473472933 / [deg] Coordinate increment at reference point
CUNIT1 = 'deg' / Units of coordinate increment and value
CUNIT2 = 'deg' / Units of coordinate increment and value
CTYPE1 = 'RA---PAR' / Right ascension, parabolic projection
CTYPE2 = 'DEC--PAR' / Declination, parabolic projection
CRVAL1 = 42.003016835441 / [deg] Coordinate value at reference point
CRVAL2 = 42.999058395073 / [deg] Coordinate value at reference point
LONPOLE = 0.0 / [deg] Native longitude of celestial pole
LATPOLE = 47.000941604927 / [deg] Native latitude of celestial pole
MJDREF = 0.0 / [d] MJD of fiducial time
RADESYS = 'ICRS' / Equatorial coordinate system
""")
wcs_sin = WCS("""WCSAXES = 2 / Number of coordinate axes
CRPIX1 = 250.5 / Pixel coordinate of reference point
CRPIX2 = 250.5 / Pixel coordinate of reference point
CDELT1 = -0.026810473472933 / [deg] Coordinate increment at reference point
CDELT2 = 0.026810473472933 / [deg] Coordinate increment at reference point
CUNIT1 = 'deg' / Units of coordinate increment and value
CUNIT2 = 'deg' / Units of coordinate increment and value
CTYPE1 = 'RA---SIN' / Right ascension, orthographic/synthesis project
CTYPE2 = 'DEC--SIN' / Declination, orthographic/synthesis projection
CRVAL1 = 42.003016835441 / [deg] Coordinate value at reference point
CRVAL2 = 42.999058395073 / [deg] Coordinate value at reference point
LONPOLE = 180.0 / [deg] Native longitude of celestial pole
LATPOLE = 42.999058395073 / [deg] Native latitude of celestial pole
MJDREF = 0.0 / [d] MJD of fiducial time
RADESYS = 'ICRS' / Equatorial coordinate system
""")
# print the first corner in sky coordinates
print(wcs_sin.pixel_to_world(*box.to_pixel(wcs_sin).corners[0]))
print(wcs_par.pixel_to_world(*box.to_pixel(wcs_par).corners[0]))
<SkyCoord (ICRS): (ra, dec) in deg
(49.82463355, 39.71836782)>
<SkyCoord (ICRS): (ra, dec) in deg
(50.17568385, 39.69382584)>
I mean: a single RectangleSkyRegion does not represent the same sky region depending on the projection that is chosen?
While potentially confusing, I think that was intentional. Pinging @astrofrog who was the original architect.
The background for the current architecture is given in https://github.com/astropy/regions/issues/107 and this was implemented in https://github.com/astropy/regions/pull/132. The current architecture mirrors how DS9 defines regions, where regions defined in world coordinates are basically the same shape as the ones defined in pixel coordinates, it's just that there are use cases where defining them with world coordinates is nicer. A good example of this is that if you are doing aperture photometry, you might want to specify the position of the circles as celestial coordinates, but you still want them to be circles in pixel space to carry out photometry.
The plan was always to add back 'proper' sky region classes defined on the celestial sphere though, but I dropped the ball on this and never had time to do it. I think we should revive this idea though!
