The following are code examples for showing how to use . They are extracted from open source Python projects. You can vote up the examples you like or vote down the exmaples you don’t like. You can also save this page to your account.
Example 1
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 Example 2
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 Example 3
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) Example 4
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) Example 5
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) Example 6
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 Example 7
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 ##
############################ Example 8
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)) Example 9
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 Example 10
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)) Example 11
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 Example 12
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 Example 13
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 Example 14
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 Example 15
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 Example 16
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 Example 17
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() Example 18
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)
########################################################################## Example 19
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 Example 20
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() Example 21
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))
# ----------------------------------------------------------------- Example 22
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))
# ----------------------------------------------------------------- Example 23
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 Example 24
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)) Example 25
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 Example 26
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']) Example 27
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 Example 28
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
Example 29
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 Example 30
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 Example 31
def setUp(self):
NX = 2
nx = np.linspace(-NX + 0.5, NX - 0.5, num=2 * NX, endpoint=True)
vx = np.linspace(-NX, NX, num=2 * NX, endpoint=True)
meshx, meshy = np.meshgrid(nx, nx)
self.cartgrid = np.dstack((meshx, meshy))
self.values = np.repeat(vx[:, np.newaxis], 2 * NX, 1)
coord = georef.sweep_centroids(4, 1, NX, 0.)
xx = coord[..., 0]
yy = np.degrees(coord[..., 1])
xxx = xx * np.cos(np.radians(90. - yy))
x = xx * np.sin(np.radians(90. - yy))
y = xxx
self.newgrid = np.dstack((x, y))
self.result = np.array([[0.47140452, 1.41421356],
[0.47140452, 1.41421356],
[-0.47140452, -1.41421356],
[-0.47140452, -1.41421356]]) Example 32
def healpixMap(nside, lon, lat, fill_value=0., nest=False):
"""
Input (lon, lat) in degrees instead of (theta, phi) in radians.
Returns HEALPix map at the desired resolution
"""
lon_median, lat_median = np.median(lon), np.median(lat)
max_angsep = np.max(ugali.utils.projector.angsep(lon, lat, lon_median, lat_median))
pix = angToPix(nside, lon, lat, nest=nest)
if max_angsep < 10:
# More efficient histograming for small regions of sky
m = np.tile(fill_value, healpy.nside2npix(nside))
pix_subset = ugali.utils.healpix.angToDisc(nside, lon_median, lat_median, max_angsep, nest=nest)
bins = np.arange(np.min(pix_subset), np.max(pix_subset) + 1)
m_subset = np.histogram(pix, bins=bins - 0.5)[0].astype(float)
m[bins[0:-1]] = m_subset
else:
m = np.histogram(pix, np.arange(hp.nside2npix(nside) + 1))[0].astype(float)
if fill_value != 0.:
m[m == 0.] = fill_value
return m
############################################################ Example 33
def galToCel(ll, bb):
"""
Converts Galactic (deg) to Celestial J2000 (deg) coordinates
"""
bb = numpy.radians(bb)
sin_bb = numpy.sin(bb)
cos_bb = numpy.cos(bb)
ll = numpy.radians(ll)
ra_gp = numpy.radians(192.85948)
de_gp = numpy.radians(27.12825)
lcp = numpy.radians(122.932)
sin_lcp_ll = numpy.sin(lcp - ll)
cos_lcp_ll = numpy.cos(lcp - ll)
sin_d = (numpy.sin(de_gp) * sin_bb) \
+ (numpy.cos(de_gp) * cos_bb * cos_lcp_ll)
ramragp = numpy.arctan2(cos_bb * sin_lcp_ll,
(numpy.cos(de_gp) * sin_bb) \
- (numpy.sin(de_gp) * cos_bb * cos_lcp_ll))
dec = numpy.arcsin(sin_d)
ra = (ramragp + ra_gp + (2. * numpy.pi)) % (2. * numpy.pi)
return numpy.degrees(ra), numpy.degrees(dec) Example 34
def celToGal(ra, dec):
"""
Converts Celestial J2000 (deg) to Calactic (deg) coordinates
"""
dec = numpy.radians(dec)
sin_dec = numpy.sin(dec)
cos_dec = numpy.cos(dec)
ra = numpy.radians(ra)
ra_gp = numpy.radians(192.85948)
de_gp = numpy.radians(27.12825)
sin_ra_gp = numpy.sin(ra - ra_gp)
cos_ra_gp = numpy.cos(ra - ra_gp)
lcp = numpy.radians(122.932)
sin_b = (numpy.sin(de_gp) * sin_dec) \
+ (numpy.cos(de_gp) * cos_dec * cos_ra_gp)
lcpml = numpy.arctan2(cos_dec * sin_ra_gp,
(numpy.cos(de_gp) * sin_dec) \
- (numpy.sin(de_gp) * cos_dec * cos_ra_gp))
bb = numpy.arcsin(sin_b)
ll = (lcp - lcpml + (2. * numpy.pi)) % (2. * numpy.pi)
return numpy.degrees(ll), numpy.degrees(bb) Example 35
def hms2dec(hms):
"""
Convert longitude from hours,minutes,seconds in string or 3-array
format to decimal degrees.
ADW: This really should be replaced by astropy
"""
DEGREE = 360.
HOUR = 24.
MINUTE = 60.
SECOND = 3600.
if isinstance(hms,basestring):
hour,minute,second = numpy.array(re.split('[hms]',hms))[:3].astype(float)
else:
hour,minute,second = hms.T
decimal = (hour + minute * 1./MINUTE + second * 1./SECOND)*(DEGREE/HOUR)
return decimal Example 36
def dms2dec(dms):
"""
Convert latitude from degrees,minutes,seconds in string or 3-array
format to decimal degrees.
"""
DEGREE = 360.
HOUR = 24.
MINUTE = 60.
SECOND = 3600.
# Be careful here, degree needs to be a float so that negative zero
# can have its signbit set:
# http://docs.scipy.org/doc/numpy-1.7.0/reference/c-api.coremath.html#NPY_NZERO
if isinstance(dms,basestring):
degree,minute,second = numpy.array(re.split('[dms]',hms))[:3].astype(float)
else:
degree,minute,second = dms.T
sign = numpy.copysign(1.0,degree)
decimal = numpy.abs(degree) + minute * 1./MINUTE + second * 1./SECOND
decimal *= sign
return decimal Example 37
def _setup_subpix(self,nside=2**16):
"""
Subpixels for random position generation.
"""
# Only setup once...
if hasattr(self,'subpix'): return
# Simulate over full ROI
self.roi_radius = self.config['coords']['roi_radius']
# Setup background spatial stuff
logger.info("Setup subpixels...")
self.nside_pixel = self.config['coords']['nside_pixel']
self.nside_subpixel = self.nside_pixel * 2**4 # Could be config parameter
epsilon = np.degrees(healpy.max_pixrad(self.nside_pixel)) # Pad roi radius to cover edge healpix
subpix = ugali.utils.healpix.query_disc(self.nside_subpixel,self.roi.vec,self.roi_radius+epsilon)
superpix = ugali.utils.healpix.superpixel(subpix,self.nside_subpixel,self.nside_pixel)
self.subpix = subpix[np.in1d(superpix,self.roi.pixels)] Example 38
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 = pl.gca()
ax.add_patch(ellipsePlot)
return ellipsePlot Example 39
def __init__(self, basis_functions, basis_weights, extra_features=None,
smoothing_kernel=None, reward=1e6, ignore_obs='dummy',
nside=default_nside, min_alt=30., max_alt=85.):
"""
min_alt : float (30.)
The minimum altitude to attempt to chace a pair to (degrees). Default of 30 = airmass of 2.
max_alt : float(85.)
The maximum altitude to attempt to chase a pair to (degrees).
"""
self.min_alt = np.radians(min_alt)
self.max_alt = np.radians(max_alt)
self.nside = nside
self.reward_val = reward
self.reward = -reward
if extra_features is None:
extra_features = {'mjd': features.Current_mjd()}
extra_features['altaz'] = features.AltAzFeature(nside=nside)
super(Scripted_survey, self).__init__(basis_functions=basis_functions,
basis_weights=basis_weights,
extra_features=extra_features,
smoothing_kernel=smoothing_kernel,
ignore_obs=ignore_obs) Example 40
def FindSkeleton(self):
rgb = cv2.cvtColor(self.ImgHSV, cv2.COLOR_HSV2BGR)
angle = 0
count = 0
gray = cv2.cvtColor(cv2.cvtColor(self.ImgHSV,cv2.COLOR_HSV2BGR), cv2.COLOR_BGR2GRAY)
edges = cv2.Canny(gray,50,150,apertureSize = 3)
lines = cv2.HoughLines(edges,1,np.pi/180,110)
#print (lines)
line_count = lines.shape[0]
for x in range(line_count):
for rho,theta in lines[x]:
a = np.cos(theta)
b = np.sin(theta)
#print(theta)
x0 = a*rho
y0 = b*rho
x1 = int(x0 + 1000*(-b))
y1 = int(y0 + 1000*(a))
x2 = int(x0 - 1000*(-b))
y2 = int(y0 - 1000*(a))
crr_angle = np.degrees(b)
if (crr_angle < 5):
#print(crr_angle)
angle = angle + crr_angle
count = count + 1
cv2.line(rgb,(x1,y1),(x2,y2),(0,0,255),2)
angle = angle / count
self.angle = angle
return (angle) Example 41
def dip_direction2strike(azimuth):
"""
Converts a planar measurment of dip direction using the dip-azimuth
convention into a strike using the right-hand-rule.
Parameters
----------
azimuth : number or string
The dip direction of the plane in degrees. This can be either a
numerical azimuth in the 0-360 range or a string representing a quadrant
measurement (e.g. N30W).
Returns
-------
strike : number
The strike of the plane in degrees following the right-hand-rule.
"""
azimuth = parse_azimuth(azimuth)
strike = azimuth - 90
if strike < 0:
strike += 360
return strike Example 42
def strike2dip_direction(strike):
"""
Converts a planar measurement of strike using the right-hand-rule into the
dip direction (i.e. the direction that the plane dips).
Parameters
----------
strike : number or string
The strike direction of the plane in degrees. This can be either a
numerical azimuth in the 0-360 range or a string representing a quadrant
measurement (e.g. N30W).
Returns
-------
azimuth : number
The dip direction of the plane in degrees (0-360).
"""
strike = parse_azimuth(strike)
dip_direction = strike + 90
if dip_direction > 360:
dip_direction -= 360
return dip_direction Example 43
def _rotate(lon, lat, theta, axis='x'):
"""
Rotate "lon", "lat" coords (in _degrees_) about the X-axis by "theta"
degrees. This effectively simulates rotating a physical stereonet.
Returns rotated lon, lat coords in _radians_).
"""
# Convert input to numpy arrays in radians
lon, lat = np.atleast_1d(lon, lat)
lon, lat = map(np.radians, [lon, lat])
theta = np.radians(theta)
# Convert to cartesian coords for the rotation
x, y, z = sph2cart(lon, lat)
lookup = {'x':_rotate_x, 'y':_rotate_y, 'z':_rotate_z}
X, Y, Z = lookup[axis](x, y, z, theta)
# Now convert back to spherical coords (longitude and latitude, ignore R)
lon, lat = cart2sph(X,Y,Z)
return lon, lat # in radians! Example 44
def plunge_bearing2pole(plunge, bearing):
"""
Converts the given `plunge` and `bearing` in degrees to a strike and dip
of the plane whose pole would be parallel to the line specified. (i.e. The
pole to the plane returned would plot at the same point as the specified
plunge and bearing.)
Parameters
----------
plunge : number or sequence of numbers
The plunge of the line(s) in degrees. The plunge is measured in degrees
downward from the end of the feature specified by the bearing.
bearing : number or sequence of numbers
The bearing (azimuth) of the line(s) in degrees.
Returns
-------
strike, dip : arrays
Arrays of strikes and dips in degrees following the right-hand-rule.
"""
plunge, bearing = np.atleast_1d(plunge, bearing)
strike = bearing + 90
dip = 90 - plunge
strike[strike >= 360] -= 360
return strike, dip Example 45
def pole2plunge_bearing(strike, dip):
"""
Converts the given *strike* and *dip* in dgrees of a plane(s) to a plunge
and bearing of its pole.
Parameters
----------
strike : number or sequence of numbers
The strike of the plane(s) in degrees, with dip direction indicated by
the azimuth (e.g. 315 vs. 135) specified following the "right hand
rule".
dip : number or sequence of numbers
The dip of the plane(s) in degrees.
Returns
-------
plunge, bearing : arrays
Arrays of plunges and bearings of the pole to the plane(s) in degrees.
"""
strike, dip = np.atleast_1d(strike, dip)
bearing = strike - 90
plunge = 90 - dip
bearing[bearing < 0] += 360
return plunge, bearing Example 46
def geographic2pole(lon, lat):
"""
Converts a longitude and latitude (from a stereonet) into the strike and dip
of the plane whose pole lies at the given longitude(s) and latitude(s).
Parameters
----------
lon : array-like
A sequence of longitudes (or a single longitude) in radians
lat : array-like
A sequence of latitudes (or a single latitude) in radians
Returns
-------
strike : array
A sequence of strikes in degrees
dip : array
A sequence of dips in degrees
"""
plunge, bearing = geographic2plunge_bearing(lon, lat)
strike = bearing + 90
strike[strike >= 360] -= 360
dip = 90 - plunge
return strike, dip Example 47
def azimuth2rake(strike, dip, azimuth):
"""
Projects an azimuth of a linear feature onto a plane as a rake angle.
Parameters
----------
strike, dip : numbers
The strike and dip of the plane in degrees following the
right-hand-rule.
azimuth : numbers
The azimuth of the linear feature in degrees clockwise from north (i.e.
a 0-360 azimuth).
Returns
-------
rake : number
A rake angle in degrees measured downwards from horizontal. Negative
values correspond to the opposite end of the strike.
"""
plunge, bearing = plane_intersection(strike, dip, azimuth, 90)
rake = project_onto_plane(strike, dip, plunge, bearing)
return rake Example 48
def vector2plunge_bearing(x, y, z):
"""
Converts a vector or series of vectors given as x, y, z in world
coordinates into plunge/bearings.
Parameters
----------
x : number or sequence of numbers
The x-component(s) of the normal vector
y : number or sequence of numbers
The y-component(s) of the normal vector
z : number or sequence of numbers
The z-component(s) of the normal vector
Returns
-------
plunge : array
The plunge of the vector in degrees downward from horizontal.
bearing : array
The bearing of the vector in degrees clockwise from north.
"""
return geographic2plunge_bearing(*xyz2stereonet(x,y,z)) Example 49
def set_azimuth_ticks(self, angles, labels=None, frac=None, **kwargs):
"""
Sets the azimuthal tick locations (Note: tick lines are not currently
drawn or supported.).
Parameters
----------
angles : sequence of numbers
The tick locations in degrees.
labels : sequence of strings
The tick label at each location. Defaults to a formatted version
of the specified angles.
frac : number
The radial location of the tick labels. 1.0 is the along the edge,
1.1 would be outside, and 0.9 would be inside.
**kwargs
Additional parameters are text properties for the labels.
"""
return self._polar.set_thetagrids(angles, labels, frac, **kwargs) Example 50
def pole(self, strike, dip, *args, **kwargs):
"""
Plot points representing poles to planes on the axes. Additional
arguments and keyword arguments are passed on to `ax.plot`.
Parameters
----------
strike, dip : numbers or sequences of numbers
The strike and dip of the plane(s) in degrees. The dip direction is
defined by the strike following the "right-hand rule".
**kwargs
Additional parameters are passed on to `plot`.
Returns
-------
A sequence of Line2D artists representing the point(s) specified by
`strike` and `dip`.
"""
lon, lat = stereonet_math.pole(strike, dip)
args, kwargs = self._point_plot_defaults(args, kwargs)
return self.plot(lon, lat, *args, **kwargs) 