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
python类radians()的实例源码
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
