Search code examples
pythontrigonometryastropyastronomycartesian

Astronomy RA and DEC to cartesian coordinates then plot in a 2D image FOV


I would like to create an artificial image which simulates the capture of a camera https://www.theimagingsource.com/products/industrial-cameras/usb-2.0-monochrome/dmk41bu02h/ and the lens https://en.ids-imaging.com/store/lens-ricoh-fl-hc1214-2m-12-mm-1-2.html.

So my question is: having a star data from the Hipparcos catalog (https://cdsarc.u-strasbg.fr/viz-bin/ReadMe/I/239?format=html&tex=true) where I can get the RA and DEC, I would like to calculate its corresponding cartesian coordinates then somehow translate this [x,y,z] coordinates to [x,y] so I can draw it in the image.

And of course, the observer position should be on the earth as if I take a picture with the camera.

I am using python as the language to code it and opencv to create images, but any other language and/or library would be valid.

Thanks.


Solution

  • I was able to load the hip_main.dat catalog into a Table like so:

    >>> from astropy.table import Table
    >>> t = Table.read('hip_main.dat', format='cds', readme='http://cdsarc.u-strasbg.fr/ftp/cats/aliases/H/Hipparcos/ReadMe')
    

    (You should be able to put the URL for the filename as well, but for some reason that wasn't working for me, but it did work if I manually downloaded the file to my local directory first--the problem seems to be that the original file on the server is gzipped, and the CDS reader gets confused by the .gz in the filename. This would be a nice and easy thing to fix in Astropy.)

    It has coordinates in hms/dms as well as degrees. Degrees are faster to parse so then I made a SkyCoord from them. It automatically detects that the formats are in degrees:

    >>> coords = SkyCoord(ra=t['RAdeg'], dec=t['DEdeg'])                                                                                                                               
    >>> coords                                                                                                                                                                         
    <SkyCoord (ICRS): (ra, dec) in deg
        [(9.11850000e-04,   1.08901332), (3.79737000e-03, -19.49883745),
         (5.00795000e-03,  38.85928608), ..., (3.59976057e+02,   5.95663786),
         (3.59978239e+02, -64.3725722 ), (3.59978792e+02, -65.57707774)]>
    

    It also correctly detects ICRS J1991.25 from the ReadMe.

    You can get an Earth-bound observer position for example like:

    >>> from astropy.coordinates import EarthLocation                                                                                                                                  
    >>> location = EarthLocation.of_site('greenwich')                                                                                                                                  
    >>> location                                                                                                                                                                       
    <EarthLocation (3980608.90246817, -102.47522911, 4966861.27310068) m>
    

    (I'm just using a known location here but you can put coordinates anywhere on Earth).

    Then if you want you can create a local alt/az frame like:

    >>> from astropy.time import Time
    >>> from astropy.coordinates import AltAz                                                                                                                                          
    >>> local_frame = AltAz(location=location, obstime=Time.now())                                                                                                                     
    >>> local_frame                                                                                                                                                                    
    <AltAz Frame (obstime=2020-08-14 12:39:58.682416, location=(3980608.90246817, -102.47522911, 4966861.27310068) m, pressure=0.0 hPa, temperature=0.0 deg_C, relative_humidity=0.0, obswl=1.0 micron)>
    
    Then transform your original ICRS coordinates to this new frame:
    
    ```python
    >>> local_coords = coords.transform_to(local_frame)                                                                                                                                
    >>> local_coords                                                                                                                                                                   
    <SkyCoord (AltAz: obstime=2020-08-14 12:39:58.682416, location=(3980608.90246817, -102.47522911, 4966861.27310068) m, pressure=0.0 hPa, temperature=0.0 deg_C, relative_humidity=0.0, obswl=1.0 micron): (az, alt) in deg
        [(327.54653574, -32.6138386 ), (316.63406529, -51.59917121),
         (339.37842396,   3.4498418 ), ..., (329.41414414, -28.02088907),
         (217.27653631, -71.09315213), (214.14690597, -70.46865069)]>
    

    You can also convert this to cartesian coordinates with local_coords.cartesian.

    Unfortunately when it comes to projecting these onto some FoV of the sky for the purposes of image simulation, that's a bit beyond me. I could probably figure it out but probably someone else will more readily know how to do this. But at least now you have part of the equation (if I'm understanding your question correctly).