WCS same within deepcoadded tracts

Support
,
#1

Hi all,

Not sure if I am doing something wrong, but when I try to get the sky position of a pixel in my deep_coadd patches they are all the same within the same tract. I made a minimal code to show here…

target_ra = 53.2
target_dec = -28.1
radius = 0.52
region = sphgeom.Region.from_ivoa_pos(f"CIRCLE {target_ra} {target_dec} {radius}")

coadd_datasetrefs_i = butler.query_datasets("deep_coadd",
                                          where="patch.region OVERLAPS region AND band='i'",
                                          bind={"region": region},
                                          with_dimension_records=True,
                                          order_by=["patch"])

num = len(coadd_datasetrefs_i)

for inum in range(num):
    coadd = butler.get(coadd_datasetrefs_i[inum])
    patch = coadd.getMetadata()['LSST BUTLER DATAID PATCH']
    tract = coadd.getMetadata()['LSST BUTLER DATAID TRACT']
    wcs = coadd.getWcs()
    p1 = wcs.pixelToSky(0.0,0.0)   
    print(patch, tract, p1)

This produces…

1 5063 (54.0342443848, -28.3506114684)
2 5063 (54.0342443848, -28.3506114684)
3 5063 (54.0342443848, -28.3506114684)
4 5063 (54.0342443848, -28.3506114684)
5 5063 (54.0342443848, -28.3506114684)
#2

Ah I think I found the solution in a tutorial video: https://www.youtube.com/watch?v=OsjxN7PfiYk

The pixel has to be translated using the bounding box information, such as…

    i = 0.0
    j = 0.0
    bbox = coadd.getBBox()
    y = i + bbox.y.min
    x = j + bbox.x.min
    p1 = wcs.pixelToSky(x,y)  
    print(patch, tract, p1)

which now produces…

1 5063 (53.8576198474, -28.3628255733)
2 5063 (53.6682598394, -28.3637587078)
3 5063 (53.4788871444, -28.3644302488)
4 5063 (53.2895058991, -28.3648401400)

So probably working now I guess.

1 Like
#3

Hi @csabiu , thanks for your question. The origin of a deep_coadd is not (1,1); it depends on the patch. When you provide pixel (0,0), it will return the bottom left of the tract. Yes, your follow up test works.

1 Like
#4

For the record you can do this loop much simpler with:

for ref in coadd_datasetrefs_i:
    wcs = butler.get(ref.makeComponentRef("wcs"))
    bbox = butler.get(ref.makeComponentRef("bbox")) 
    patch = ref.dataId["patch"]
    tract = ref.dataId["tract"]
    ...

it’s also faster because it doesn’t have to read in all the coadds. The code you have takes 2.5 minutes because it has to download the files and read them all in but the component gets go straight to the server. My version takes 55s.

1 Like