skyfielders / python-skyfield

Elegant astronomy for Python
MIT License
1.43k stars 213 forks source link

Is there a way to check if a RA/Dec is within a specific interval? #247

Closed victordomingos closed 5 years ago

victordomingos commented 5 years ago

I am trying to filter a list of objects according to their coordinates (RA/Dec), depending if they are within the boundaries of an "imaginary rectangle" defined two RA/Dec positions. Is there an easy way to do it using Skyfield?

Right now, what I need is something along these lines:

def is_inside_window(obj, min_ra, min_dec, max_ra, max_dec):
    ra_angle = Angle(f'{ra}d')
    dec_angle = Angle(f'{dec}d')
    if ra_angle.is_within_bounds(lower=min_ra, upper=max_ra) \
            and dec_angle.is_within_bounds(lower=min_dec, max=max_dec):
        return True
    else:
        print(f'Out: {obj}')
        return False

I would prefer to stay with Skyfield, since it has the level of abstracion I want, instead of going back and forth between different packages and their respective types/classes implementations. Having said that, I believe that it could be done using astropy.coordinates.Angle.is_within_bounds, but again, I couldn't find a straightforward way to convert between the Skyfield and Astropy Angle types. Am I missing something?

brandon-rhodes commented 5 years ago

For very large collections of objects (stars), I tend to put them in a Pandas dataframe and then use its filtering features to quickly remove objects not in the region of the sky I'm going to plot. You can find a few examples here:

https://rhodesmill.org/skyfield/stars.html#filtering-the-star-catalog

For other objects, I would probably just use simple inequality comparison on an attribute like ra.degrees. I cannot quite understand your code — it uses ra and dec but I don't see where they come from? Could you make your example more complete so that it will run if I put it in a .py file? Thanks!

victordomingos commented 5 years ago

Actually it was meant to be read as pseudo-code. I have been exploring astronomy stuff as a compete outsider and I have a little project I am trying do that does not require a lot of precision.

Basically what I am trying to do is to get data from PyOngc (includes objects from NGC/Messier/IC catalogs) and generate a list of interesting stuff (including nebulas and galaxies) for amateur sky watchers, allowing to filter the catalog by a few criteria, and display coordinates for a given place/time. One of the features I wanted to include is the ability to only display objects that are within certain altitude and azimuth ranges. Like a virtual rectangle in the sky.

I have not explored pandas yet, but I am curious if it could be used here to access the pyongc data, or any other similar database that is not limited to stars.

Sorry if my explanation does not comply with the right terminology or any other error.

brandon-rhodes commented 5 years ago

Would this routine work, and eliminate the need for another library with its Angle object?

def is_inside_window(obj, min_ra, min_dec, max_ra, max_dec):
    if not (min_ra <= ra_angle <= max_ra):
        return False
    if not (min_dec <= dec_angle <= max_dec):
        return False
    return True
victordomingos commented 5 years ago

Well, that's what i intended to do, but while fiddling with different data types or classes (most data comes from PyONGC as strings), I am finding it a bit hard, especially when trying to test it with external data for comparison. I must confess I am still reading through the skyfield/astropy/pyongc docs along the way, so I could have missed a lot of important stuff.

Can Skyfield Angle objects be compared as in your example, just like that, with operators? That would be very nice. I think that was one of my first attempts, but I started getting errors like this:

TypeError: '<=' not supported between instances of 'Angle' and 'Angle'

At the time, I even did a little search in Skyfield's Angle class code, and it seemed to me that those operators were not implemented yet…

brandon-rhodes commented 5 years ago

No, to compare angles you would do something like a.degrees < b.degrees. I did think about adding dozens of numeric operations to the angle objects so they could pretend to be numbers, but decided against it; it would be many lines of code to simply obscure what's really going on, which is comparing the degrees. (Or radians, or whatever.)

victordomingos commented 5 years ago

Ok, I understand that. In my exploratory code experiments, I managed to get something to run, but I keep getting weird results while testing. It could very well be that I still don't understand coordinate systems as well as I should.

In this file that I am using for testing code snippets, I try to extract the Messier objects, with variables for max/min alt/az. I have set those variables to -90/90 and 0/360, respectively.

_test_angles.py.zip

I would expect to get all objects this way, but it does not happen. I am suspecting right now that there could be some issue with the way I am trying to get RA/DEC corresponding to those values, in lines 86-87.

brandon-rhodes commented 5 years ago

I'm not familiar enough with how that other library operates to predict what kind of objects its methods return. Are you using Python 2 or Python 3?

victordomingos commented 5 years ago

I am using Python 3. PyONGC’s listObjects() returns a list of objects of their Dso class, which represents a single row from OpenNGC database:

https://github.com/mattiaverga/PyOngc/wiki/Dso-class-reference

It includes .getRA() and .gerDec() methods, which return Right Ascension and Declination in J2000 Epoch in string format. It also has a .getCoords That outputs a numpy array instead.

brandon-rhodes commented 5 years ago

Good, I'm glad you're using a recent version; Python 2 would let people compare strings and numbers by accident and so I wondered if the comparison was failing for you because of that.

My approach is always to add print statements until I understand. If a comparison a < b isn't behaving, I try:

print((a, b, a < b))

which generally shows me what's going on.

brandon-rhodes commented 5 years ago

Pending a reply revealing what you learned when you tried some print() calls, I'm going to temporarily close this issue to clear it off of my dashboard for now. Simple re-open it once you have some further ideas to share. Thanks!