[SOLVED] Astropy: determe if a point (RA, DEC) is inside a squared region given the coordinates of the corners of this region

Astropy: determe if a point (RA, DEC) is inside a squared region given the coordinates of the corners of this region

I have the coordinates of the 4 corners of a sky region obtained from the header fits image with astropy.wcs and the projection used to obtain the coordinates.

import astropy.wcs as wcs
import astropy.fits as fits
hdu=fits.open('filenmae.fits')
w=wcs.WCS(hdu[0].header)
print(w)

Which outputs:

[out]: WCS Keywords
Number of WCS axes: 2
CTYPE : 'RA---ZPN'  'DEC--ZPN'  
CRVAL : 308.45901  41.424847  
CRPIX : 6010.0186  -1881.9392  
CD1_1 CD1_2  : 1.0754576e-07  -5.5698074e-05  
CD2_1 CD2_2  : 5.5690351e-05  6.8120784e-08  
NAXIS : 4119  4119

From which I got the coordinates of the corners of the field with:

corners=w.calc_footprint()
print(corners)


[out]: [[308.318759    41.08966578]
 [308.01327548  41.08869031]
 [308.01293347  41.31890048]
 [308.31905451  41.31954629]]

How can I check if a random coordinate is inside the region defined by the coordinates of the coorners while respecting the effects of projection (as the output shows, the projection used was ZPN)?

Solution

Generally speaking, if you have a local coordinate system (say image pixels) that can be transformed via WCS to Sky coordinates, then one strategy would be to transform the "random coordinate" into the local pixel coordinates and do your region test in that flat rectilinear coordinate system.