def cov_ellipse(cov, q=None, nsig=None, **kwargs):
""" Code is slightly modified, but essentially borrowed from:
https://stackoverflow.com/questions/18764814/make-contour-of-scatter
"""
if q is not None:
q = np.asarray(q)
elif nsig is not None:
q = 2 * norm.cdf(nsig) - 1
else:
raise ValueError('Either `q` or `nsig` should be specified')
r2 = chi2.ppf(q, 2)
val, vec = np.linalg.eigh(cov)
width, height = 2 * np.sqrt(val[:, None] * r2)
rotation = np.degrees(np.arctan2(*vec[::-1, 0]))
return width, height, rotation
python类degrees()的实例源码
def decimal_to_dms(decimal_value):
'''
This converts from decimal degrees to DD:MM:SS, returned as a tuple.
'''
if decimal_value < 0:
negative = True
dec_val = fabs(decimal_value)
else:
negative = False
dec_val = decimal_value
degrees = trunc(dec_val)
minutes_deg = dec_val - degrees
minutes_mm = minutes_deg * 60.0
minutes_out = trunc(minutes_mm)
seconds = (minutes_mm - minutes_out)*60.0
if negative:
degrees = degrees
return '-', degrees, minutes_out, seconds
else:
return '+', degrees, minutes_out, seconds
def plot_ellipse(ax, mu, sigma, color="k"):
"""
Based on
http://stackoverflow.com/questions/17952171/not-sure-how-to-fit-data-with-a-gaussian-python.
"""
# Compute eigenvalues and associated eigenvectors
vals, vecs = np.linalg.eigh(sigma)
# Compute "tilt" of ellipse using first eigenvector
x, y = vecs[:, 0]
theta = np.degrees(np.arctan2(y, x))
# Eigenvalues give length of ellipse along each eigenvector
w, h = 2 * np.sqrt(vals)
ax.tick_params(axis='both', which='major', labelsize=20)
ellipse = Ellipse(mu, w, h, theta, color=color) # color="k")
ellipse.set_clip_box(ax.bbox)
ellipse.set_alpha(0.2)
ax.add_artist(ellipse)
def plot_ellipse(ax, mu, sigma, color="b"):
"""
Based on
http://stackoverflow.com/questions/17952171/not-sure-how-to-fit-data-with-a-gaussian-python.
"""
# Compute eigenvalues and associated eigenvectors
vals, vecs = np.linalg.eigh(sigma)
# Compute "tilt" of ellipse using first eigenvector
x, y = vecs[:, 0]
theta = np.degrees(np.arctan2(y, x))
# Eigenvalues give length of ellipse along each eigenvector
w, h = 2 * np.sqrt(vals)
ellipse = Ellipse(mu, w, h, theta, color=color) # color="k")
ellipse.set_clip_box(ax.bbox)
ellipse.set_alpha(0.2)
ax.add_artist(ellipse)
def plot_ellipse(ax, mu, sigma, color="b"):
"""
Based on
http://stackoverflow.com/questions/17952171/not-sure-how-to-fit-data-with-a-gaussian-python.
"""
# Compute eigenvalues and associated eigenvectors
vals, vecs = np.linalg.eigh(sigma)
# Compute "tilt" of ellipse using first eigenvector
x, y = vecs[:, 0]
theta = np.degrees(np.arctan2(y, x))
# Eigenvalues give length of ellipse along each eigenvector
w, h = 2 * np.sqrt(vals)
ellipse = Ellipse(mu, w, h, theta, color=color) # color="k")
ellipse.set_clip_box(ax.bbox)
ellipse.set_alpha(0.2)
ax.add_artist(ellipse)
def angle_wrap(angle,radians=False):
'''
Wraps the input angle to 360.0 degrees.
if radians is True: input is assumed to be in radians, output is also in
radians
'''
if radians:
wrapped = angle % (2.0*PI)
if wrapped < 0.0:
wrapped = 2.0*PI + wrapped
else:
wrapped = angle % 360.0
if wrapped < 0.0:
wrapped = 360.0 + wrapped
return wrapped
def dms_to_decimal(sign, degrees, minutes, seconds):
'''
Converts from DD:MM:SS to a decimal value. Returns decimal degrees.
'''
dec_deg = fabs(degrees) + fabs(minutes)/60.0 + fabs(seconds)/3600.0
if sign == '-':
return -dec_deg
else:
return dec_deg
############################
## DISTANCE AND XMATCHING ##
############################
def ecliptic_longitude(hUTC, dayofyear, year):
""" Ecliptic longitude
Args:
hUTC: fractional hour (UTC time)
dayofyear (int):
year (int):
Returns:
(float) the ecliptic longitude (degrees)
Details:
World Meteorological Organization (2006).Guide to meteorological
instruments and methods of observation. Geneva, Switzerland.
"""
jd = julian_date(hUTC, dayofyear, year)
n = jd - 2451545
# mean longitude (deg)
L = numpy.mod(280.46 + 0.9856474 * n, 360)
# mean anomaly (deg)
g = numpy.mod(357.528 + 0.9856003 * n, 360)
return L + 1.915 * numpy.sin(numpy.radians(g)) + 0.02 * numpy.sin(
numpy.radians(2 * g))
def hour_angle(hUTC, dayofyear, year, longitude):
""" Sun hour angle
Args:
hUTC: fractional hour (UTC time)
dayofyear (int):
year (int):
longitude (float): the location longitude (degrees, east positive)
Returns:
(float) the hour angle (hour)
Details:
World Meteorological Organization (2006).Guide to meteorological
instruments and methods of observation. Geneva, Switzerland.
"""
jd = julian_date(hUTC, dayofyear, year)
n = jd - 2451545
gmst = numpy.mod(6.697375 + 0.0657098242 * n + hUTC, 24)
lmst = numpy.mod(gmst + longitude / 15., 24)
ra = right_ascension(hUTC, dayofyear, year)
ha = numpy.mod(lmst - ra / 15. + 12, 24) - 12
return ha
def sun_elevation(hUTC, dayofyear, year, latitude, longitude):
""" Sun elevation
Args:
hUTC: fractional hour (UTC time)
dayofyear (int):
year (int):
latitude (float): the location latitude (degrees)
longitude (float): the location longitude (degrees)
Returns:
(float) the sun elevation (degrees)
Details:
World Meteorological Organization (2006).Guide to meteorological
instruments and methods of observation. Geneva, Switzerland.
"""
dec = declination(hUTC, dayofyear, year)
lat = numpy.radians(latitude)
ha = numpy.radians(hour_angle(hUTC, dayofyear, year, longitude) * 15)
sinel = numpy.sin(dec) * numpy.sin(lat) + numpy.cos(dec) * numpy.cos(
lat) * numpy.cos(ha)
return numpy.degrees(numpy.arcsin(sinel))
def _vlines(lines, ctrs=None, lengths=None, vecs=None, angle_lo=20, angle_hi=160, ransac_options=RANSAC_OPTIONS):
ctrs = ctrs if ctrs is not None else lines.mean(1)
vecs = vecs if vecs is not None else lines[:, 1, :] - lines[:, 0, :]
lengths = lengths if lengths is not None else np.hypot(vecs[:, 0], vecs[:, 1])
angles = np.degrees(np.arccos(vecs[:, 0] / lengths))
points = np.column_stack([ctrs[:, 0], angles])
point_indices, = np.nonzero((angles > angle_lo) & (angles < angle_hi))
points = points[point_indices]
if len(points) > 2:
model_ransac = linear_model.RANSACRegressor(**ransac_options)
model_ransac.fit(points[:, 0].reshape(-1, 1), points[:, 1].reshape(-1, 1))
inlier_mask = model_ransac.inlier_mask_
valid_lines = lines[point_indices[inlier_mask], :, :]
else:
valid_lines = []
return valid_lines
def _hlines(lines, ctrs=None, lengths=None, vecs=None, angle_lo=20, angle_hi=160, ransac_options=RANSAC_OPTIONS):
ctrs = ctrs if ctrs is not None else lines.mean(1)
vecs = vecs if vecs is not None else lines[:, 1, :] - lines[:, 0, :]
lengths = lengths if lengths is not None else np.hypot(vecs[:, 0], vecs[:, 1])
angles = np.degrees(np.arccos(vecs[:, 1] / lengths))
points = np.column_stack([ctrs[:, 1], angles])
point_indices, = np.nonzero((angles > angle_lo) & (angles < angle_hi))
points = points[point_indices]
if len(points) > 2:
model_ransac = linear_model.RANSACRegressor(**ransac_options)
model_ransac.fit(points[:, 0].reshape(-1, 1), points[:, 1].reshape(-1, 1))
inlier_mask = model_ransac.inlier_mask_
valid_lines = lines[point_indices[inlier_mask], :, :]
else:
valid_lines = []
return valid_lines
def screw_axis(self):
""" The rotation, translation and screw axis from the dual quaternion. """
rotation = 2. * np.degrees(np.arccos(self.q_rot.w))
rotation = np.mod(rotation, 360.)
if (rotation > 1.e-12):
translation = -2. * self.q_dual.w / np.sin(rotation / 2. * np.pi / 180.)
screw_axis = self.q_rot.q[0:3] / np.sin(rotation / 2. * np.pi / 180.)
else:
translation = 2. * np.sqrt(np.sum(np.power(self.q_dual.q[0:3], 2.)))
if (translation > 1.e-12):
screw_axis = 2. * self.q_dual.q[0:3] / translation
else:
screw_axis = np.zeros(3)
# TODO(ntonci): Add axis point for completeness
return screw_axis, rotation, translation
def getProjectedAngleInXYPlane(self, z=0, ref_axis=[0,1], centre=[0,0], inDeg=True):
'''
Project the OA vector to z=z, calculate the XY position, construct a
2D vector from [centre] to this XY and measure the angle subtended by
this vector from [ref_axis] (clockwise).
'''
ref_axis = np.array(ref_axis)
centre = np.array(centre)
point_vector_from_fit_centre = np.array(self.getXY(z=z)) - centre
dotP = np.dot(ref_axis, point_vector_from_fit_centre)
crossP = np.cross(ref_axis, point_vector_from_fit_centre)
angle = np.arccos(dotP/(np.linalg.norm(ref_axis)*np.linalg.norm(point_vector_from_fit_centre)))
if np.sign(crossP) > 0:
angle = (np.pi-angle) + np.pi
if inDeg:
dir_v = self._eval_direction_vector()
return np.degrees(angle)
else:
return angle
def _box_vectors_to_lengths_angles(box_vectors):
unitcell_lengths = []
for basis in box_vectors:
unitcell_lengths.append(np.array([np.linalg.norm(frame_v) for frame_v in basis]))
unitcell_angles = []
for vs in box_vectors:
angles = np.array([np.degrees(
np.arccos(np.dot(vs[i], vs[j])/
(np.linalg.norm(vs[i]) * np.linalg.norm(vs[j]))))
for i, j in [(0,1), (1,2), (2,0)]])
unitcell_angles.append(angles)
unitcell_angles = np.array(unitcell_angles)
return unitcell_lengths, unitcell_angles
def angle_map(self):
'''Returns a map of the angle for each pixel (w.r.t. origin).
0 degrees is vertical, +90 degrees is right, -90 degrees is left.'''
if self.angle_map_data is not None:
return self.angle_map_data
x = (np.arange(self.width) - self.x0)
y = (np.arange(self.height) - self.y0)
X,Y = np.meshgrid(x,y)
#M = np.degrees(np.arctan2(Y, X))
# Note intentional inversion of the usual (x,y) convention.
# This is so that 0 degrees is vertical.
#M = np.degrees(np.arctan2(X, Y))
# TODO: Lookup some internal parameter to determine direction
# of normal. (This is what should befine the angle convention.)
M = np.degrees(np.arctan2(X, -Y))
self.angle_map_data = M
return self.angle_map_data
def plot_motion(motion_mat):
time = motion_mat[:,0]
plt.figure(figsize=(15,5))
plt.subplot(1,2,1)
plt.plot(time,motion_mat[:,1]* 1000,label='x')
plt.plot(time,motion_mat[:,2]* 1000,label='y')
plt.plot(time,motion_mat[:,3]* 1000,label='z')
plt.xlabel('Time / s')
plt.ylabel('Translation / mm')
plt.legend()
plt.subplot(1,2,2)
plt.plot(time,np.degrees(motion_mat[:,4]),label='x')
plt.plot(time,np.degrees(motion_mat[:,5]),label='y')
plt.plot(time,np.degrees(motion_mat[:,6]),label='z')
plt.ylabel('Rotations / degrees')
plt.xlabel('Time / s')
plt.legend()
plt.show()
def random_points(self, n=1):
"""
Generate uniformly distributed random points within the sky
(i.e., all sky; on an unit sphere).
Returns
-------
lon : float, or 1D `~numpy.ndarray`
Longitudes (Galactic/equatorial), [0, 360) [deg].
lat : float, or 1D `~numpy.ndarray`
Latitudes (Galactic/equatorial), [-90, 90] [deg].
"""
theta, phi = spherical_uniform(n)
lon = np.degrees(phi)
lat = 90.0 - np.degrees(theta)
return (lon, lat)
##########################################################################
def threshold_test(self):
mx_adj, my_adj, mz_adj = self.mag_adj()
m_normal = np.sqrt(np.square(mx_adj)+np.square(my_adj)+np.square(mz_adj))
heading = np.degrees(np.arctan2(mx_adj/m_normal, my_adj/m_normal))
heading_diff = np.diff(heading)
rotate_index = np.insert(np.where(np.absolute(heading_diff)>20.0), 0, 0)
plt.plot(heading_diff)
plt.show()
angle_lst = []
for i in range(rotate_index.size):
try:
angle_onestep = np.mean(heading[rotate_index[i]: rotate_index[i+1]])
angle_lst.append(angle_onestep)
except:
pass
print angle_lst
def raw_mag_heading(self):
mx = self.raw_data['mx'].astype(np.float32)
my = self.raw_data['my'].astype(np.float32)
mz = self.raw_data['mz'].astype(np.float32)
m_normal = np.sqrt(np.square(mx)+np.square(my)+np.square(mz))
heading = np.arctan2(mx/m_normal, my/m_normal)
roll = np.arctan2(my/m_normal, mz/m_normal)
pitch = np.arctan2(mx/m_normal, mz/m_normal)
plt.plot(np.degrees(heading), "red", label="heading")
#plt.plot(np.degrees(roll), "green", label="roll")
#plt.plot(np.degrees(pitch), "blue", label="pitch")
#plt.plot(m_normal, "yellow", label='m_normal')
plt.legend(loc='upper left')
plt.show()
def ra_distance(declination, ra_a, ra_b):
"""
This function ...
:param declination: dec in degrees
:param ra_a: ra in degrees
:param ra_b: ra in degrees
:return:
"""
cos_ra_distance = np.sin(np.radians(declination))**2 + np.cos(np.radians(declination))**2 * np.cos(np.radians(ra_b-ra_a))
if cos_ra_distance > 1.0 and np.isclose(cos_ra_distance, 1.0): cos_ra_distance = 1.0 # Avoid crashes of np.arcos
# Return ...
return np.degrees(np.arccos(cos_ra_distance))
# -----------------------------------------------------------------
def ra_distance(declination, ra_a, ra_b):
"""
This function ...
:param declination: dec in degrees
:param ra_a: ra in degrees
:param ra_b: ra in degrees
:return:
"""
cos_ra_distance = np.sin(np.radians(declination))**2 + np.cos(np.radians(declination))**2 * np.cos(np.radians(ra_b-ra_a))
if cos_ra_distance > 1.0 and np.isclose(cos_ra_distance, 1.0): cos_ra_distance = 1.0 # Avoid crashes of np.arcos
# Return ...
return np.degrees(np.arccos(cos_ra_distance))
# -----------------------------------------------------------------
def tiltFactor(self, midpointdepth=None,
printAvAngle=False):
'''
get tilt factor from inverse distance law
https://en.wikipedia.org/wiki/Inverse-square_law
'''
# TODO: can also be only def. with FOV, rot, tilt
beta2 = self.viewAngle(midpointdepth=midpointdepth)
try:
angles, vals = getattr(
emissivity_vs_angle, self.opts['material'])()
except AttributeError:
raise AttributeError("material[%s] is not in list of know materials: %s" % (
self.opts['material'], [o[0] for o in getmembers(emissivity_vs_angle)
if isfunction(o[1])]))
if printAvAngle:
avg_angle = beta2[self.foreground()].mean()
print('angle: %s DEG' % np.degrees(avg_angle))
# use averaged angle instead of beta2 to not overemphasize correction
normEmissivity = np.clip(
InterpolatedUnivariateSpline(
np.radians(angles), vals)(beta2), 0, 1)
return normEmissivity
def _moveruler(self, evt):
x, y = self.mouseCoord(evt)
txtPosX = (self.rulersStartX + x) * 0.5
txtPosY = (self.rulersStartY + y) * 0.5
dx = x - self.rulersStartX
dy = y - self.rulersStartY
lenruler = (dx**2 + dy**2)**0.5
lenruler *= self.scale
self.rulersLen[-1].setPos(txtPosX, txtPosY)
if lenruler > 1:
txt = '%.3f' % lenruler
else:
txt = '%s' % lenruler
if self.pAngle.value():
txt += '; angle=%.2f DEG' % np.degrees(np.arctan2(-dy, dx))
self.rulersLen[-1].setText(txt)
self.rulers[-1].setData(
(self.rulersStartX, x),
(self.rulersStartY, y))
def __init__(self, img, batch, usage, ID, p, a=0, v=60, l=4, ):
pyglet.sprite.Sprite.__init__(self, img=img, batch=batch, usage=usage)
self.a = a
self.v = v/3.6 # convert to m/s
self.p = p
self.l = l # length
self.ID = ID
self.scale = 0.05
self.image.anchor_x = self.image.width / 2
self.image.anchor_y = self.image.height / 2
self.length = self.image.width
window.pixel_unit = self.l / self.width
self.central_radian = window.unit_to_screen(self.p)/window.centre_radius
dx = window.centre_radius * np.cos(self.central_radian)
dy = window.centre_radius * np.sin(self.central_radian)
self.position = window.region_centre + np.array([dx, dy])
self.rotation = -np.degrees([self.central_radian + np.pi/2])
self.isCollide = False
self.reward = 0
def test_calc_ocb_vec_sign(self):
""" Test the calculation of the OCB vector signs
"""
# Set the initial values
self.vdata.ocb_aacgm_mlt = self.ocb.phi_cent[self.vdata.ocb_ind] / 15.0
self.vdata.ocb_aacgm_lat = 90.0 - self.ocb.r_cent[self.vdata.ocb_ind]
(self.vdata.ocb_lat,
self.vdata.ocb_mlt) = self.ocb.normal_coord(self.vdata.aacgm_lat,
self.vdata.aacgm_mlt)
self.vdata.calc_vec_pole_angle()
self.vdata.define_quadrants()
vmag = np.sqrt(self.vdata.aacgm_n**2 + self.vdata.aacgm_e**2)
self.vdata.aacgm_naz = np.degrees(np.arccos(self.vdata.aacgm_n / vmag))
# Calculate the vector data signs
vsigns = self.vdata.calc_ocb_vec_sign(north=True, east=True)
self.assertTrue(vsigns['north'])
self.assertTrue(vsigns['east'])
def spherical_to_cartesian(s,degrees=True,normalize=False):
'''
Takes a vector in spherical coordinates and converts it to cartesian.
Assumes the input vector is given as [radius,colat,lon]
'''
if degrees:
s[1] = np.radians(s[1])
s[2] = np.radians(s[2])
x1 = s[0]*np.sin(s[1])*np.cos(s[2])
x2 = s[0]*np.sin(s[1])*np.sin(s[2])
x3 = s[0]*np.cos(s[1])
x = [x1,x2,x3]
if normalize:
x /= np.linalg.norm(x)
return x
def rotate_delays(lat_r,lon_r,lon_0=0.0,lat_0=0.0,degrees=0):
'''
Rotates the source and receiver of a trace object around an
arbitrary axis.
'''
alpha = np.radians(degrees)
colat_r = 90.0-lat_r
colat_0 = 90.0-lat_0
x_r = lon_r - lon_0
y_r = colat_0 - colat_r
#rotate receivers
lat_rotated = 90.0-colat_0+x_r*np.sin(alpha) + y_r*np.cos(alpha)
lon_rotated = lon_0+x_r*np.cos(alpha) - y_r*np.sin(alpha)
return lat_rotated, lon_rotated
def test_planets(date):
p = np.degrees(_planets(date))
assert abs(p[0] - 314.9122873) < 1e-7
assert abs(p[1] - 91.9393769) < 1e-7
assert abs(p[2] - 169.0970043) < 1e-7
assert abs(p[3] - 196.7516428) < 1e-7
assert abs(p[4] - 42.6046467) < 1e-7
assert abs(p[5] % 360. - 143.319167) < 1e-6
assert abs(p[6] % 360. - 156.221635) < 1e-6
assert abs(p[7] % 360. - 194.890465) < 1e-6
assert abs(p[8] % 360. - 91.262347) < 1e-6
assert abs(p[9] % 360. - 163.710186) < 1e-6
assert abs(p[10] % 360. - 102.168400) < 1e-5 # <== I don't know why but this one is not precise enought
assert abs(p[11] % 360. - 332.317825) < 1e-6
assert abs(p[12] % 360. - 313.661341) < 1e-6
assert abs(p[13] % 360. - 0.059545) < 1e-6
def error_ellipse(mu, cov, ax=None, factor=1.0, **kwargs):
"""
Plot the error ellipse at a point given its covariance matrix.
"""
# some sane defaults
facecolor = kwargs.pop('facecolor', 'none')
edgecolor = kwargs.pop('edgecolor', 'k')
x, y = mu
U, S, V = np.linalg.svd(cov)
theta = np.degrees(np.arctan2(U[1, 0], U[0, 0]))
ellipsePlot = Ellipse(xy=[x, y],
width=2 * np.sqrt(S[0]) * factor,
height=2 * np.sqrt(S[1]) * factor,
angle=theta,
facecolor=facecolor, edgecolor=edgecolor, **kwargs)
if ax is None:
ax = plt.gca()
ax.add_patch(ellipsePlot)
return ellipsePlot
