boxes and mocpy pixels do not match

[samotracio]

This might be an issue with mocpy, but also from astropy.wcs or my own ignorance. Please excuse me if its related to the later.

I generated a rectangular region and investigated the positions of their flattened pixels. There seems to be an offset between these two.

from mocpy import MOC
import matplotlib.pyplot as plt
import numpy as np
from astropy.coordinates import Angle, Latitude, Longitude, SkyCoord
from astropy.table import Table
import regions

# define box
ra_center = Angle(131.0, 'deg')
dec_center = Angle(3.8, 'deg')
a = Angle(5.0, 'deg')
b = Angle(0.7, 'deg')
angle = Angle(90.0, 'deg')

box = regions.RectangleSkyRegion(SkyCoord(ra_center, dec_center, frame='icrs'), 2*a, 2*b, Angle(0, 'deg'))

# get moc from boxes
moc_box = MOC.from_box(lon=ra_center, lat=dec_center, a=a, b=b, angle=angle, max_depth=15)
hp_index = moc_box.flatten()

Now, the top-right corner of the box should be at (ra,dec)=(126, 4.5), but it is not (some wcs transformation?). Mocpy pixels seem to also follow this "shifted" box.

# plot moc and rectangle
fig = plt.figure(figsize=(7, 7))
wcs = moc_box.wcs(fig, projection='SIN')
ax = fig.add_subplot(projection=wcs)
moc_box.border(ax, wcs, color='r')
moc_box.fill(ax, wcs, alpha=0.2)
box.to_pixel(wcs).plot()
corner = [ra_center.value - a.value, dec_center.value + b.value]   # choose right-top corner
dwin = 0.1                                                         # half window size for plotting
ra_min, ra_max = corner[0]-dwin, corner[0]+dwin 
dec_min, dec_max = corner[1]-dwin, corner[1]+dwin
ax.set_xlim(wcs.wcs_world2pix(ra_max, dec_min, 1)[0], wcs.wcs_world2pix(ra_min, dec_max, 1)[0])
ax.set_ylim(wcs.wcs_world2pix(ra_max, dec_min, 1)[1], wcs.wcs_world2pix(ra_min, dec_max, 1)[1])
ax.coords['ra'].set_major_formatter('d.dddddd') ; ax.coords['dec'].set_major_formatter('d.dddddd')

Please correct me, but border pixels should be at order=15, so have a size of ~6.4", which is clearly not the case. Is there a way to increase the resolution while plotting?

After exporting the flattened list of pixels at max_depth and plotting using topcat, the shift is still present so there is little chance this could be an issue with matplotlib. Note, however that pixels are shown at the correct max_depth.

Converting the moc to a healsparse map and generating random points within the rectangle pixels suggests there are pixels outside the original box.

[ManonMarchand]

Thanks for opening the issue,

I see three different points so far:

1. Representation of the astropy region box

This is something that we do not have our hands on, but it looks like the conversion from world to pixels of the astropy-region box has some issues that are more or less apparent depending on the projection. See for example:

import matplotlib.pyplot as plt
import regions
from astropy import units as u
from astropy.coordinates import SkyCoord
from mocpy import MOC

center = SkyCoord(42, 43, unit="deg", frame="fk5")
box = regions.RectangleSkyRegion(center, 12 * u.deg, 6 * u.deg, angle)
moc_box = MOC.from_astropy_regions(box, max_depth=15)

fig = plt.figure(figsize=(5, 5))
wcs = moc_box.wcs(fig, projection="PAR")
ax = fig.add_subplot(projection=wcs)
moc_box.border(ax, wcs)
box.to_pixel(wcs).plot()

ax.set_ylim(120, 150)
ax.set_xlim(15, 45)

ax.coords['ra'].set_major_formatter('d.ddd') ; ax.coords['dec'].set_major_formatter('d.ddd')

plt.show()
plt.close()

fig = plt.figure(figsize=(5, 5))
wcs = moc_box.wcs(fig, projection="SIN")
ax = fig.add_subplot(projection=wcs)
moc_box.border(ax, wcs)
box.to_pixel(wcs).plot()

ax.set_ylim(120, 150)
ax.set_xlim(15, 45)

ax.coords['ra'].set_major_formatter('d.ddd') ; ax.coords['dec'].set_major_formatter('d.ddd')

Which highlights non conservation of the longitude of the corner between these two projections. It might be related to their conversion method here https://github.com/astropy/regions/blob/f807d23ec92e5d2b70af4f4febbfca378f822f72/regions/_utils/wcs_helpers.py#L13 that clearly states that the north/south axis should be good but not the other one. Maybe it is worth opening an issue in astropy-regions?

2. We silently degrade the MOC for low resolution plots.

You don't see order 15 cells in the plot because there was a choice in our library to degrade the MOC before plotting if the size of the healpix cell is smaller than a pixel defined in the WCS. Clearly, we did not think about svg exports or taking zoom ins of a WCS like you do. This is a bug on our side. Or a least we should add an option not to optimize the plot when requested.

3. The MOC does not contain the corners when defined with from_box

For this, I'll let @fxpineau confirm but I think it has something to do with box (conservation of width and height) versus quadrangle?

If you want a quadrangle (conservation of the corners) then you can use the method from_polygon

from astropy.coordinates import Longitude, Latitude

moc_polygon = MOC.from_polygon(
    Longitude([ra_center.value - a.value, ra_center.value - a.value,
               ra_center.value + a.value, ra_center.value + a.value], unit="deg"),
    Latitude([dec_center.value - b.value, dec_center.value + b.value,
              dec_center.value + b.value, dec_center.value - b.value], unit="deg"), max_depth=15)

[samotracio]

Thank you very much for your reply. I still can't find origin of the discrepancy and if its related to astropy-regions or mocpy. Perhaps is a matter of definition. As you well suggested the from_polygon(), fixes well the vertices and its sides are all great circles, as expected. The from_box() however is harder to follow. Below there is a plot of the pixels (order 13, so very dense) returned by mocpy in both cases (polygon in blue, box in red) for input coordinates (10,0), (10,20), (50,20), (50,0). From reading the docs, I would naively say a (non-rotated) RectangleSkyRegion() in astropy-regions should represent a rectangle in the sky, bounded by great circles and parallels, and from_box(), should retrieve its corresponding mocpy pixels. If this is not the way they work, is there a primitive that represents just that, a region bounded by great circles in ra and parallels in dec?

[fxpineau]

@samotracio If we understand correctly, you are looking for a Zone defined by two points of coordinates (RAmin, DEmin) and (RAmax, DEmax), in which edges of RA=cte follow great circle arcs while edges of DE=cte follow small circle arcs. We do have such a method in the Rust library, see here and we can definitely add it to MOCPy. Manon and FX

[samotracio]

Thanks @fxpineau and @ManonMarchand . That is indeed the case. I was trying to build a rectangle surrounding HSC's Hectomap field (high declination) and it certainly doesn't work with the primitives as they are. No experience with Rust, but from_zone() seems to do the trick. It would be great if you can add it to MOCPy. In any case, I am intrigued on why box() looks like this. Perhaps I am looking at the wrong documentation.

[fxpineau]

@samotracio What we call 'box' here (in MOCPy) is actually the polygon resulting from the intersection, on the unit-sphere, between two orthogonal spherical wedges having a same center. We may explicit this definition in the doc.

[fxpineau]

@ManonMarchand I just pushed from_zone in the Rust wrapper, see commit 916f8871c14b8600431cf8a0c336fdc1b9ac765a.