[AstroPy] Object within region

Mon Jun 13 11:44:56 EDT 2022

Hi Miguel, 

Maybe another way to say it that the regions code only deals with polygons drawn on flat maps of the sky, not with polygons on the actual, spherical sky.  So to make a flat map (of the surface of the Earth or of the sky) you have to make some choices, most notably about where the center of the map is going to be. (The same is true in Stellarium, btw - to view the sky on your flat computer screen, you have to choose which point in the sky is going to be at the center of your current view.) 

Yet another way to say the same thing is that, while routines for both SkyRegion and PixelRegion exist, I think the underlying calculations are always done as a PixelRegion.  (I’m not 100% sure that’s true, but that seems consistent with how it’s described in the documentation.) 

> I'm not sure I fully get the point. I'm referring to absolute coordinates, not points (ie, countries) which I have to figure out their coordinates. I have a set of 7 coordinates and want to know if another one falls into the corresponding polygon. I understand your example, but I'm not resolving/transforming between celestial objects and coordinates.

It doesn’t matter whether or not you have a specific object in mind - a coordinate gives you a spot on the map.  Lat/long gives you a point on the Earth, whether there is a city/country there or not.  RA/Dec gives you a point on the sky, whether there is a star there or not.  And either one can be (and indeed, has to be) transformed to an x,y coordinate to give you a position on a flat map.

To get a WCS that is behaving more like what you expect, i.e. where the point you are considering is at or near the center of your map, I would try something like this: 

import astropy.wcs as wcs
from astropy.wcs import WCS
import numpy as np

from regions import PolygonSkyRegion
from astropy.coordinates import SkyCoord

# Create a new WCS with two axes: 
w = wcs.WCS(naxis=2)
# Use RA and Dec as coords
w.wcs.ctype = ["RA---AIT", "DEC--AIT"]
# Set the pixel size to be whatever you want;
# units are degrees, so 1/60 is 1 arcminute.
pix_size = 1/60.
w.wcs.cdelt = np.array([pix_size, pix_size])
# How many pixels in x and y, given that we
# want to span the whole sky:
width_x = 360/pix_size
width_y = 180/pix_size
# Set the key reference pixel to be at the center pixel
# of the image: 
w.wcs.crpix = [int(width_x/2), int(width_y/2)]
# Choose your coordinates of interest - replace this with
# the RA and Dec of the point you're interested in testing, so
# it will be at the center of the image: 
myra = 21 * 15 # degrees, not hours
mydec = 57 # degrees
# Put that as the RA and Dec of the WCS reference pixel
w.wcs.crval = [myra, mydec]

reg=PolygonSkyRegion(SkyCoord([SkyCoord("2h","86d"),SkyCoord("7h","40d"),SkyCoord("5h","8d"),SkyCoord("4h","-5d"),SkyCoord("2h","-12d"),SkyCoord("2h","38d"),SkyCoord("2h","86d")]))
reg.contains(SkyCoord("21h","57d"),w)
reg.contains(SkyCoord(myra, mydec, unit="deg"),w)

(prints “False”)

Note that if you want to test other points on the sky, you will need to change the ‘crval’ in the wcs object again.  You don’t have to recreate the whole wcs, just change that one value.  So a cycle that looks like this: 

myra = 21 * 15 # degrees, not hours
mydec = 57 # degrees
# Put that as the RA and Dec of the WCS reference pixel
w.wcs.crval = [myra, mydec]
reg.contains(SkyCoord(myra, mydec, unit="deg"),w)

with different values of myra and mydec would allow you to test a bunch of points again the same region. 

Eric

> On Jun 13, 2022, at 10:40 AM, Miguel Gutiérrez Páez <mgutierrez at gmail.com> wrote:
> 
> Hi Eric,
> 
> Thanks again for your help.
> I'm not sure I fully get the point. I'm referring to absolute coordinates, not points (ie, countries) which I have to figure out their coordinates. I have a set of 7 coordinates and want to know if another one falls into the corresponding polygon. I understand your example, but I'm not resolving/transforming between celestial objects and coordinates. In any case, I agree that wcs could be involved and affect the results. So, as suspected, it's mainly an usability issue from my side. But please, since no image is involved in this process, do you know how to define/use/create the wcs I need?
> 
> Miguel
> 
> El lun, 13 jun 2022 a las 15:48, Eric Jensen (<ejensen1 at swarthmore.edu <mailto:ejensen1 at swarthmore.edu>>) escribió:
> Hi Miguel, 
> 
> Thanks for the example code.  To see what’s going on with your example, it might be helpful to think about a corresponding example with a world map: 
> 
> If I draw a polygon with vertices in the United Kingdom, Spain, Turkey, and Poland, is any point in North America inside that polygon? 
> 
> Your immediate answer is probably “no”, but the real answer is, “it depends on the map.”  If the world map is projected so that Europe is at or near the center of the map, then the answer is no.  But if the map is projected so that the center is in the Pacific Ocean, then it’s possible to have the UK and Spain on the far right-hand side of the map, and Turkey and Poland on the far left.  So a polygon with those four vertices, *if each edge stays on the map*, will indeed included parts of North America.  
> 
> That’s why the WCS is needed - it tells the routine how your map (image) is oriented, i.e. what coordinates are in the center, and what coordinates are nearer the edges.  When you look at Stellarium, presumably you are orienting your view of the sky so that your point of interest is somewhere near the center of your view.   But the WCS that you’re giving to the regions routine probably isn’t centered in that same way, so it may be that the polygon has some edges that have to go the “long way” (equivalent to a line from Spain to Turkey that cuts across North America) in order to stay on the map. 
> 
> Hope this helps, 
> 
> Eric
> 
> 
>> On Jun 13, 2022, at 6:10 AM, Miguel Gutiérrez Páez <mgutierrez at gmail.com <mailto:mgutierrez at gmail.com>> wrote:
>> 
>> Hi Eric,
>> 
>> Thanks for the reply.
>> Yes, I saw that note, but since I'm only working with RA/DEC coordinates, I didn't understand why I needed a wcs. But ok, no problem if I can use the example one.
>> 
>> I've tried to simplify the code and this is the simplest one:
>> 
>> reg=PolygonSkyRegion(SkyCoord([SkyCoord("2h","86d"),SkyCoord("7h","40d"),SkyCoord("5h","8d"),SkyCoord("4h","-5d"),SkyCoord("2h","-12d"),SkyCoord("2h","38d"),SkyCoord("2h","86d")]))
>> reg.contains(SkyCoord("21h","57d"),wcs)
>> 
>> which returns array([ True]). Checking with stellarium, it should return false.
>> 
>> My guess is the problem is with the polygon definition, which is quite irregular. Removing 3th and 4th coordinates makes a similar polygon and it does work.
>> 
>> Regards
>> 
>> El dom, 12 jun 2022 a las 15:56, Eric Jensen (<ejensen1 at swarthmore.edu <mailto:ejensen1 at swarthmore.edu>>) escribió:
>> Hi Miguel,
>> 
>>> One of my questions is why the "contains" method requires a wcs. If my only inputs are skycoords in icrs, why is needed to perform wcs transformations? This is a concept that I don't fully understand.
>> 
>> 
>> I think this is answered elsewhere in the docs, here: https://astropy-regions.readthedocs.io/en/stable/getting_started.html <https://astropy-regions.readthedocs.io/en/stable/getting_started.html>.  There it says: “Sky regions are regions that are defined using celestial coordinates. Note that these are not defined as regions on the celestial sphere, but rather are meant to represent shapes on an image, but simply defined using celestial coordinates as opposed to pixel coordinates.” 
>> 
>> Given this definition (that regions - even sky regions - always exist in the context of some image), then the WCS is needed to be able to convert between sky and pixel coords for the associated image.  
>> 
>> For your case, where maybe you only care about regions on the celestial sphere, it seems like you should be able to make this work using the WCS from the sample image they use in the docs.  (Caveat: I’ve never used this, so there could be some subtleties I’m missing here.)
>> 
>> If you can provide some sample code, and what results you get vs. what you expect to get, we can go from there.  
>> 
>> Eric
>> 
>> 
>>> On Jun 12, 2022, at 7:18 AM, Miguel Gutiérrez Páez <mgutierrez at gmail.com <mailto:mgutierrez at gmail.com>> wrote:
>>> 
>>> 
>>> Hi all,
>>> 
>>> I'm trying to figure out how to determine if a particular object falls within a certain sky region.
>>> First, please note I'm a complete noob with astropy, so I guess I'm probably missing some concept.
>>> 
>>> I've read examples about regions and sky regions: https://astropy-regions.readthedocs.io/en/stable/contains.html <https://astropy-regions.readthedocs.io/en/stable/contains.html>
>>> 
>>> I've been trying some examples by myself but I get some wrong results. That is, the script reports an object falls within the region, but after checking with stellarium for example, clearly is not.
>>> 
>>> One of my questions is why the "contains" method requires a wcs. If my only inputs are skycoords in icrs, why is needed to perform wcs transformations? This is a concept that I don't fully understand.
>>> The wcs I'm passing to the contains method, is the one described in the link above, the mentioned example.
>>> 
>>> As stated, I'm pretty sure I'm failing due to some missing concepts. Glad to hear from your for any clarification.
>>> 
>>> Thanks in advance.
>>> 
>>> Regards
>>> 
>>> PS: if you with I can paste some lines of my code
> 
> 
> _______________________________________________
> AstroPy mailing list
> AstroPy at python.org <mailto:AstroPy at python.org>
> https://mail.python.org/mailman/listinfo/astropy <https://mail.python.org/mailman/listinfo/astropy>
> _______________________________________________
> AstroPy mailing list
> AstroPy at python.org
> https://mail.python.org/mailman/listinfo/astropy

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://mail.python.org/pipermail/astropy/attachments/20220613/a2e65fe1/attachment-0001.html>