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