A body of black

that carried no reflection

defying its own room

un-earthly eggs of decreation

(In Flames)


Analysis of potentially hazardous asteroids' orbital parameters

Using data on already discovered asteroid's we'll try to analyse orbital parameters of near-Earth asteroids and find out familes of orbits with high probability of yet undiscovered potentially hazardous asteroids existance.


In [1]:
import pickle
from copy import copy, deepcopy
import numpy as np
import pandas as pd
import scipy.stats as ss
from sklearn import neighbors, svm
import matplotlib as mpl

# Import Asterion modules
import read_database as rdb
import generate_orbits as go
import learn_data as ld
import visualize_data as vd
import asterion_learn as al

# Plotting settings
%matplotlib inline
# font = {'size': 32}
# font = {'size': 25}
font = {'size': 14}
mpl.rc('font', **font)

Asteroid database

Load asteroid database:


In [2]:
# columns = ['a', 'e', 'i', 'w', 'om', 'q', 'H', 'neo', 
#                'pha', 'moid', 'per', 'n', 'ma', 'epoch']

columns = ['a', 'e', 'i', 'w', 'om', 'q', 'H', 'neo', 'pha', 'moid']
database = rdb.load_database(columns, jobtime=True)
database.head(n=5)


Loading asteroid database...
Asteroid database loaded in 6.207178 seconds.
Out[2]:
neo pha H e a q i om w moid
0 N N 3.34 0.075705 2.768134 2.558572 10.591705 80.314273 72.814757 1.59307
1 N N 4.13 0.230835 2.772779 2.132725 34.841162 173.087844 309.999215 1.23093
2 N N 5.33 0.256601 2.668831 1.984005 12.990102 169.862997 248.236797 1.03424
3 N N 3.20 0.089067 2.361348 2.151031 7.140406 103.842512 151.111635 1.13866
4 N N 6.85 0.191552 2.574343 2.081222 5.367954 141.591962 358.878896 1.09307

Extract NEAs


In [3]:
neas, num_neas = rdb.get_neas(database, columns)

Exclude extremal NEAs


In [4]:
print len(neas[neas.a > 4])
print len(neas[neas.i > 90])


12
2

In [5]:
neas = neas[neas.a < 4]
neas = neas[neas.i < 90]

Extract close-approaching NEAs


In [6]:
neas_close = neas[neas.q <= 1.08]

Extract groups of NEAs


In [7]:
# apollos, num_apollos = rdb.get_apollos(neas)
# atens, num_atens = rdb.get_atens(neas)
# amors, num_amors = rdb.get_amors(neas)

# print "num_neas:", num_neas
# print "apollos:", num_apollos, "ratio:", num_apollos/float(num_neas)
# print "atens:", num_atens, "ratio:", num_atens/float(num_neas)
# print "amors:", num_amors, "ratio:", num_amors/float(num_neas)

Extract PHAs and NHAs


In [8]:
haz_real, nohaz_real = rdb.get_hazMOID(neas_close)
print "Number of real PHAs", len(haz_real)
print "Number of real NHAs:", len(nohaz_real)


Number of real PHAs 6606
Number of real NHAs: 4476

In [9]:
# vd.display_allparams([neas_close, neas_close], vd.combs, vd.colnames)

Plot orbital distributions of rel PHAs and NHAs in the w-q plane


In [10]:
cutcol=['w', 'q']
vd.plot_distribution2d(cutcol, haz=haz_real, nohaz=nohaz_real, invertaxes=[0,1], labels=True)


/usr/lib/pymodules/python2.7/matplotlib/collections.py:548: FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison
  if self._edgecolors == 'face':

Generate virtual asteroids

Specify a full subset of orbital parameters


In [11]:
names = ['a', 'i', 'w', 'om', 'q']
full_names = [vd.colnames[nm] for nm in names]
data_subset = neas_close[names]

Fit distributions of asteroids' orbital parameters with continuous distributions


In [12]:
statdists = [ss.rayleigh, ss.lognorm, ss.uniform, go.HarmonicDistribution(), ss.johnsonsu]
distlist = go.get_param_distributions(data_subset, names, statdists, n=30, verbose=True)
go.plot_param_distributions(distlist, full_names, figsize=(16, 12))


Data cdf(xmax): 1.000000 	rayleigh_cdf(xmax): 0.999225
Data cdf(xmax): 1.000000 	lognorm_cdf(xmax): 0.982480
Data cdf(xmax): 1.000000 	uniform_cdf(xmax): 1.000000
Data cdf(xmax): 1.000000 	harmonic_cdf(xmax): 0.999970
Data cdf(xmax): 1.000000 	johnsonsu_cdf(xmax): 0.923356

Generate random asteroids on the basis of the fitted continuous distributions


In [13]:
rnum_str = '2e5'
rnum = int(float(rnum_str))
distdict = {name:dist for name, dist in zip(names, distlist)}
randdata = go.gen_orbits(distdict, num=rnum)

Check generated dataset to contain only asteroids with elliptical orbits


In [14]:
e_rand = randdata['e']
print "Total number of generated asteroids:", len(randdata)
print "Number of failed asteroids:", len(e_rand[e_rand < 0])
print "Number of asteroids with hyperbolic orbits:", len(e_rand[e_rand > 1])


Total number of generated asteroids: 200000
Number of failed asteroids: 0
Number of asteroids with hyperbolic orbits: 0

Calculate MOID for generated virtual asteroids or load previously computed results


In [15]:
haz_dump = './asteroid_data/neas_haz_gen_%s.p' % rnum_str
nohaz_dump = './asteroid_data/neas_nohaz_gen_%s.p' % rnum_str
try:
    haz_gen = rdb.loadObject(haz_dump)
    nohaz_gen = rdb.loadObject(nohaz_dump)
except:
    rdb.calc_moid(randdata, jobtime=True)
    haz_gen, nohaz_gen = rdb.get_hazMOID(randdata)
    
    rdb.dumpObject(haz_gen, haz_dump)
    rdb.dumpObject(nohaz_gen, nohaz_dump)

In [16]:
print "Number of generated non-uniform PHAs", len(haz_gen)
print "Number of generated non-uniform NHAs:", len(nohaz_gen)


Number of generated non-uniform PHAs 97817
Number of generated non-uniform NHAs: 102183

Generate additional set of uniformly distributed virtual asteroids


In [17]:
names = ['a', 'i', 'w', 'om', 'q']
data_subset = neas_close[names]
unum_str = '3e4'
unum = int(float(unum_str))
statdists_u = [ss.uniform]*len(names)

distlist_u = go.get_param_distributions(data_subset, names, statdists_u, n=30, verbose=True)
distdict_u = {name:dist for name, dist in zip(names, distlist_u)}
randdata_u = go.gen_orbits(distdict_u, num=unum)


Data cdf(xmax): 1.000000 	uniform_cdf(xmax): 1.000000
Data cdf(xmax): 1.000000 	uniform_cdf(xmax): 1.000000
Data cdf(xmax): 1.000000 	uniform_cdf(xmax): 1.000000
Data cdf(xmax): 1.000000 	uniform_cdf(xmax): 1.000000
Data cdf(xmax): 1.000000 	uniform_cdf(xmax): 1.000000

Check generated dataset to contain only asteroids with elliptical orbits


In [18]:
e_randu = randdata_u['e']
print "Total number of generated uniform asteroids:", len(randdata_u)
print "Number of failed asteroids:", len(e_randu[e_randu < 0])
print "Number of asteroids with hyperbolic orbits:", len(e_randu[e_randu > 1])


Total number of generated uniform asteroids: 30000
Number of failed asteroids: 0
Number of asteroids with hyperbolic orbits: 0

Calculate MOID for uniformly distributed virtual asteroids or load previously computed results


In [19]:
haz_dump = './asteroid_data/neas_haz_genu_%s.p' % unum_str
nohaz_dump = './asteroid_data/neas_nohaz_genu_%s.p' % unum_str
try:
    haz_genu = rdb.loadObject(haz_dump)
    nohaz_genu = rdb.loadObject(nohaz_dump)
except:
    rdb.calc_moid(randdata_u, jobtime=True)
    haz_genu, nohaz_genu = rdb.get_hazMOID(randdata_u)

    rdb.dumpObject(haz_genu, haz_dump)
    rdb.dumpObject(nohaz_genu, nohaz_dump)

In [20]:
print "Number of generated uniform PHAs", len(haz_genu)
print "Number of generated uniform NHAs:", len(nohaz_genu)


Number of generated uniform PHAs 7173
Number of generated uniform NHAs: 22827

Analyze virtual asteroids' orbital parameters


In [21]:
datasets_real = [haz_real, nohaz_real]
datasets_gen = [haz_gen, nohaz_gen]
datasets_genu = [haz_genu, nohaz_genu]

Investigate NEA distributions for pairs of orbital parameters


In [22]:
# vd.plot_alldistcombs(haz_real, nohaz_real)

Plot distributions in a w-q plane for the non-uniform generated virtual asteroids


In [23]:
cols = ['w', 'q']
vd.plot_distributions2d(cols, haz_gen, nohaz_gen, labels=True, invertaxes=[0,1])