Table of Contents
- 1 Useful Links
- 2 Imports
- 3 Retrieve information about Europa relating to our location and given time frame:
- 4 Format equatorial coordinates
- 5 Get center coordinates
- 6 Retrieve DSS survey images at FOV=1deg and FOV=11min
- 7 Image Data and WCS conversion
- 8 Image FOV and Overplot Asteroid Coordinates
- 9 Query Objects in FOV and Plot
Seminar 7 - Asteroids and Minor Planets¶
Useful Links¶
WCS and astronomic coordinate handling:
- https://docs.astropy.org/en/stable/wcs/
- http://docs.astropy.org/en/stable/visualization/wcsaxes/overlays.html
NASA JPS webpage to retrieve info about observable minor solar system objects: https://ssd.jpl.nasa.gov/sbwobs.cgi
MPC: Minor solar system object database: https://www.minorplanetcenter.net/data
astroplan
Imports¶
In [1]:
import numpy as np
from matplotlib import pylab as plt
import datetime as dt
### for minor solar system objects
import astroquery.jplhorizons as jp
import astroquery.skyview as SW
from astropy.wcs import WCS
from astropy import units as u
from astropy.coordinates import SkyCoord
from astroquery.simbad import Simbad
from astroquery.sdss import SDSS
Set observer coorinates for jp.Horizons
In [2]:
nhao = {"lon" : 134.33556,
"lat" : 35.02528,
"elevation" : 449}
In [3]:
nhao["lat"]
Out[3]:
Retrieve information about Europa relating to our location and given time frame:¶
In [4]:
Europa = jp.Horizons(id="Europa", location=nhao,
epochs={"start" : "2019-07-25",
"stop" : "2019-08-03",
"step" : "1h"})
In [5]:
Eur_data = Europa.ephemerides()
In [6]:
Eur_data
Out[6]:
In [7]:
Eur_data["datetime_str"][0]
Out[7]:
Format the time list to be accessable to matplotlib
In [8]:
dt.datetime.now()
Out[8]:
In [9]:
def get_time(data):
return dt.datetime.strptime(data,"%Y-%b-%d %H:%M")
In [10]:
TIME = []
for i in Eur_data["datetime_str"]:
TIME.append(get_time(i))
In [11]:
plt.figure(figsize = (12,8))
plt.plot(Eur_data["RA"][::10], Eur_data["DEC"][::10], "o")
for i in range(len(TIME))[::10]:
plt.text(Eur_data["RA"][i], Eur_data["DEC"][i], TIME[i])
plt.show()
In [12]:
Eur_data["RA"][0]
Out[12]:
Format equatorial coordinates¶
In [13]:
RAS_DECS = SkyCoord(Eur_data["RA"], Eur_data["DEC"], unit=(u.deg, u.deg),frame='icrs')
RAS = RAS_DECS.ra.deg
DECS= RAS_DECS.dec.deg
In [ ]:
Get center coordinates¶
In [14]:
RA_mid = RAS[RAS.shape[0]/2]
DEC_mid= DECS[RAS.shape[0]/2]
print "mid RA: ",RA_mid,"deg ---- mid DEC: ",DEC_mid, "deg"
Retrieve DSS survey images at FOV=1deg and FOV=11min¶
In [15]:
EU_im = SW.SkyView.get_images(position=str(RA_mid)+" "+str(DEC_mid),
radius = 1.*u.deg,
survey = ["DSS"])
EU_im2 = SW.SkyView.get_images(position=str(RAS[0])+" "+str(DECS[0]),
radius = 11./60.*u.deg,
survey = ["DSS"])
Image Data and WCS conversion¶
In [16]:
IM = EU_im[0][0].data
wcs = WCS(EU_im[0][0].header)
IM2 = EU_im2[0][0].data
wcs2 = WCS(EU_im2[0][0].header)
Image FOV and Overplot Asteroid Coordinates¶
In [17]:
plt.figure(figsize = (12,8))
### left image
ax = plt.subplot(121, projection = wcs)
ax.imshow(IM, cmap = "Greys")
ax.plot(Eur_data["RA"][::10], Eur_data["DEC"][::10], ".",
transform = ax.get_transform("world"), label = "Europa")
for i in range(len(TIME))[::50]:
ax.text(Eur_data["RA"][i]-0.02, Eur_data["DEC"][i]-0.03, TIME[i],
transform = ax.get_transform("world"))
### right image
ax2 = plt.subplot(122, projection = wcs2)
ax2.imshow(IM2,vmax = IM2.mean()+8.*IM2.std())
ax2.plot(Eur_data["RA"][::10][:3], Eur_data["DEC"][::10][:3], "--o",
transform = ax2.get_transform("world"), color = "w", label = "Europa")
ax2.text(Eur_data["RA"][0]+0.01, Eur_data["DEC"][0]-0.015, TIME[0],
transform = ax2.get_transform("world"), color = "w", fontsize = 16)
ax.legend()
ax2.legend()
ax.set_ylabel(r"DEC $[^{\circ}]$")
ax.set_xlabel(r"RA $[^{\circ}]$")
ax2.set_ylabel(r"DEC $[^{\circ}]$")
ax2.set_xlabel(r"RA $[^{\circ}]$")
ax.set_title( r"FOV : $1^{\circ}$" , loc="right")
ax2.set_title(r"FOV : $11^{\prime}$", loc="right")
plt.grid(True)
plt.show()
Query Objects in FOV and Plot¶
In [18]:
result_table = Simbad.query_region(str(RA_mid)+" "+str(DEC_mid), radius=0.7 * u.deg)
result_table
Out[18]:
In [19]:
result_table
TCOORD = RAS_DECS = SkyCoord(result_table["RA"], result_table["DEC"], unit=(u.hour, u.deg),frame='icrs')
SRA = TCOORD.ra.deg
SDEC= TCOORD.dec.deg
In [20]:
SRA[0]
Out[20]:
In [21]:
plt.figure(figsize = (5,4))
ax = plt.subplot(111, projection = wcs)
ax.imshow(IM, cmap = "Greys")
ax.scatter(SRA, SDEC, facecolors='none' ,
edgecolor="red",
transform = ax.get_transform("world"),
label ="Simbad Sources",
alpha = 0.4)
ax.plot(Eur_data["RA"][::10], Eur_data["DEC"][::10],
".",
transform = ax.get_transform("world"),
label = "Europa")
for i in range(len(TIME))[::50]:
ax.text(Eur_data["RA"][i]-0.02, Eur_data["DEC"][i]-0.03, TIME[i],
transform = ax.get_transform("world"))
ax.legend(loc = "lower right")
ax.set_ylabel(r"DEC $[^{\circ}]$")
ax.set_xlabel(r"RA $[^{\circ}]$")
ax.set_title(r"19/07/35 to 08/02 - FOV : $11^{\prime}$", loc="right")
plt.xlim([0,IM.shape[1]])
plt.ylim([0,IM.shape[0]*0.7])
plt.grid(True)
plt.show()
In [ ]:
In [ ]: