Python numpy 模块，radians() 实例源码
我们从Python开源项目中，提取了以下50个代码示例，用于说明如何使用numpy.radians()。
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 great_circle_dist(p1, p2):
"""Return the distance (in km) between two points in
geographical coordinates.
"""
lon0, lat0 = p1
lon1, lat1 = p2
EARTH_R = 6372.8
lat0 = np.radians(float(lat0))
lon0 = np.radians(float(lon0))
lat1 = np.radians(float(lat1))
lon1 = np.radians(float(lon1))
dlon = lon0 - lon1
y = np.sqrt(
(np.cos(lat1) * np.sin(dlon)) ** 2
+ (np.cos(lat0) * np.sin(lat1)
- np.sin(lat0) * np.cos(lat1) * np.cos(dlon)) ** 2)
x = np.sin(lat0) * np.sin(lat1) + \
np.cos(lat0) * np.cos(lat1) * np.cos(dlon)
c = np.arctan2(y, x)
return EARTH_R * c
def get_data(filename,headers,ph_units):
# Importation des données .DAT
dat_file = np.loadtxt("%s"%(filename),skiprows=headers,delimiter=',')
labels = ["freq", "amp", "pha", "amp_err", "pha_err"]
data = {l:dat_file[:,i] for (i,l) in enumerate(labels)}
if ph_units == "mrad":
data["pha"] = data["pha"]/1000 # mrad to rad
data["pha_err"] = data["pha_err"]/1000 # mrad to rad
if ph_units == "deg":
data["pha"] = np.radians(data["pha"]) # deg to rad
data["pha_err"] = np.radians(data["pha_err"]) # deg to rad
data["phase_range"] = abs(max(data["pha"])-min(data["pha"])) # Range of phase measurements (used in NRMS error calculation)
data["Z"] = data["amp"]*(np.cos(data["pha"]) + 1j*np.sin(data["pha"]))
EI = np.sqrt(((data["amp"]*np.cos(data["pha"])*data["pha_err"])**2)+(np.sin(data["pha"])*data["amp_err"])**2)
ER = np.sqrt(((data["amp"]*np.sin(data["pha"])*data["pha_err"])**2)+(np.cos(data["pha"])*data["amp_err"])**2)
data["Z_err"] = ER + 1j*EI
# Normalization of amplitude
data["Z_max"] = max(abs(data["Z"])) # Maximum amplitude
zn, zn_e = data["Z"]/data["Z_max"], data["Z_err"]/data["Z_max"] # Normalization of impedance by max amplitude
data["zn"] = np.array([zn.real, zn.imag]) # 2D array with first column = real values, second column = imag values
data["zn_err"] = np.array([zn_e.real, zn_e.imag]) # 2D array with first column = real values, second column = imag values
return data
def get_data(filename,headers,ph_units):
# Importation des données .DAT
dat_file = np.loadtxt("%s"%(filename),skiprows=headers,delimiter=',')
labels = ["freq", "amp", "pha", "amp_err", "pha_err"]
data = {l:dat_file[:,i] for (i,l) in enumerate(labels)}
if ph_units == "mrad":
data["pha"] = data["pha"]/1000 # mrad to rad
data["pha_err"] = data["pha_err"]/1000 # mrad to rad
if ph_units == "deg":
data["pha"] = np.radians(data["pha"]) # deg to rad
data["pha_err"] = np.radians(data["pha_err"]) # deg to rad
data["phase_range"] = abs(max(data["pha"])-min(data["pha"])) # Range of phase measurements (used in NRMS error calculation)
data["Z"] = data["amp"]*(np.cos(data["pha"]) + 1j*np.sin(data["pha"]))
EI = np.sqrt(((data["amp"]*np.cos(data["pha"])*data["pha_err"])**2)+(np.sin(data["pha"])*data["amp_err"])**2)
ER = np.sqrt(((data["amp"]*np.sin(data["pha"])*data["pha_err"])**2)+(np.cos(data["pha"])*data["amp_err"])**2)
data["Z_err"] = ER + 1j*EI
# Normalization of amplitude
data["Z_max"] = max(abs(data["Z"])) # Maximum amplitude
zn, zn_e = data["Z"]/data["Z_max"], data["Z_err"]/data["Z_max"] # Normalization of impedance by max amplitude
data["zn"] = np.array([zn.real, zn.imag]) # 2D array with first column = real values, second column = imag values
data["zn_err"] = np.array([zn_e.real, zn_e.imag]) # 2D array with first column = real values, second column = imag values
return data
def plotArc(start_angle, stop_angle, radius, width, **kwargs):
""" write a docstring for this function"""
numsegments = 100
theta = np.radians(np.linspace(start_angle+90, stop_angle+90, numsegments))
centerx = 0
centery = 0
x1 = -np.cos(theta) * (radius)
y1 = np.sin(theta) * (radius)
stack1 = np.column_stack([x1, y1])
x2 = -np.cos(theta) * (radius + width)
y2 = np.sin(theta) * (radius + width)
stack2 = np.column_stack([np.flip(x2, axis=0), np.flip(y2,axis=0)])
#add the first values from the first set to close the polygon
np.append(stack2, [[x1[0],y1[0]]], axis=0)
arcArray = np.concatenate((stack1,stack2), axis=0)
return patches.Polygon(arcArray, True, **kwargs), ((x1, y1), (x2, y2))
def orthogonalization_matrix(lengths, angles):
"""Return orthogonalization matrix for crystallographic cell coordinates.
Angles are expected in degrees.
The de-orthogonalization matrix is the inverse.
>>> O = orthogonalization_matrix([10, 10, 10], [90, 90, 90])
>>> numpy.allclose(O[:3, :3], numpy.identity(3, float) * 10)
True
>>> O = orthogonalization_matrix([9.8, 12.0, 15.5], [87.2, 80.7, 69.7])
>>> numpy.allclose(numpy.sum(O), 43.063229)
True
"""
a, b, c = lengths
angles = numpy.radians(angles)
sina, sinb, _ = numpy.sin(angles)
cosa, cosb, cosg = numpy.cos(angles)
co = (cosa * cosb - cosg) / (sina * sinb)
return numpy.array([
[ a*sinb*math.sqrt(1.0-co*co), 0.0, 0.0, 0.0],
[-a*sinb*co, b*sina, 0.0, 0.0],
[ a*cosb, b*cosa, c, 0.0],
[ 0.0, 0.0, 0.0, 1.0]])
def computeRotMatrix(self,Phi=False):
#######################################
# COMPUTE ROTATION MATRIX SUCH THAT m(t) = A*L(t)*A'*Hp
# Default set such that phi1,phi2 = 0 is UXO pointed towards North
if Phi is False:
phi1 = np.radians(self.phi[0])
phi2 = np.radians(self.phi[1])
phi3 = np.radians(self.phi[2])
else:
phi1 = np.radians(Phi[0]) # Roll (CCW)
phi2 = np.radians(Phi[1]) # Inclination (+ve is nose pointing down)
phi3 = np.radians(Phi[2]) # Declination (degrees CW from North)
# A1 = np.r_[np.c_[np.cos(phi1),-np.sin(phi1),0.],np.c_[np.sin(phi1),np.cos(phi1),0.],np.c_[0.,0.,1.]] # CCW Rotation about z-axis
# A2 = np.r_[np.c_[1.,0.,0.],np.c_[0.,np.cos(phi2),np.sin(phi2)],np.c_[0.,-np.sin(phi2),np.cos(phi2)]] # CW Rotation about x-axis (rotates towards North)
# A3 = np.r_[np.c_[np.cos(phi3),-np.sin(phi3),0.],np.c_[np.sin(phi3),np.cos(phi3),0.],np.c_[0.,0.,1.]] # CCW Rotation about z-axis (direction of head of object)
A1 = np.r_[np.c_[np.cos(phi1),np.sin(phi1),0.],np.c_[-np.sin(phi1),np.cos(phi1),0.],np.c_[0.,0.,1.]] # CW Rotation about z-axis
A2 = np.r_[np.c_[1.,0.,0.],np.c_[0.,np.cos(phi2),np.sin(phi2)],np.c_[0.,-np.sin(phi2),np.cos(phi2)]] # CW Rotation about x-axis (rotates towards North)
A3 = np.r_[np.c_[np.cos(phi3),np.sin(phi3),0.],np.c_[-np.sin(phi3),np.cos(phi3),0.],np.c_[0.,0.,1.]] # CW Rotation about z-axis (direction of head of object)
return np.dot(A3,np.dot(A2,A1))
def orthogonalization_matrix(lengths, angles):
"""Return orthogonalization matrix for crystallographic cell coordinates.
Angles are expected in degrees.
The de-orthogonalization matrix is the inverse.
>>> O = orthogonalization_matrix([10, 10, 10], [90, 90, 90])
>>> numpy.allclose(O[:3, :3], numpy.identity(3, float) * 10)
True
>>> O = orthogonalization_matrix([9.8, 12.0, 15.5], [87.2, 80.7, 69.7])
>>> numpy.allclose(numpy.sum(O), 43.063229)
True
"""
a, b, c = lengths
angles = numpy.radians(angles)
sina, sinb, _ = numpy.sin(angles)
cosa, cosb, cosg = numpy.cos(angles)
co = (cosa * cosb - cosg) / (sina * sinb)
return numpy.array([
[ a*sinb*math.sqrt(1.0-co*co), 0.0, 0.0, 0.0],
[-a*sinb*co, b*sina, 0.0, 0.0],
[ a*cosb, b*cosa, c, 0.0],
[ 0.0, 0.0, 0.0, 1.0]])
def test_projection(self):
"""Tests the electric field projection."""
projection = self.field.projection
# Top-right quadrant
a = radians(45)
self.assertTrue(isclose(projection([0, 0], a), -4*cos(a)))
self.assertTrue(isclose(projection([3, 0], a), 0.375*cos(a)))
self.assertTrue(isclose(projection([0, 1], a), -sqrt(2)*cos(a)))
self.assertTrue(isclose(projection([[0, 0], [3, 0], [0, 1]], a),
array([-4, 0.375, -sqrt(2)])*cos(a)).all())
# Bottom-left quadrant
a1 = radians(-135)
a2 = radians(45)
self.assertTrue(isclose(projection([0, 0], a1), 4*cos(a2)))
self.assertTrue(isclose(projection([3, 0], a1), -0.375*cos(a2)))
self.assertTrue(isclose(projection([0, 1], a1),
sqrt(2)*cos(a2)))
self.assertTrue(isclose(projection([[0, 0], [3, 0], [0, 1]], a1),
array([4, -0.375, sqrt(2)])*cos(a2)).all())
def cie_relative_luminance(sky_elevation, sky_azimuth=None, sun_elevation=None,
sun_azimuth=None, type='soc'):
""" cie relative luminance of a sky element relative to the luminance
at zenith
angle in radians
type is one of 'soc' (standard overcast sky), 'uoc' (uniform radiance)
or 'clear_sky' (standard clear sky low turbidity)
"""
if type == 'clear_sky' and (
sun_elevation is None or sun_azimuth is None or sky_azimuth is None):
raise ValueError, 'Clear sky requires sun position'
if type == 'soc':
return cie_luminance_gradation(sky_elevation, 4, -0.7)
elif type == 'uoc':
return cie_luminance_gradation(sky_elevation, 0, -1)
elif type == 'clear_sky':
return cie_luminance_gradation(sky_elevation, -1,
-0.32) * cie_scattering_indicatrix(
sun_azimuth, sun_elevation, sky_azimuth, sky_elevation, 10, -3,
0.45)
else:
raise ValueError, 'Unknown sky type'
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 actual_sky_irradiances(dates, ghi, dhi=None, Tdew=None, longitude=_longitude, latitude=_latitude, altitude=_altitude, method='dirint'):
""" derive a sky irradiance dataframe from actual weather data"""
df = sun_position(dates=dates, latitude=latitude, longitude=longitude, altitude=altitude, filter_night=False)
df['am'] = air_mass(df['zenith'], altitude)
df['dni_extra'] = sun_extraradiation(df.index)
if dhi is None:
pressure = pvlib.atmosphere.alt2pres(altitude)
dhi = diffuse_horizontal_irradiance(ghi, df['elevation'], dates, pressure=pressure, temp_dew=Tdew, method=method)
df['ghi'] = ghi
df['dhi'] = dhi
el = numpy.radians(df['elevation'])
df['dni'] = (df['ghi'] - df['dhi']) / numpy.sin(el)
df['brightness'] = brightness(df['am'], df['dhi'], df['dni_extra'])
df['clearness'] = clearness(df['dni'], df['dhi'], df['zenith'])
return df.loc[(df['elevation'] > 0) & (df['ghi'] > 0) , ['azimuth', 'zenith', 'elevation', 'am', 'dni_extra', 'clearness', 'brightness', 'ghi', 'dni', 'dhi' ]]
def rotation_matrix(alpha, beta, gamma):
"""
Return the rotation matrix used to rotate a set of cartesian
coordinates by alpha radians about the z-axis, then beta radians
about the y'-axis and then gamma radians about the z''-axis.
"""
ALPHA = np.array([[np.cos(alpha), -np.sin(alpha), 0],
[np.sin(alpha), np.cos(alpha), 0],
[0, 0, 1]])
BETA = np.array([[np.cos(beta), 0, np.sin(beta)],
[0, 1, 0],
[-np.sin(beta), 0, np.cos(beta)]])
GAMMA = np.array([[np.cos(gamma), -np.sin(gamma), 0],
[np.sin(gamma), np.cos(gamma), 0],
[0, 0, 1]])
R = ALPHA.dot(BETA).dot(GAMMA)
return(R)
def __init__(self, lat0, lon0, depth0, nlat, nlon, ndepth, dlat, dlon, ddepth):
# NOTE: Origin of spherical coordinate system and geographic coordinate
# system is not the same!
# Geographic coordinate system
self.lat0, self.lon0, self.depth0 =\
seispy.coords.as_geographic([lat0, lon0, depth0])
self.nlat, self.nlon, self.ndepth = nlat, nlon, ndepth
self.dlat, self.dlon, self.ddepth = dlat, dlon, ddepth
# Spherical/Pseudospherical coordinate systems
self.nrho = self.ndepth
self.ntheta = self.nlambda = self.nlat
self.nphi = self.nlon
self.drho = self.ddepth
self.dtheta = self.dlambda = np.radians(self.dlat)
self.dphi = np.radians(self.dlon)
self.rho0 = seispy.constants.EARTH_RADIUS\
- (self.depth0 + (self.ndepth - 1) * self.ddepth)
self.lambda0 = np.radians(self.lat0)
self.theta0 = ?/2 - (self.lambda0 + (self.nlambda - 1) * self.dlambda)
self.phi0 = np.radians(self.lon0)
def _add_triangular_sides(self, xy_mask, angle, y_top_right, y_bot_left,
x_top_right, x_bot_left, n_material):
angle = np.radians(angle)
trap_len = (y_top_right - y_bot_left) / np.tan(angle)
num_x_iterations = round(trap_len/self.x_step)
y_per_iteration = round(self.y_pts / num_x_iterations)
lhs_x_start_index = int(x_bot_left/ self.x_step + 0.5)
rhs_x_stop_index = int(x_top_right/ self.x_step + 1 + 0.5)
for i, _ in enumerate(xy_mask):
xy_mask[i][:lhs_x_start_index] = False
xy_mask[i][lhs_x_start_index:rhs_x_stop_index] = True
if i % y_per_iteration == 0:
lhs_x_start_index -= 1
rhs_x_stop_index += 1
self.n[xy_mask] = n_material
return self.n
def orthogonalization_matrix(lengths, angles):
"""Return orthogonalization matrix for crystallographic cell coordinates.
Angles are expected in degrees.
The de-orthogonalization matrix is the inverse.
>>> O = orthogonalization_matrix([10, 10, 10], [90, 90, 90])
>>> numpy.allclose(O[:3, :3], numpy.identity(3, float) * 10)
True
>>> O = orthogonalization_matrix([9.8, 12.0, 15.5], [87.2, 80.7, 69.7])
>>> numpy.allclose(numpy.sum(O), 43.063229)
True
"""
a, b, c = lengths
angles = numpy.radians(angles)
sina, sinb, _ = numpy.sin(angles)
cosa, cosb, cosg = numpy.cos(angles)
co = (cosa * cosb - cosg) / (sina * sinb)
return numpy.array([
[ a*sinb*math.sqrt(1.0-co*co), 0.0, 0.0, 0.0],
[-a*sinb*co, b*sina, 0.0, 0.0],
[ a*cosb, b*cosa, c, 0.0],
[ 0.0, 0.0, 0.0, 1.0]])
def semiMajAmp(m1,m2,inc,ecc,p):
"""
K = [(2*pi*G)/p]^(1/3) [m2*sin(i)/m2^(2/3)*sqrt(1-e^2)]
units:
K [m/s]
m1 [Msun]
m2 [Mj]
p [yrs]
inc [deg]
"""
pSecs = p*sec_per_year
m1KG = m1*const.M_sun.value
A = ((2.0*np.pi*const.G.value)/pSecs)**(1.0/3.0)
B = (m2*const.M_jup.value*np.sin(np.radians(inc)))
C = m1KG**(2.0/3.0)*np.sqrt(1.0-ecc**2.0)
print('Resulting K is '+repr(A*(B/C)))
#return A*(B/C)
def inc_prior_fn(self, inc):
ret = 1.0
if (self.inc_prior==False) or (self.inc_prior=="uniform"):
ret = self.uniform_fn(inc,7)
else:
if inc not in [0.0,90.0,180.0]:
mn = self.mins_ary[7]
mx = self.maxs_ary[7]
# Account for case where min = -(max), as it causes error in denominator
if mn == -1*mx:
mn = mn-0.1
mn_rad = np.radians(mn)
mx_rad = np.radians(mx)
inc_rad = np.radians(inc)
if (self.inc_prior == True) or (self.inc_prior == 'sin'):
ret = np.abs(np.sin(inc_rad)) / np.abs(np.cos(mn_rad)-np.cos(mx_rad))
elif self.inc_prior == 'cos':
ret = np.abs(np.cos(inc_rad)) / np.abs(np.cos(mn_rad)-np.cos(mx_rad))
#if ret==0: ret=-np.inf
return ret
def inc_prior_fn(self, inc):
ret = 1.0
if (self.inc_prior==False) or (self.inc_prior=="uniform"):
ret = self.uniform_fn(inc,7)
else:
if inc not in [0.0,90.0,180.0]:
mn = self.mins_ary[7]
mx = self.maxs_ary[7]
# Account for case where min = -(max), as it causes error in denominator
if mn == -1*mx:
mn = mn-0.1
mn_rad = np.radians(mn)
mx_rad = np.radians(mx)
inc_rad = np.radians(inc)
if (self.inc_prior == True) or (self.inc_prior == 'sin'):
ret = np.abs(np.sin(inc_rad)) / np.abs(np.cos(mn_rad)-np.cos(mx_rad))
elif self.inc_prior == 'cos':
ret = np.abs(np.cos(inc_rad)) / np.abs(np.cos(mn_rad)-np.cos(mx_rad))
#if ret==0: ret=-np.inf
return ret
def hav_dist(locs1, locs2):
"""
Return a distance matrix between two set of coordinates.
Use geometric distance (default) or haversine distance (if longlat=True).
Parameters
----------
locs1 : numpy.array
The first set of coordinates as [(long, lat), (long, lat)].
locs2 : numpy.array
The second set of coordinates as [(long, lat), (long, lat)].
Returns
-------
mat_dist : numpy.array
The distance matrix between locs1 and locs2
"""
locs1 = np.radians(locs1)
locs2 = np.radians(locs2)
cos_lat1 = np.cos(locs1[..., 0])
cos_lat2 = np.cos(locs2[..., 0])
cos_lat_d = np.cos(locs1[..., 0] - locs2[..., 0])
cos_lon_d = np.cos(locs1[..., 1] - locs2[..., 1])
return 6367000 * np.arccos(
cos_lat_d - cos_lat1 * cos_lat2 * (1 - cos_lon_d))
def _get_rotation_matrix(axis, angle):
"""
Helper function to generate a rotation matrix for an x, y, or z axis
Used in get_major_angles
"""
cos = np.cos
sin = np.sin
angle = np.radians(angle)
if axis == 2:
# z axis
return np.array([[cos(angle), -sin(angle), 0], [sin(angle), cos(angle), 0], [0, 0, 1]])
if axis == 1:
# y axis
return np.array([[cos(angle), 0, sin(angle)], [0, 1, 0], [-sin(angle), 0, cos(angle)]])
else:
# x axis
return np.array([[1, 0, 0], [0, cos(angle), -sin(angle)], [0, sin(angle), cos(angle)]])
def euler(xyz, order='xyz', units='deg'):
if not hasattr(xyz, '__iter__'):
xyz = [xyz]
if units == 'deg':
xyz = np.radians(xyz)
r = np.eye(3)
for theta, axis in zip(xyz, order):
c = np.cos(theta)
s = np.sin(theta)
if axis == 'x':
r = np.dot(np.array([[1, 0, 0], [0, c, -s], [0, s, c]]), r)
if axis == 'y':
r = np.dot(np.array([[c, 0, s], [0, 1, 0], [-s, 0, c]]), r)
if axis == 'z':
r = np.dot(np.array([[c, -s, 0], [s, c, 0], [0, 0, 1]]), r)
return r
def get_q_per_pixel(self):
'''Gets the delta-q associated with a single pixel. This is computed in
the small-angle limit, so it should only be considered approximate.
For instance, wide-angle detectors will have different delta-q across
the detector face.'''
if self.q_per_pixel is not None:
return self.q_per_pixel
c = (self.pixel_size_um/1e6)/self.distance_m
twotheta = np.arctan(c) # radians
self.q_per_pixel = 2.0*self.get_k()*np.sin(twotheta/2.0)
return self.q_per_pixel
# Maps
########################################
def add_motion_spikes(motion_mat,frequency,severity,TR):
time = motion_mat[:,0]
max_translation = 5 * severity / 1000 * np.sqrt(2*np.pi)#Max translations in m, factor of sqrt(2*pi) accounts for normalisation factor in norm.pdf later on
max_rotation = np.radians(5 * severity) *np.sqrt(2*np.pi) #Max rotation in rad
time_blocks = np.floor(time[-1]/TR).astype(np.int32)
for i in range(time_blocks):
if np.random.uniform(0,1) < frequency: #Decide whether to add spike
for j in range(1,4):
if np.random.uniform(0,1) < 1/6:
motion_mat[:,j] = motion_mat[:,j] \
+ max_translation * random.uniform(-1,1) \
* norm.pdf(time,loc = (i+0.5)*TR,scale = TR/5)
for j in range(4,7):
if np.random.uniform(0,1) < 1/6:
motion_mat[:,j] = motion_mat[:,j] \
+ max_rotation * random.uniform(-1,1) \
* norm.pdf(time,loc = (i+0.5 + np.random.uniform(-0.25,-.25))*TR,scale = TR/5)
return motion_mat
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 setPose(self, obj_center=None, distance=None,
rotation=None, tilt=None, pitch=None):
tvec, rvec = self.pose()
if distance is not None:
tvec[2, 0] = distance
if obj_center is not None:
tvec[0, 0] = obj_center[0]
tvec[1, 0] = obj_center[1]
if rotation is None and tilt is None:
return rvec
r, t, p = rvec2euler(rvec)
if rotation is not None:
r = np.radians(rotation)
if tilt is not None:
t = np.radians(tilt)
if pitch is not None:
p = np.radians(pitch)
rvec = euler2rvec(r, t, p)
self._pose = tvec, rvec
def hav(alpha):
""" Formula for haversine
Parameters
----------
alpha : (float)
Angle in radians
Returns
--------
hav_alpha : (float)
Haversine of alpha, equal to the square of the sine of half-alpha
"""
hav_alpha = np.sin(alpha * 0.5)**2
return hav_alpha
def archav(hav):
""" Formula for the inverse haversine
Parameters
-----------
hav : (float)
Haversine of an angle
Returns
---------
alpha : (float)
Angle in radians
"""
alpha = 2.0 * np.arcsin(np.sqrt(hav))
return alpha
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 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 rotate_setup(self,lon_0=60.0,colat_0=90.0,degrees=0):
alpha = np.radians(degrees)
lon_s = self.sy
lon_r = self.ry
colat_s = self.sx
colat_r = self.rx
x_s = lon_s - lon_0
y_s = colat_0 - colat_s
x_r = lon_r - lon_0
y_r = colat_0 - colat_r
#rotate receiver
self.rx = colat_0+x_r*np.sin(alpha) + y_r*np.cos(alpha)
self.ry = lon_0+x_r*np.cos(alpha) - y_r*np.sin(alpha)
#rotate source
self.sx = colat_0+x_s*np.sin(alpha) + y_s*np.cos(alpha)
self.sy = lon_0+x_s*np.cos(alpha) - y_s*np.sin(alpha)
#########################################################################
# Plot map of earthquake and station
#########################################################################
def _calcArea_(self, v1, v2):
"""
Private method to calculate the area covered by a spherical
quadrilateral with one corner defined by the normal vectors
of the two intersecting great circles.
INPUTS:
v1, v2: float array(3), the normal vectors
RETURNS:
area: float, the area given in square radians
"""
angle = self.calcAngle(v1, v2)
area = (4*angle - 2*np.math.pi)
return area
def _rotVector_(self, v, angle, axis):
"""
Rotate a vector by an angle around an axis
INPUTS:
v: 3-dim float array
angle: float, the rotation angle in radians
axis: string, 'x', 'y', or 'z'
RETURNS:
float array(3): the rotated vector
"""
# axisd = {'x':[1,0,0], 'y':[0,1,0], 'z':[0,0,1]}
# construct quaternion and rotate...
rot = cgt.quat()
rot.fromAngleAxis(angle, axis)
return list(rot.rotateVec(v))
def _geodetic_to_cartesian(lat, lon, alt):
"""Conversion from latitude, longitue and altitude coordinates to
cartesian with respect to an ellipsoid
Args:
lat (float): Latitude in radians
lon (float): Longitue in radians
alt (float): Altitude to sea level in meters
Return:
numpy.array: 3D element (in meters)
"""
C = Earth.r / np.sqrt(1 - (Earth.e * np.sin(lat)) ** 2)
S = Earth.r * (1 - Earth.e ** 2) / np.sqrt(1 - (Earth.e * np.sin(lat)) ** 2)
r_d = (C + alt) * np.cos(lat)
r_k = (S + alt) * np.sin(lat)
norm = np.sqrt(r_d ** 2 + r_k ** 2)
return norm * np.array([
np.cos(lat) * np.cos(lon),
np.cos(lat) * np.sin(lon),
np.sin(lat)
])
def orthogonalization_matrix(lengths, angles):
"""Return orthogonalization matrix for crystallographic cell coordinates.
Angles are expected in degrees.
The de-orthogonalization matrix is the inverse.
>>> O = orthogonalization_matrix([10, 10, 10], [90, 90, 90])
>>> numpy.allclose(O[:3, :3], numpy.identity(3, float) * 10)
True
>>> O = orthogonalization_matrix([9.8, 12.0, 15.5], [87.2, 80.7, 69.7])
>>> numpy.allclose(numpy.sum(O), 43.063229)
True
"""
a, b, c = lengths
angles = numpy.radians(angles)
sina, sinb, _ = numpy.sin(angles)
cosa, cosb, cosg = numpy.cos(angles)
co = (cosa * cosb - cosg) / (sina * sinb)
return numpy.array([
[ a*sinb*math.sqrt(1.0-co*co), 0.0, 0.0, 0.0],
[-a*sinb*co, b*sina, 0.0, 0.0],
[ a*cosb, b*cosa, c, 0.0],
[ 0.0, 0.0, 0.0, 1.0]])
def orthogonalization_matrix(lengths, angles):
"""Return orthogonalization matrix for crystallographic cell coordinates.
Angles are expected in degrees.
The de-orthogonalization matrix is the inverse.
>>> O = orthogonalization_matrix([10, 10, 10], [90, 90, 90])
>>> numpy.allclose(O[:3, :3], numpy.identity(3, float) * 10)
True
>>> O = orthogonalization_matrix([9.8, 12.0, 15.5], [87.2, 80.7, 69.7])
>>> numpy.allclose(numpy.sum(O), 43.063229)
True
"""
a, b, c = lengths
angles = numpy.radians(angles)
sina, sinb, _ = numpy.sin(angles)
cosa, cosb, cosg = numpy.cos(angles)
co = (cosa * cosb - cosg) / (sina * sinb)
return numpy.array([
[ a*sinb*math.sqrt(1.0-co*co), 0.0, 0.0, 0.0],
[-a*sinb*co, b*sina, 0.0, 0.0],
[ a*cosb, b*cosa, c, 0.0],
[ 0.0, 0.0, 0.0, 1.0]])
def orthogonalization_matrix(lengths, angles):
"""Return orthogonalization matrix for crystallographic cell coordinates.
Angles are expected in degrees.
The de-orthogonalization matrix is the inverse.
>>> O = orthogonalization_matrix([10, 10, 10], [90, 90, 90])
>>> numpy.allclose(O[:3, :3], numpy.identity(3, float) * 10)
True
>>> O = orthogonalization_matrix([9.8, 12.0, 15.5], [87.2, 80.7, 69.7])
>>> numpy.allclose(numpy.sum(O), 43.063229)
True
"""
a, b, c = lengths
angles = numpy.radians(angles)
sina, sinb, _ = numpy.sin(angles)
cosa, cosb, cosg = numpy.cos(angles)
co = (cosa * cosb - cosg) / (sina * sinb)
return numpy.array([
[ a*sinb*math.sqrt(1.0-co*co), 0.0, 0.0, 0.0],
[-a*sinb*co, b*sina, 0.0, 0.0],
[ a*cosb, b*cosa, c, 0.0],
[ 0.0, 0.0, 0.0, 1.0]])
def set_radmask(self, cluster, mpcscale):
"""
Assign mask (0/1) values to maskgals for a given cluster
parameters
----------
cluster: Cluster object
mpcscale: float
scaling to go from mpc to degrees (check units) at cluster redshift
results
-------
sets maskgals.mark
"""
# note this probably can be in the superclass, no?
ras = cluster.ra + self.maskgals.x/(mpcscale*SEC_PER_DEG)/np.cos(np.radians(cluster.dec))
decs = cluster.dec + self.maskgals.y/(mpcscale*SEC_PER_DEG)
self.maskgals.mark = self.compute_radmask(ras,decs)
def _calc_bkg_density(self, r, chisq, refmag):
"""
Internal method to compute background filter
parameters
----------
bkg: Background object
background
cosmo: Cosmology object
cosmology scaling info
returns
-------
bcounts: float array
b(x) for the cluster
"""
mpc_scale = np.radians(1.) * self.cosmo.Dl(0, self.z) / (1 + self.z)**2
sigma_g = self.bkg.sigma_g_lookup(self.z, chisq, refmag)
return 2 * np.pi * r * (sigma_g/mpc_scale**2)
def orthogonalization_matrix(lengths, angles):
"""Return orthogonalization matrix for crystallographic cell coordinates.
Angles are expected in degrees.
The de-orthogonalization matrix is the inverse.
>>> O = orthogonalization_matrix([10, 10, 10], [90, 90, 90])
>>> numpy.allclose(O[:3, :3], numpy.identity(3, float) * 10)
True
>>> O = orthogonalization_matrix([9.8, 12.0, 15.5], [87.2, 80.7, 69.7])
>>> numpy.allclose(numpy.sum(O), 43.063229)
True
"""
a, b, c = lengths
angles = numpy.radians(angles)
sina, sinb, _ = numpy.sin(angles)
cosa, cosb, cosg = numpy.cos(angles)
co = (cosa * cosb - cosg) / (sina * sinb)
return numpy.array([
[ a*sinb*math.sqrt(1.0-co*co), 0.0, 0.0, 0.0],
[-a*sinb*co, b*sina, 0.0, 0.0],
[ a*cosb, b*cosa, c, 0.0],
[ 0.0, 0.0, 0.0, 1.0]])
def apply_grad_cartesian_tensor(grad_X, zmat_dist):
"""Apply the gradient for transformation to cartesian space onto zmat_dist.
Args:
grad_X (:class:`numpy.ndarray`): A ``(3, n, n, 3)`` array.
The mathematical details of the index layout is explained in
:meth:`~chemcoord.Cartesian.get_grad_zmat()`.
zmat_dist (:class:`~chemcoord.Zmat`):
Distortions in Zmatrix space.
Returns:
:class:`~chemcoord.Cartesian`: Distortions in cartesian space.
"""
columns = ['bond', 'angle', 'dihedral']
C_dist = zmat_dist.loc[:, columns].values.T
try:
C_dist = C_dist.astype('f8')
C_dist[[1, 2], :] = np.radians(C_dist[[1, 2], :])
except (TypeError, AttributeError):
C_dist[[1, 2], :] = sympy.rad(C_dist[[1, 2], :])
cart_dist = np.tensordot(grad_X, C_dist, axes=([3, 2], [0, 1])).T
from chemcoord.cartesian_coordinates.cartesian_class_main import Cartesian
return Cartesian(atoms=zmat_dist['atom'],
coords=cart_dist, index=zmat_dist.index)
def find_lines(img, acc_threshold=0.25, should_erode=True):
if len(img.shape) == 3 and img.shape[2] == 3: # if it's color
img = cv2.cvtColor(img, cv2.COLOR_RGB2GRAY)
img = cv2.GaussianBlur(img, (11, 11), 0)
img = cv2.adaptiveThreshold(
img,
255,
cv2.ADAPTIVE_THRESH_MEAN_C,
cv2.THRESH_BINARY,
5,
2)
img = cv2.bitwise_not(img)
# thresh = 127
# edges = cv2.threshold(img, thresh, 255, cv2.THRESH_BINARY)[1]
# edges = cv2.Canny(blur, 500, 500, apertureSize=3)
if should_erode:
element = cv2.getStructuringElement(cv2.MORPH_RECT, (4, 4))
img = cv2.erode(img, element)
theta = np.pi/2000
angle_threshold = 2
horizontal = cv2.HoughLines(
img,
1,
theta,
int(acc_threshold * img.shape[1]),
min_theta=np.radians(90 - angle_threshold),
max_theta=np.radians(90 + angle_threshold))
vertical = cv2.HoughLines(
img,
1,
theta,
int(acc_threshold * img.shape[0]),
min_theta=np.radians(-angle_threshold),
max_theta=np.radians(angle_threshold),
)
horizontal = list(horizontal) if horizontal is not None else []
vertical = list(vertical) if vertical is not None else []
horizontal = [line[0] for line in horizontal]
vertical = [line[0] for line in vertical]
horizontal = np.asarray(horizontal)
vertical = np.asarray(vertical)
return horizontal, vertical
def drawSkymapCatalog(ax,lon,lat,**kwargs):
mapping = {
'ait':'aitoff',
'mol':'mollweide',
'lam':'lambert',
'ham':'hammer'
}
kwargs.setdefault('proj','aitoff')
kwargs.setdefault('s',2)
kwargs.setdefault('marker','.')
kwargs.setdefault('c','k')
proj = kwargs.pop('proj')
projection = mapping.get(proj,proj)
#ax.grid()
# Convert from
# [0. < lon < 360] -> [-pi < lon < pi]
# [-90 < lat < 90] -> [-pi/2 < lat < pi/2]
lon,lat= np.radians([lon-360.*(lon>180),lat])
ax.scatter(lon,lat,**kwargs)
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
############################################################
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)
def ang2const(lon,lat,coord='gal'):
import ephem
scalar = np.isscalar(lon)
lon = np.array(lon,copy=False,ndmin=1)
lat = np.array(lat,copy=False,ndmin=1)
if coord.lower() == 'cel':
ra,dec = lon,lat
elif coord.lower() == 'gal':
ra,dec = gal2cel(lon,lat)
else:
msg = "Unrecognized coordinate"
raise Exception(msg)
x,y = np.radians([ra,dec])
const = [ephem.constellation(coord) for coord in zip(x,y)]
if scalar: return const[0]
return const
def orthogonalization_matrix(lengths, angles):
"""Return orthogonalization matrix for crystallographic cell coordinates.
Angles are expected in degrees.
The de-orthogonalization matrix is the inverse.
>>> O = orthogonalization_matrix([10, 10, 10], [90, 90, 90])
>>> numpy.allclose(O[:3, :3], numpy.identity(3, float) * 10)
True
>>> O = orthogonalization_matrix([9.8, 12.0, 15.5], [87.2, 80.7, 69.7])
>>> numpy.allclose(numpy.sum(O), 43.063229)
True
"""
a, b, c = lengths
angles = numpy.radians(angles)
sina, sinb, _ = numpy.sin(angles)
cosa, cosb, cosg = numpy.cos(angles)
co = (cosa * cosb - cosg) / (sina * sinb)
return numpy.array([
[ a*sinb*math.sqrt(1.0-co*co), 0.0, 0.0, 0.0],
[-a*sinb*co, b*sina, 0.0, 0.0],
[ a*cosb, b*cosa, c, 0.0],
[ 0.0, 0.0, 0.0, 1.0]])
def __init__(self, survey_features=None, condition_features=None, time_lag=10.,
filtername='r', twi_change=-18.):
"""
Paramters
---------
time_lag : float (10.)
If there is a gap between observations longer than this, let the filter change (minutes)
twi_change : float (-18.)
The sun altitude to consider twilight starting/ending
"""
self.time_lag = time_lag/60./24. # Convert to days
self.twi_change = np.radians(twi_change)
self.filtername = filtername
if condition_features is None:
self.condition_features = {}
self.condition_features['Current_filter'] = features.Current_filter()
self.condition_features['Sun_moon_alts'] = features.Sun_moon_alts()
self.condition_features['Current_mjd'] = features.Current_mjd()
if survey_features is None:
self.survey_features = {}
self.survey_features['Last_observation'] = features.Last_observation()
super(Strict_filter_basis_function, self).__init__(survey_features=self.survey_features,
condition_features=self.condition_features)
def treexyz(ra, dec):
"""
Utility to convert RA,dec postions in x,y,z space, useful for constructing KD-trees.
Parameters
----------
ra : float or array
RA in radians
dec : float or array
Dec in radians
Returns
-------
x,y,z : floats or arrays
The position of the given points on the unit sphere.
"""
# Note ra/dec can be arrays.
x = np.cos(dec) * np.cos(ra)
y = np.cos(dec) * np.sin(ra)
z = np.sin(dec)
return x, y, z