In [1]:
%run drop_test_simulator.ipynb
import pandas as pd
from bokeh.io import output_file, output_notebook, show
from bokeh.models import (
GMapPlot, GMapOptions, ColumnDataSource, Circle, LogColorMapper, BasicTicker, ColorBar,
PanTool, WheelZoomTool, BoxSelectTool, Range1d, Arrow, Segment, Ray
)
from bokeh.models.mappers import ColorMapper, LinearColorMapper
from bokeh.palettes import Viridis5, Inferno256
In [2]:
# linearization of globe based on our latitude
# I used http://www.csgnetwork.com/degreelenllavcalc.html and it seems to work out fine
def convert_to_lat(y):
lat_in_ft = 364618.1186706298
return y/lat_in_ft
def convert_to_long(x):
long_in_ft = 257742.06079152823
return x/long_in_ft
# convert from degrees on unit circle to degrees on a compass
def deg(theta):
theta = -(theta - 90) % 360
if theta == 0:
return str(360)
elif theta < 10:
return '00' + str(theta)
elif theta < 100:
return '0' + str(theta)
else:
return str(theta)
def polar_coord(df):
df.drop('norm', axis='columns', inplace=True)
df.drop('theta', axis='columns', inplace=True)
# magnitude of landing coordinate
df['norm'] = df.loc[:,['x','y']].apply(np.linalg.norm, axis=1, raw=True)
# azimuth angle of landing (note that this is on the unit circle)
df['theta'] = df.loc[:,['x','y']].apply(lambda pos: np.arctan2(pos[1], pos[0]), axis=1, raw=True)
# therefore, Y axis is latitude and X axis is longitude
def lat_long(df):
df.drop('latitude', axis='columns', inplace=True)
df.drop('longitude', axis='columns', inplace=True)
df['latitude'] = df.loc[:, 'y'].apply(lambda tmp: convert_to_lat(tmp))#, axis=1, raw=True)
df['longitude'] = df.loc[:, 'x'].apply(lambda tmp: convert_to_long(tmp))#, axis=1, raw=True)
# make an angry face, rawr! (mean impact point)
def get_mean(df):
mean_x = df.loc[:,['norm','theta']].apply(
lambda pos: pos[0] * np.cos(pos[1]), axis=1, raw=True).describe()['mean']
mean_y = df.loc[:,['norm','theta']].apply(
lambda pos: pos[0] * np.sin(pos[1]), axis=1, raw=True).describe()['mean']
return mean_x, mean_y
# slice off any data points where the wind is too strong
def wind_envelope(df, cutoff):
bool_list = df['v_w'] <= cutoff
df = df[bool_list]
return df
# shift the whole data set as a rigid motion
def translate_cluster(df, offset, release):
array = df.loc[:, ['x', 'y']]
X = [x for x in array['x'].apply(lambda x: x + offset[0])]
Y = [y for y in array['y'].apply(lambda y: y + offset[1])]
df.drop('x', axis='columns', inplace=True)
df.drop('y', axis='columns', inplace=True)
df['x'] = X
df['y'] = Y
polar_coord(df)
lat_long(df)
return df, release + offset
# i don't want to talk about how ugly this code block is lol
# rotates the whole data set as a rigid motion
def rotate_cluster(df, theta, release_angle, offset):
theta = - theta # we're rotating clockwise!
array = df.loc[:, ['x', 'y']]
X = [x for x in array['x']]
Y = [y for y in array['y']]
new = [rotate(np.array([x, Y[i]]), np.radians(theta)) for i, x in enumerate(X)]
X = [x[0] for x in new]
Y = [y[1] for y in new]
df.drop('x', axis='columns', inplace=True)
df.drop('y', axis='columns', inplace=True)
df['x'] = X
df['y'] = Y
polar_coord(df)
lat_long(df)
return df, theta + release_angle, rotate(offset, np.radians(theta))
In [3]:
def visualize(results_df, rls_pt, rls_ngl):
# these are not the means of the coordinates, taken individually
# this is the coordinate of the mean impact point, hence the polar coordinates
# technically i don't think there's a difference
mean_x, mean_y = get_mean(results_df)
print('Mean wind from ' + deg(rls_ngl % 360) + ' degrees true')
print('Release point: (' + str(rls_pt[0]) + ' ft, ' + str(rls_pt[1]) + ' ft)')
print('Release heading: ' + deg(rls_ngl) + ' degrees true')
print('Mean Impact: (' + str(mean_x) + ' ft, ' + str(mean_y) + ' ft)')
# assume that the bias is 0 until we have an actual target point
# so CEP = sqrt(MSE) = sqrt(var_norm + var_theta + cov_norm_theta)
# this is the mean distance from the mean impact point
variance = results_df.var(axis=0)[2:4].sum()
covariance = results_df.cov().iloc[2, 3]
# GPS is 95% inside a 2x2 m box, and 95% of values on one axis are within 2*RMS of ideal
# thus RMS = 1/2, since 2* 1/2 =1, and 1 on each side gives a 2 m bound. Therefore MSE = 1/4, which we double
# because there are two axes to account for, and also convert m^2 to ft^2
GPS_MSE = m_to_ft(m_to_ft(1/4))
bias = 0
total_MSE = variance + covariance + bias + 2*GPS_MSE
CEP = np.sqrt(total_MSE)
print('CEP:', CEP.round(3), 'ft')
plt.figure(1)
plt.hist2d(results_df['v_w'], results_df['deg_w'], 25)
plt.title('wind vectors')
plt.xlabel('v (knots)')
plt.ylabel('direction (degrees)')
plt.show()
plt.figure(2)
plt.hist(results_df['v_w'], 30)
plt.xlabel('wind magnitude (knots)')
plt.ylabel('count')
plt.title('wind speed distribution')
plt.show()
plt.figure(3)
plt.hist(results_df['norm'], 30)
plt.title('distances from origin')
plt.xlabel('radius (ft)')
plt.ylabel('count')
plt.show()
plt.figure(4)
plt.hist2d(results_df['x'],results_df['y'], 30)
plt.plot(mean_x, mean_y, 'r^')
plt.title('landing sites')
plt.xlabel('x (ft)')
plt.ylabel('y (ft)')
plt.show()
arrow_x = 200 * np.cos(np.radians(rls_ngl))
arrow_y = 200 * np.sin(np.radians(rls_ngl))
plt.figure(5)
results_df.plot.scatter('x', 'y', s=1, c='v_plane', colormap='inferno')
plt.axis('equal')
circle1 = plt.Circle((mean_x, mean_y), CEP, color='r', fill=False)
circle2 = plt.Circle((mean_x, mean_y), 2*CEP, color='tab:orange', fill=False)
circle3 = plt.Circle((mean_x, mean_y), 3*CEP, color='y', fill=False)
plt.gca().add_artist(circle1)
plt.gca().add_artist(circle2)
plt.gca().add_artist(circle3)
plt.title('impact points (and initial helicopter velocities)')
plt.arrow(rls_pt[0], rls_pt[1], arrow_x, arrow_y, width=.1, head_width=50, Fill=False)
plt.show()
plt.figure(6)
results_df.plot.scatter('x', 'y', s=1, c='v_w', colormap='inferno')
plt.axis('equal')
circle1 = plt.Circle((mean_x, mean_y), CEP, color='r', fill=False)
circle2 = plt.Circle((mean_x, mean_y), 2*CEP, color='tab:orange', fill=False)
circle3 = plt.Circle((mean_x, mean_y), 3*CEP, color='y', fill=False)
plt.gca().add_artist(circle1)
plt.gca().add_artist(circle2)
plt.gca().add_artist(circle3)
plt.title('impact points (and initial wind velocities)')
plt.arrow(rls_pt[0], rls_pt[1], arrow_x, arrow_y, width=.1, head_width=50, Fill=False)
plt.show()
return ((mean_x, mean_y), CEP)
# this requires a google api key. ask me on slack if you want to use mine, please keep it secret, keep it safe
def plot_on_map(df, release_pt, release_theta, CEP):
# release in the dead center of field unless shifted some
lat_0 = 45.208778
long_0 = -123.138722
mean_x, mean_y = get_mean(df)
release_long = convert_to_long(release_pt[0]) + long_0
release_lat = convert_to_lat(release_pt[1]) + lat_0
mean_long = convert_to_long(mean_x) + long_0
mean_lat = convert_to_lat(mean_y) + lat_0
print('Mean wind from ' + deg(release_theta % 360) + ' degrees true')
print("Release point: Lat " + str(release_lat) + ", Long " + str(release_long))
print("Release angle: Fly heading " + deg(release_theta % 360) + ' degrees true')
print('Mean impact point: Lat ' + str(mean_lat) + ', Long ' + str(mean_long))
print('Circular Error Probable: ' + str(CEP) + ' ft')
arrow_len = (convert_to_long(200) + convert_to_lat(200)) / 2
arrow_long = arrow_len * np.cos(np.radians(release_theta)) + release_long
arrow_lat = arrow_len * np.sin(np.radians(release_theta)) + release_lat
CEP *= 0.3048 # convert from feet to meters
map_options = GMapOptions(lat=lat_0, lng=long_0, map_type="satellite", zoom=16)
plot = GMapPlot(x_range=Range1d(), y_range=Range1d(), map_options=map_options)
plot.title.text = "Impact Points"
plot.api_key = "AIzaSyCYsl3277rp6nR3BUUx2q3Z9R4fyI2J1qo"
circle1 = Circle(x=mean_long, y=mean_lat, radius=CEP, line_width=2, line_color='red', fill_alpha=0)
circle2 = Circle(x=mean_long, y=mean_lat, radius=CEP*2, line_width=2, line_color='green', fill_alpha=0)
circle3 = Circle(x=mean_long, y=mean_lat, radius=CEP*3, line_width=2, line_color='blue', fill_alpha=0)
plot.add_glyph(circle1)
plot.add_glyph(circle2)
plot.add_glyph(circle3)
source = ColumnDataSource(
data=dict(
lat=df.latitude.apply(lambda y: y + lat_0).tolist(),
lon=df.longitude.apply(lambda x: x + long_0).tolist(),
v_h=df.v_plane.tolist(),
v_w=df.v_w.tolist()))
color_mapper = LinearColorMapper(palette=Inferno256, low=min(df['v_plane']), high=max(df['v_plane']))
circle = Circle(x="lon", y="lat", size="v_w", fill_color={'field': 'v_h', 'transform': color_mapper},
fill_alpha=0.5, line_color=None)
plot.add_glyph(source, circle)
line = Segment(x0=release_long, y0=release_lat, x1=arrow_long, y1=arrow_lat,
line_color="red", line_width=2)
arrow1 = Ray(x=arrow_long, y=arrow_lat, length=20, angle=release_theta-180+30,
angle_units="deg", line_color="red", line_width=2)
arrow2 = Ray(x=arrow_long, y=arrow_lat, length=20, angle=release_theta-180-30,
angle_units="deg", line_color="red", line_width=2)
plot.add_glyph(line)
plot.add_glyph(arrow1)
plot.add_glyph(arrow2)
color_bar = ColorBar(color_mapper=color_mapper, ticker=BasicTicker(),
label_standoff=12, border_line_color=None, location=(0,0))
plot.add_layout(color_bar, 'right')
plot.add_tools(PanTool(), WheelZoomTool())
output_notebook()
show(plot)
In [4]:
# data was generated by monte carlo sims, see other document for details
data_S = pd.read_csv('./drop_test_sample_data/succ_data.csv')
data_S.drop('Unnamed: 0', axis='columns', inplace=True)
data_A = pd.read_csv('./drop_test_sample_data/data_A.csv')
data_A.drop('Unnamed: 0', axis='columns', inplace=True)
data_B = pd.read_csv('./drop_test_sample_data/data_B.csv')
data_B.drop('Unnamed: 0', axis='columns', inplace=True)
data_C = pd.read_csv('./drop_test_sample_data/data_C.csv')
data_C.drop('Unnamed: 0', axis='columns', inplace=True)
In [5]:
default_release = np.array([0, 0])
default_angle = 230 # on unit circle, NOT compass
rls_S = default_release.copy()
ngl_S = default_angle
rls_A = default_release.copy()
ngl_A = default_angle
rls_B = default_release.copy()
ngl_B = default_angle
rls_C = default_release.copy()
ngl_C = default_angle
In [6]:
# logic added so that we can use just this block to control each of our scenarios in tandem
max_wind = 8.
do_cut_off = False
shift = np.array([-200, -200]) # +X = East, +Y = North
do_translate = False
spin = 60 # clockwise rotation
do_rotate = False
In [7]:
if do_cut_off:
data_S = wind_envelope(data_S, max_wind)
if do_rotate:
data_S, ngl_S, rls_S = rotate_cluster(data_S, spin, ngl_S, rls_S)
if do_translate:
data_S, rls_S = translate_cluster(data_S, shift, rls_S)
pos_S, CEP_S = visualize(data_S, rls_S, ngl_S)
plot_on_map(data_S, rls_S, ngl_S, CEP_S)
In [8]:
if do_cut_off:
data_A = wind_envelope(data_A, max_wind)
if do_rotate:
data_A, ngl_A, rls_A = rotate_cluster(data_A, spin, ngl_A, rls_A)
if do_translate:
data_A, rls_A = translate_cluster(data_A, shift, rls_A)
pos_A, CEP_A = visualize(data_A, rls_A, ngl_A)
plot_on_map(data_A, rls_A, ngl_A, CEP_A)
In [9]:
if do_cut_off:
data_B = wind_envelope(data_B, max_wind)
if do_rotate:
data_B, ngl_B, rls_B = rotate_cluster(data_B, spin, ngl_B, rls_B)
if do_translate:
data_B, rls_B = translate_cluster(data_B, shift, rls_B)
pos_B, CEP_B = visualize(data_B, rls_B, ngl_B)
plot_on_map(data_B, rls_B, ngl_B, CEP_B)
In [10]:
if do_cut_off:
data_C = wind_envelope(data_C, max_wind)
if do_rotate:
data_C, ngl_C, rls_C = rotate_cluster(data_C, spin, ngl_C, rls_C)
if do_translate:
data_C, rls_C = translate_cluster(data_C, shift, rls_C)
pos_C, CEP_C = visualize(data_C, rls_C, ngl_C)
plot_on_map(data_C, rls_C, ngl_C, CEP_C)
In [11]:
aggregate = pd.concat(map(pd.read_csv, ['./drop_test_sample_data/succ_data.csv',
'./drop_test_sample_data/data_A.csv',
'./drop_test_sample_data/data_B.csv',
'./drop_test_sample_data/data_C.csv']))
aggregate.drop('Unnamed: 0', axis='columns', inplace=True)
rls_agg = default_release.copy()
ngl_agg = default_angle
In [12]:
if do_cut_off:
aggregate = wind_envelope(aggregate, max_wind)
if do_rotate:
aggregate, ngl_agg, rls_agg = rotate_cluster(aggregate, spin, ngl_agg, rls_agg)
if do_translate:
aggregate, rls_agg = translate_cluster(aggregate, shift, rls_agg)
pos_agg, CEP_agg = visualize(aggregate, rls_agg, ngl_agg)
plot_on_map(aggregate, rls_agg, ngl_agg, CEP_agg)
In [ ]: